D.机器学习

机器学习概论

假设在一项100人的研究中，我们对根据运动员的鞋码、身高和体重来预测他的跑步时间感兴趣。我们可以用一个简单的线性回归模型

Y = beta0 + beta1 *身高+ beta2 *体重+ beta3 *鞋长

这里y是响应变量(运行时间)，n是观察数(100人)，p是变量/特征/预测数(3个IE身高，体重，鞋码)，X是一个nxp矩阵

该数据集是一个低维数据，其中n >> p，但大多数来自现代生物技术的生物数据集是高维IE n << p。这带来了统计挑战，简单的线性回归不再能帮助我们。

例如,

• 根据基因表达数据识别前列腺癌的危险因素(基因)
• 预测乳腺癌患者的存活几率。
• 识别不同亚型乳腺癌的基因表达模式

在上面列出的3个例子中，n，观察的数量，是30-40个病人，而p，特征的数量，大约是30000个基因。试着在以上三种情况中为结果变量y写一个线性回归公式。

下面列出的是高维数据可能出错的地方——其中一些预测器是有用的，一些则不是——如果我们包含太多的预测器，我们就会过度拟合数据

这就是为什么我们需要机器学习。让我们首先介绍一些基本概念，然后深入到例子和实验环节。

监督式学习—使用数据集X来预测与响应变量y的关联。响应变量可以是连续的，也可以是分类的。例如:预测乳腺癌患者的存活几率。

无监督学习-发现x中的关联或模式，没有响应变量。例如:把相似的基因聚在一起。

训练和测试数据集-通常我们将观察分为测试数据集和训练数据集。在训练数据集上拟合模型，在测试数据集上进行评估。测试集错误率是对模型在未来数据集上性能的估计。

模型选择-对于一个给定的问题，我们通常会考虑许多模型。例如,我们试图确定特定疾病的基因通过基因表达数据集,我们可以有以下模型)模型1 -使用所有30000个基因从数组中构建一个模型b)模型2 -我们只包括基因相关疾病的途径,我们知道的是调节在建立一个模型c)模型3 -包括已知基因中发现文学这种疾病影响强烈建议使用测试集只有在我们最后的模型,看看我们的模型将与新的，看不见的数据。那么，我们如何选择可以在测试数据集中进行测试的最佳模型呢?

交叉验证我们可以用不同的方法找到最好的模型。让我们看看常用的方法，也就是，验证集，省略一个交叉验证，k倍交叉验证。

简单来讲,验证集方法将全部数据集分成3组——训练集、验证集和测试集。我们在训练集上训练模型，在验证集上评估它们的性能，然后选择最佳模型来拟合测试集。

的省略一个交叉验证首先拟合n个模型(其中n为训练数据集中的观察数)，每个模型在n-1个观察上，在遗漏的观察上评估每个模型。最佳模型是总测试误差最小的模型，然后用它来预测测试集。

最后,5倍交叉验证(这里k=5)，将训练数据集分成5个集合，在另外4个集合上重复训练模型，在第5个集合上评估性能。

偏差，方差，过拟合-偏差是指实际贝塔值和预测贝塔值之间的平均差异，方差是指实验中贝塔值的差异量。随着模型复杂性(无变量)的增加，偏差减小，方差增大。这就是所谓的偏差-方差权衡，方差过大的模型被称为过拟合。

数据集

为无监督学习，我们将使用bioconductor包中的RNA-Seq计数数据，气道．摘要中，对气道平滑肌(ASM)细胞系的RNA-Seq实验进行了简要描述:“使用RNA-Seq(一种高通量测序方法)，我们对四种原代人ASM细胞系进行了地塞米松(一种高效合成糖皮质激素(1微磨牙18小时)处理后的转录组变化进行了表征。”

库(气道)数据("气道")se <-气道colData(se)

## 8行9列的DataFrame ## SampleName cell dex albut运行avgLength ##       ## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 ## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 ## SRR1039513 GSM1275866 N052611 untrt untrt SRR1039513 87 ## SRR1039516 GSM1275867 N052611 untrt untrt SRR1039516 120 ## srr1275871 N080611 trt untrt SRR1039517 126 ## SRR1039517SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 ## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521生物样品## <因子> <因子> <因子> ## SRR1039508 SRX384345 SRS508568 SAMN02422675 ## SRR1039509 SRX384346 SRS508567 SAMN02422675 ## SRR1039513 SRX384349 SRS508571 SAMN02422678 ## SRR1039513 SRX384353 SRS508575 SAMN02422682 ## SRR1039517 SRX384354 SRS508576 SAMN02422673 ## SRR1039521 SRX384354 SRS508576 SRS508575 ## SRR1039521 srx1039522 SRX384353 SRS508575 SAMN02422670 ## SRR1039517 SRX384354 SRS508576 srs508557 SRS508579Samn02422683 ## srr1039521 srx384358 srs508580 samn02422677

library("DESeq2") dds <- DESeqDataSet(se, design = ~ cell + dex)

为监督式学习，我们将使用bioconductor包中的宫颈计数数据，MLSeq．该数据集包含人类样本的714个miRNA的表达。宫颈肿瘤样本29例，非肿瘤样本29例。出于学习的目的，我们可以将它们视为两个独立的组，并运行各种分类算法。

library(MLSeq) filepath = system.file("extdata/宫颈.txt"， package = "MLSeq")table(filepath, header = TRUE)

无监督学习

无监督学习(Unsupervised Learning)是一套统计工具，适用于这样一种情况:我们只有一组“p”个特征，通过“n”个观察结果进行测量。我们主要感兴趣的是发现关于“p”特征的有趣的东西。

无监督学习通常是探索性数据分析的一部分。这些工具帮助我们更好地了解数据集。与监督学习问题不同，在监督学习问题中，我们可以使用预测来获得对我们的学习算法的一些信心，但没有办法检查我们的模型。因此，这种学习算法被恰当地命名为“无监督”。

RLOG转换

许多用于多维数据探索性分析的常用统计方法，特别是聚类和排序方法(如主成分分析等)，对(至少近似地)同方差数据最有效;这意味着观察量的方差(在这里，基因的表达强度)不依赖于平均值。

在RNA-Seq数据中，方差随均值增长而增长。如果一个人直接对归一化读计数矩阵进行主成分分析(PCA)，结果通常只依赖于少数表达最强的基因，因为它们在样本之间显示出最大的绝对差异。

作为一种解决方案，DESeq2提供了正则化对数转换，简称rlog。有关更多信息和选项，请参阅?rlog的帮助。

函数rlog返回一个summarizeexperiment对象，其中包含在其化验槽中经过rlog转换的值:

RLD <- rlog(dds)头(化验(RLD))

## srr1039508 srr1039509 srr1039513 srr1039516 ## ensg00000000003 9.399151 9.142478 9.501695 9.320796 9.757212 ## ensg00000000005 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 9.032567 9.063925 8.981930 ## ensg00000000457 7.949897 9.113976 9.032567 9.063925 7.834273 7.916459 7.773819 ## ensg00000000457 5.849597 5.882363 5.486937 5.770334 5.940407 ## ensg00000000938 -1.638084 -1.637483 - 1.638248 -1.636072 -1.597606 ## srr1039517 srr1039520 srr1039521 ##Ensg00000000003 9.512183 9.617378 9.315309 ## ensg00000000005 0.000000 0.000000 0.000000 0.000000 ## ensg00000000419 9.108531 8.894830 9.052303 ## ensg00000000457 7.886645 7.946411 7.908338 ## ensg00000000460 5.663847 6.107733 5.907824 ## ensg00000000938 -1.639362 -1.637608 -1.637724

评估样本之间的整体相似性:哪些样本彼此相似，哪些不同?这符合实验设计的预期吗?我们用R函数距离来计算样本之间的欧氏距离。为了避免距离度量被几个高度可变的基因所主导，并从所有基因中获得大致相等的贡献，我们在rlog转换数据上使用它

<- dist(t(化验(rld))

# # # # SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039509 40.89060 # # SRR1039512 37.35231 - 50.07638 # # SRR1039513 55.74569 41.49280 43.61052 # # SRR1039516 # # SRR1039517 57.69438 47.59326 53.52310 57.10447 41.98797 53.58929 40.99513 37.06633 51.80994 34.86653 52.54968 43.21786 46.13742 - 42.10583 # # SRR1039520 # # SRR1039521 56.04254 41.46514 - 51.90045 34.82975 - 58.40428 # # SRR1039517 SRR1039520 # # SRR1039509 # # SRR1039512 # # SRR1039513 # # SRR1039516 # # SRR1039517 # # SRR1039520 57.13688## srr1039521 47.90244 44.78171

注意使用函数t来转置数据矩阵。我们需要它是因为dist计算数据行和构成列的样本之间的距离。

的热图

我们使用函数热图，在热图中可视化样品到样品的距离。2从gplot包。注意，我们更改了距离矩阵的行名，以包含治疗类型和患者编号，而不是样本ID，以便在查看热图时查看所有这些信息。

库("gplots")库("RColorBrewer")matrix(sampleaders) rownames(sampleDistMatrix) <- paste(rld$dex, rld$cell, sep="-") colors <- colorrampalette (rev(brewer.))pal(9， "Blues")))(255) hc <- hclust(sampledist)热图。2(sampleDistMatrix, Rowv=as.dendrogram(hc)， symm=TRUE, trace="none"， col=colors, margin =c(2,10)， labCol=FALSE)

主成分分析

另一种可视化样本间距离的方法是主成分分析(PCA)。在这种排序方法中，数据点(即这里的样本)被投影到2D平面上，以便它们在两个方向上展开，这解释了数据中的大部分差异。x轴是分离数据点最多的方向(或主成分)。包含在方向上的总方差的数量印在轴标签上。在这里，我们使用了DESeq2附带的plotPCA函数。intgroup指定的两个术语是用于标记样本的有趣组;它们告诉函数使用它们来选择颜色。

plotPCA(rld, intgroup = c("dex"， "cell"))

从这两种可视化图中，我们可以看到细胞之间的差异是相当大的，但并不比地塞米松治疗引起的差异更大。这说明了为什么在差异检测中使用配对设计(“配对”，因为每个dex处理的样本与来自同一细胞系的一个未处理的样本配对)将这一点考虑到是很重要的。在一开始设置数据对象时，我们已经为此设置了设计公式~ cell + dex。

MDS另一个图，非常类似于PCA图，可以使用以r为基数的多维缩放(MDS)函数来制作。当我们没有原始数据，但只有一个距离矩阵时，这很有用。这里我们有从rlog转换计数计算出的距离的MDS图:

库(ggplot2) mds <- data.frame(cmdscale(sampleDistMatrix)) mds <- cbind(mds, colData(rld)) qplot(X1,X2,color=dex,shape=cell,data=as.data.frame(mds))

suppressPackageStartupMessages({library(limma) library(DESeq2) library(气道)})plotMDS(assay(rld)， col=as.integer(dds$dex)， pch=as.integer(dds$cell))

监督式学习

在监督学习中，除了“p”特征，我们还在同样的n次观察中测量了a反应Y。然后，目标是使用X (n xp矩阵)对新的观测结果预测Y。

对于宫颈数据，我们知道前29个是非肿瘤样本，而后29个是肿瘤样本。我们将把它们分别编码为0和1。我们将随机抽取30%的数据作为测试集。剩余70%的数据将用作培训数据

set.seed(9) class = data.frame(condition = factor(rep(c(0,1)， c(29,29))))) nTest = ceiling(ncol(cervical) * 0.2) ind = sample(ncol(cervical)， nTest, FALSE)宫颈的，颈部的Train = as.matrix(matrix);Train + 1) classstr = data.frame(condition = class[-ind，])宫颈的，宫颈的Test = as.matrix(matrix);Test + 1) class = data.frame(condition = class[ind，])

MLSeq的目标是让计算对用户来说不那么复杂，并允许用户使用单一函数学习使用各种分类器的模型。

这个包的主要功能是classification，它需要DESeqDataSet实例形式的数据。DESeqDataSet是summarizeexperiment的一个子类，用于存储输入值、中间计算和差分表达式分析的结果。

因此，让我们为训练集和测试集创建DESeqDataSet对象，并在其上运行DESeq。

颈椎。trainS4 = DESeqDataSetFromMatrix(countData = cervical.)训练，colData = classtr，公式(~condition))

将计数转换为整数模式

颈椎。颈部;颈部;trainS4, fitType = "local")

##估计尺寸因素##估计分散度##基因分散估计##平均分散关系##最终分散估计##拟合模型和测试##—替换异常值并对72个基因进行重新设计##—DESeq参数'minReplicatesForReplace' = 7 ##—原始计数保存在计数(dds)中##估计分散度##拟合模型和测试

颈椎。testS4 = DESeqDataSetFromMatrix(countData =宫颈数据。test, colData = classts，公式(~condition))

将计数转换为整数模式

颈椎。testS4 = DESeq(宫颈;testS4, fitType = "local")

##估计尺寸因素##估计分散度##基因分散估计##平均分散关系##最终分散估计##拟合模型和测试##—替换异常值并对55个基因进行重新设计##—DESeq参数'minReplicatesForReplace' = 7 ##—原始计数保存在计数(dds)中##估计分散度##拟合模型和测试

使用支持向量机进行分类。

支持向量机=分类(数据=颈椎。trainS4, method = "svm"， normalize = "deseq"， deseqTransform = "vst"， cv = 5, rpt = 3, ref = "1")

##发现已经估计的分散，替换这些##基因分散估计##平均分散关系##最终分散估计##加载所需的包:kernlab ## ##附加包:'kernlab' ## ##以下对象从'package:Biostrings': ## ##类型被屏蔽

支持向量机

## ##类MLSeq的对象## ##方法:svm ## ##准确率(%):97.83 ##灵敏度(%):100 ##特异性(%):95 ## ##参考类:1

它返回一个类为' MLseq '的对象，我们观察到它成功地拟合了一个模型，准确率为97.8%。我们可以通过。来访问这个S4对象的插槽

getSlots(“MLSeq”)

## method deseqTransform normalization confusionMat ##“character”“character”“character”“confusionMatrix”##训练ref ##“train”“character”

还有，询问训练的模特。

训练(支持向量机)

##支持向量机径向基函数核## ## 46个样本## 714预测器## 2类:'1'，'0' ## ##没有预处理##重采样:交叉验证(5次，重复3次)## ##样本大小总结:37,37,36,37,37,37,37，…## ##跨调优参数的重采样结果:## ## C精度Kappa精度SD Kappa SD ## 0.25 0.5644444 0.0000000 0.01840175 0.0000000 ## 0.50 0.7955556 0.5553664 0.10836284 0.2471095 ## 1.00 0.8770370 0.7329642 0.11396855 0.2544595 ## ##调优参数“sigma”保持在0.001113943的值不变##精度用于使用最大值选择最佳模型。该模型的最终值为sigma = 0.001113943和C = 1。

我们可以使用predict来预测测试数据的类标签

pred。支持向量机= predictClassify(svm, cervical.testS4)

已经发现的估计色散，取代这些基因方面的色散估计#平均色散关系#最终色散估计

表(pred。支持向量机, relevel(cervical.testS4$condition, 2))

## ## pred。支持向量机1 0 ## 1 3 1 ## 0 0 8

其他可用的分类方法是“随机森林”，“推车”和“袋式支持向量机”。

Rf =分类(数据=颈椎。trainS4, method = "randomforest"， normalize = "deseq"， deseqTransform = "vst"， cv = 5, rpt = 3, ref = "1")

已经发现的估计色散，取代这些基因方面的色散估计#平均色散关系#最终色散估计

训练(rf)

##随机森林## ## 46个样本## 714预测器## 2类:'1'，'0' ## ##没有预处理##重采样:交叉验证(5次，重复3次)## ##样本大小总结:37,37,36,37,37,37,37，…## ##跨调优参数的重采样结果:## ## mtry精度Kappa精度SD Kappa SD ## 2 0.8044444 0.5877801 0.11061795 0.2326000 ## 37 0.8614815 0.7105296 0.11399949 0.2379795 ## 713 0.8474074 0.6871242 0.09186167 0.1882186 ## ##精度用于选择使用最大值的最优模型。模型的最终值是mtry = 37。

pred。rf = predictclassification (rf, cervical.testS4)

已经发现的估计色散，取代这些基因方面的色散估计#平均色散关系#最终色散估计

表(pred。射频,relevel(颈。testS4 $条件,2))

## ## pred。Rf 1 0 ## 1 3 1 ## 0 0 8

SessionInfo

sessionInfo ()

## R版本3.2.0 alpha (2015-03-25 r68090) ##平台:x86_64-unknown-linux-gnu(64位)##运行在Ubuntu 14.04.2 LTS下## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_US。UTF-8 LC_COLLATE= c# [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# [9] LC_ADDRESS=C LC_TELEPHONE= c# [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 parallel stats graphics grDevices utils datasets ## [8] methods base ## ##其他附加的包:# # # # [1] kernlab_0.9-20 [2] RColorBrewer_1.1-2 # # [3] gplots_2.16.0 # # [4] MLSeq_1.3.0 # # [5] edgeR_3.9.15 # # [6] randomForest_4.6-10 # # [7] limma_3.23.12 # # [8] caret_6.0-41 # # [9] ggplot2_1.0.1 # # [10] lattice_0.20-31 # # [11] DESeq2_1.7.50 # # [12] RcppArmadillo_0.4.650.1.1 # # [13] Rcpp_0.11.5 # # [14] airway_0.101.1 # # [15] genefilter_1.49.2 # # [16] TxDb.Hsapiens.UCSC.hg19.knownGene_3.1.2 # # [17] GenomicFeatures_1.19.36 # # [18] AnnotationDbi_1.29.21 # # [19] Biobase_2.27.3 # # [20]## [25] Biostrings_2.35.12 ## [26] XVector_0.7.4 ## [27] genomic icranges_1.19.52 ## [28] GenomeInfoDb_1.3.16 ## [29] IRanges_2.1.43 ## [30] S4Vectors_0.5.22 ## [31] BiocGenerics_0.13.11 ## [32] BiocStyle_1.5.3 ## ##通过命名空间加载(并没有附加):# [7] colorspace_1.2-6 nnet_7.3-9 bradleyterry2_1 .8-6 ## [13] compiler_3.2.0 quantreg_5.11 formatR_1.1 ## [16] SparseM_1.6 labeling_0.3 caTools_1.17.1 ## [19] scales_0.2.4 brglm_0.5-9 stringgr_0.6.2 ## [22] digest_0.6.8 foreign_0.8-63 minqa_1.2.4 ## [25] rmarkdown_0.5.1 htmltools_0.2.6 lme4_1.1-7 ## [28] RSQLite_1.0.0 BiocParallel_1.1.21 gtools_3.4.1 ## [7] Hmisc_3.15-0 DBI_0.3.1 kkernsmooth_2 - 23-14 ## [7] Hmisc_3.15-0 DBI_0.3.1 mgcv_1.8-6 ## [10] colorspace_1.2-6 ## [10] colorspace_1.2-6 ## [13] compiler_3.2.0 quantreg_5.11 formatR_1.1 ## [13]acepack_1.3-3.3 car_2.0-25 RCurl_1.95-4.5 ## [34] Formula_1.2-0 futile.logger_1.4 Matrix_1.2-0 ## [37] munsell_0.4.2 proto_0.3-10 yaml_2.1.13 ## [40] MASS_7.3-40 zlibbioc_1.13.3 plyr_1.8.1 ## [43] grid_3.2.0 gdata_2.13.3 splines_3.2.0 ## [46] annotate_1.45.4 locfit_1.5-9.1 knitr_1.9 ## [49] geneplotter_1.45.0 reshape2_1.4.1 codetools_0.2-11 ## [52] biomaRt_2.23.5 futile.options_1.0.0 XML_3.98-1.1 ## [55] evaluate_0.5.5 latticeExtra_0.6-26 lambda.r_1.1.7 ## [58] nloptr_1.0.4 foreach_1.4.2 gtable_0.1.2 ## [61] xtable_1.7-4 e1071_1.6-4 class_7.3-12 ## [64] survival_2.38-1 snow_0.3-13 iterators_1.0.7 ## [67] cluster_2.0.1

参考文献

1. Zararsiz G, Goksuluk D, Korkmaz S, Eldem V, Duru IP, Unver T和Ozturk A(2014)。MLSeq:用于RNA-Seq数据的机器学习接口。R包版本1.3.0。
2. 希姆斯，E. B，江，X.，瓦格纳，P.，胡，R.，王，Q.， Klanderman, B.，惠特克，M. R.，段，Q.， Lasky-Su, J.，尼古拉斯，C.，杰斯特，W.，约翰逊，M.， Panettieri, A. R.，坦蒂希拉，G. K，韦斯，T. S，卢和Q.(2014)。“RNA-Seq转录组分析鉴定CRISPLD2为糖皮质激素响应基因，调节气道平滑肌细胞中的细胞因子功能。”科学通报，29 (6)，pp. e99625。http://www.ncbi.nlm.nih.gov/pubmed/24926665．
3. 《统计学习与R语言应用导论》，加雷思·詹姆斯，丹妮拉·威滕，特雷弗·哈斯蒂和罗伯特·蒂布希拉尼
4. 统计学习的要素:数据挖掘，推断和预测。特雷弗·哈斯蒂，罗伯特·蒂布希拉尼，杰罗姆·弗里德曼