edgeR-DESeq2分析RNA-seq差异表达

edgeR 包的安装

edgeR 包是基于 Bioconductor 平台发布的，所以安装不能直接用 install.packages() 命令从 CRAN 上来下载 ? 安装:

?

# try http:// if https:// URLs are not supported >source(\) >biocLite(\)

数据导入

由于 edgeR 对测序结果的下游分析是依赖 count 计数来进行基因差异表达分析的，在这里使用的是featureCounts 来进行统计 ｀.bam｀ 文件中 Map 的结果 ? count 结果如下:

?

>library(edgeR)

>mydata<-read.table(\,header=TRUE,quote='\\t',skip=1) >sampleNames<-c(\,\,\,\,\,\) >names(mydata)[7:12]<-sampleNames >head(mydata)

GeneidChrStartEndStrandLengthCA_1CA_2CA_3CC_1CC_2CC_3 1gene1314NW_139421.112571745+489000000 2gene1315NW_139421.121153452+1338000000 3gene1316NW_139421.138564680+825000000 4gene1317NW_139421.148665435-570000000 5gene1318NW_139421.160666836-771000000 6gene1319NW_139421.172949483+2190000000

?

在这里我们只是需要 Geneid 和后 6 列的样本的 count 信息来组成矩阵，所以要处理下

>countMatrix<-as.matrix(mydata[7:12]) >rownames(countMatrix)<-mydata$Geneid >head(countMatrix)

CA_1CA_2CA_3CC_1CC_2CC_3 gene1314000000 gene1315000000 gene1316000000 gene1317000000 gene1318000000 gene1319000000

*要导入的矩阵由3v3样本组成(三组生物学重复)

创建 DEGlist

>group<-factor(c(\,\,\,\,\,\)) >y<-DGEList(counts=countMatrix,group=group) >y

Anobjectofclass\ $counts

CA_1CA_2CA_3CC_1CC_2CC_3 gene1314000000 gene1315000000 gene1316000000 gene1317000000 gene1318000000 14212morerows...

$samples

grouplib.sizenorm.factors CA_1CA_117885371 CA_2CA_218255461 CA_3CA_319030171 CC_1CC_118260421 CC_2CC_221244681 CC_3CC_320250631

过滤

?

过滤掉那些 count 结果都为0的数据，这些没有表达的基因对结果的分析没有用，过滤又两点好处:

1 可以减少内存的压力 2 可以减少计算的压力

>keep<-rowSums(cpm(y)>1)>=2 >y<-y[keep,,keep.lib.sizes=FALSE] >y

Anobjectofclass\ $counts

CA_1CA_2CA_3CC_1CC_2CC_3 gene1321161138129218194220 gene1322231133 gene1323202733475146 gene132460877986100132 gene1325322921587556 3877morerows...

$samples

grouplib.sizenorm.factors CA_1CA_117883621 CA_2CA_218253081 CA_3CA_319027961 CC_1CC_118258891 CC_2CC_221241551 CC_3CC_320247861

标准化处理

?

edgeR采用的是 TMM 方法进行标准化处理,只有标准化处理后的数据才又可比性

>y<-calcNormFactors(y) >y

Anobjectofclass\ $counts

CA_1CA_2CA_3CC_1CC_2CC_3

gene1321161138129218194220 gene1322231133 gene1323202733475146 gene132460877986100132 gene1325322921587556 3877morerows...

$samples

grouplib.sizenorm.factors CA_1CA_117883620.9553769 CA_2CA_218253080.9052539 CA_3CA_319027960.9686232 CC_1CC_118258890.9923455 CC_2CC_221241551.1275178 CC_3CC_320247861.0668754

设计矩阵

?

为什么要一个设计矩阵呢，道理很简单，有了一个设计矩阵才能够更好的分组分析

>subGroup<-factor(substring(colnames(countMatrix),4,4)) >design<-model.matrix(~subGroup+group) >rownames(design)<-colnames(y) >design

(Intercept)subGroup2subGroup3groupCC CA_11000 CA_21100 CA_31010 CC_11001 CC_21101 CC_31011 attr(,\) [1]0112

attr(,\)

attr(,\)$subGroup [1]\

attr(,\)$group

[1]\

评估离散度

>y<-estimateDisp(y,design,robust=TRUE) >y$common.dispersion

[1]0.02683622 #plot >plotBCV(y)

差异表达基因

>fit<-glmQLFit(y,design,robust=TRUE) >qlf<-glmQLFTest(fit) >topTags(qlf) Coefficient:groupCC logFClogCPMFPValueFDR

gene7024-5.5156489.612809594.92326.431484e-442.496702e-40 gene66125.1302828.451143468.20601.557517e-393.023140e-36

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库edgeR-DESeq2分析RNA-seq差异表达在线全文阅读。