基本的问题

一些活动与RNA-seq数据的聚类分析和分类有关

你将组成5个4-5人的小组。你将花5-7分钟来概述以下问题的书面攻击计划，然后提交你的计划以供讨论和演示。

给定来自129个个体的12.9K基因水平读数，代表一个不确定数量的群体

我们应该修改数据吗?
定义一种计算总体数量和分配样本的方法。如何比较各种方法?
事实:从每个总体中大约有相同数量的样本。

假设您得到了人口标签。阐述提取区分人群的表达谱的原则。建议如何发展生物解释的剖面。
定义一个分类程序，可以推断人口来源的基础上一个或多个转录谱。这个算法有多好?

解决的尝试

的基本数据。推导了以后。

负载(“montpick_eset.RData”)

montpick.eset

##错误:没有找到montpick.eset对象

删除纯0。

pz = apply(exprs(montpick.eset)，1, function(x) all(x==0)) mnz = montpick.eset[-which(pz)，]

还需要一些额外的擦洗…

暗(mnz)

##功能样本## 12160 129

简单的聚类方法。

hc = hclust (dist (t (exprs (mnz))))

表(cutree (hc, k = 2))

## ## 12 ## 128 1

表(cutree (hc, k = 3))

## ## 1 2 3 ## 84 1

表(cutree (hc, k = 4))

## ## 1 2 3 4 ## 84 33 11 1

在移位一之后尝试简单的对数。

hc。l = hclust (dist (t(日志(exprs (mnz) + 1))))

表(cutree (hc。l, k = 2))

## ## 12 ## 122 7

表(cutree (hc。l, k = 3)

## ## 1 2 3 ## 30 92 7

表(cutree (hc。l, k = 4))

## ## 1 2 3 4 ## 30 92 6 1

试试PAM:

pam库(集群)。2 = pam(dist(t(log(exprs(mnz)+1))))， k=2)3 = pam(dist(t(log(exprs(mnz)+1))))， k=3)4 = pam (dist (t(日志(exprs (mnz) + 1))), k = 4)表(pam.4 clust美元)

## ## 1 2 3 4 ## 45 20 49 15

试试k - means:

公里。2 = kmeans(t(log(exprs(mnz)+1))， 2) km。3 = kmeans(t(log(exprs(mnz)+1))， 3) table(km.2$clust)

## ## 1 2 ## 77 52

表(km.3 clust美元)

## ## 1 2 3 ## 33 60 36

试试kmeans的音量转换计数。

library(limma) vmnz = voom(exprs(mnz)) txd = vmnz$E lkm2v = kmeans(t(txd)，2) lkm4v = kmeans(t(txd)，4) table(lkm2v$clust)

## ## 1 2 ## 63 66

表(lkm4v clust美元)

## ## 1 2 3 4 ## 15 43 60 11

考虑“clusterStability”措施。

尝试DESeq分析。

se = summarizeexperimental (assays=list(counts=exprs(mnz))) colData(se) = DataFrame(pData(mnz)) de = DESeqDataSet(se， ~population) derun = DESeq(de)

# # # # # #估计大小因素估计分散gene-wise分散估计# # # #最后分散估计平均分散关系# # # #拟合模型和测试——271个基因替换异常值和改装# #——DESeq论证“minReplicatesForReplace”= 7 # #,保留原始计数(dds)) # #估计拟合模型和测试

res =结果(derun)

summarizeResults (res)

# # # #与非零读总数12160 # # # #在调整p值< 0.1:# # # #利物浦> 0(上):3143年,26% # #利物浦< 0(下降):2975年,24% # # # #过滤*:1216年,10% # #标记* *:0,0 % ## ## * 低意味着计数,看到independentFiltering争论的结果()# # * *高库克的距离,看到cooksCutoff争论的结果()

res ((res padj < .005),美元)

## log2 fold change (MAP): population YRI vs CEU ## Wald test p-value:与3804年人口套vs CEU # # DataFrame行6列# # baseMean log2FoldChange lfcSE stat pvalue # # <数字> <数字> <数字> <数字> <数字> # # ENSG00000000419 e-08 # # ENSG00000000457 41.785 0.5559 0.09824 5.658 1.528 48.692 0.4202 0.09047 4.645 - 3.400 e-06 # # ENSG00000002330 65.562 -0.6426 0.08949 -7.180 6.966 e-13 # #ENSG00000002834 155.210 -0.6368 0.10974 -5.803 6.516e-09 ## ... ... ... ... ... ...# # ENSG00000253873 e-05 # # ENSG00000253954 8.008 -0.6856 0.1749 -3.920 8.849 9.913 -0.6413 0.1505 -4.262 2.023 e-05 # # ENSG00000254004 e-09 # # ENSG00000254128 4.520 1.0417 0.1723 6.045 1.492 9.615 1.0443 0.1772 5.894 3.774 e-09 # # ENSG00000254206 1.843 2.2559 0.3519 6.410 1.452平台以及# # padj # # <数字> 1.151 e-07 # # # # ENSG00000000419 ENSG000000004571.732e-05 ## ENSG00000002330 1.015e-11 ## ENSG00000002745 1.007e-04 ## ENSG00000002834 5.190e-08 ## ... ...## ENSG00000253873 3.470e-04 ## ENSG00000253954 8.945e-05 ## ENSG00000254004 1.331e-08 ## ENSG00000254206 3.120e-08 ## ENSG00000254206 1.511e-09

这对技术工件敏感吗?对日志转换后的数据使用SVA并进行检查。

(1) . sv1 = sva(log(exprs(mnz)+1)， mod=model.matrix(log(exprs(mnz))+1)， n.sv=num)， n.sv=num

##迭代(5):1 2 3 4 5

打印(num)

# # 1 [1]

总结(sv1 sv)美元

##最小1区。中位数3区。最大值。## -0.1600 -0.0633 -0.0220 0.0000 0.0503 0.3120

secorr = se secorr$sv = sv1$sv

derun2 = DESeq(decorr， ~sv+population)

##估计尺寸因子##估计离散度##基因离散度估计##均值-离散关系##最终离散度估计##拟合模型和检验

相关系数=结果(derun2)

summarizeResults (corr)

# # # #与非零读总数12160 # # # #在调整p值< 0.1:# # # #利物浦> 0(上):3161年,26% # #利物浦< 0(下降):3006年,25% # # # #过滤*:1216年,10% # #标记* *:0,0 % ## ## * 低意味着计数,看到independentFiltering争论的结果()# # * *高库克的距离,看到cooksCutoff争论的结果()

收款银行((corr padj < .005),美元)

## log2 fold change (MAP): population YRI vs CEU ## Wald test p-value:与3818年人口套vs CEU # # DataFrame行6列# # baseMean log2FoldChange lfcSE stat pvalue # # <数字> <数字> <数字> <数字> <数字> # # ENSG00000000419 e-09 # # ENSG00000000457 41.785 0.5780 0.09948 5.810 6.247 48.692 0.3842 0.09019 4.260 - 2.040 e-05 # # ENSG00000002330 65.562 -0.6635 0.09199 -7.213 5.487 e-13 # #ENSG00000002834 155.210 -0.7126 0.11082 -6.430 1.277e-10 ## ... ... ... ... ... ...# # ENSG00000253873 e-04 # # ENSG00000253954 8.008 -0.6432 0.1780 -3.612 3.034 9.913 -0.6010 0.1526 -3.938 8.220 e-05 # # ENSG00000254004 e-09 # # ENSG00000254128 4.520 1.0304 0.1752 5.880 4.095 9.615 1.0214 0.1800 5.673 1.400 e-08 # # ENSG00000254206 1.843 2.3030 0.3464 6.647 2.983 e-11 # # padj # # <数字> 4.947 e-08 # # # # ENSG00000000419 ENSG000000004578.866e-05 ## ENSG00000002330 7.980e-12 ## ENSG00000002745 8.639e-05 ## ENSG00000002834 1.302e-09 ## ... ...## ENSG00000253873 1.045e-03 ## ENSG00000253954 3.181e-04 ## ENSG00000254004 3.337e-08 ## ENSG00000254128 1.048e-07 ## ENSG00000254206 3.281e-10

方法分类

库(randomGLM)

##加载required package: MASS ## ##加载package: 'MASS' ## ##以下对象是从'package:genefilter': ## ## area ## ##以下对象是从'package:AnnotationDbi': ## ## select ## ##加载required package: foreach ## foreach:使用Revolution R实现可伸缩性、容错等功能。## http://www.revolutionanalytics.com ## # Loading required package: doparlesle# #加载required package: iterators

test = mnz[，-sort(mysamp)] y = 1*(train$pop == "YRI") *(train$pop == "YRI")

##警告:名称在数据帧中部分匹配

rglm1 = randomGLM(t(exprs(train))， y, xtest=t(exprs(test))))

测试表(rglm1 predictedTest美元,美元流行)

##警告:名称在数据帧中部分匹配

## ## CEU YRI ## 0 22 3 ## 1 0 29

集成学习的另一种方法

库(bigrf)

##加载所需的包:bigmemory。##加载所需的包:BH ## ## bigmemory >= 4.0是自3.1.2以来的主要修订;更多信息请参见软件包## biganalytics和bigtabulate以及http://www.bigmemory.org。

isyri = as.integer(1 *(火车流行美元= =“套”)+ 1)

##警告:名称在数据帧中部分匹配

bb = bigrfc(data.frame(t(exprs(train)))， isyri)

## OOB错误:## Tree Overall error error by class ## 12 ## 10 12.00 5.26 18.92 ## 20 10.67 7.89 13.51 ## 30 6.67 5.26 8.11 ## 40 6.67 2.63 10.81 ## 50 5.33 2.63 8.11

bigpr = predict(bb, x=data.frame(t(exprs(test))))))

##处理树编号:## 10 ## 20 ## 30 ## 40 ## 50

表(bigpr、测试流行美元)

##警告:名称在数据帧中部分匹配

## ## bigpr CEU YRI ## 1 22 3 ## 2 0 29