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

edgeR包的安装edgeR包是基于Bioconductor平台发布的，所以安装不能直接用install.packages()命令从CRAN上来下载安装:#try()biocLite(edgeR)数据导入由于edgeR对测序结果的下游分析是依赖count计数来进行基因差异表达分析的，在这里使用的是featureCounts来进行统计｀.bam｀文件中Map的结果count结果如下:library(edgeR)mydata-read.table(counts.txt,header=TRUE,quote='\t',skip=1)sampleNames-c(CA_1,CA_2,CA_3,CC_1,CC_2,CC_3)names(mydata)[7:12]-sampleNameshead(mydata)GeneidChrStartEndStrandLengthCA_1CA_2CA_3CC_1CC_2CC_31gene1314NW_139421.112571745+4890000002gene1315NW_139421.121153452+13380000003gene1316NW_139421.138564680+8250000004gene1317NW_139421.148665435-5700000005gene1318NW_139421.160666836-7710000006gene1319NW_139421.172949483+2190000000在这里我们只是需要Geneid和后6列的样本的count信息来组成矩阵，所以要处理下countMatrix-as.matrix(mydata[7:12])rownames(countMatrix)-mydata$Geneidhead(countMatrix)CA_1CA_2CA_3CC_1CC_2CC_3gene1314000000gene1315000000gene1316000000gene1317000000gene1318000000gene1319000000*要导入的矩阵由3v3样本组成(三组生物学重复)创建DEGlistgroup-factor(c(CA,CA,CA,CC,CC,CC))y-DGEList(counts=countMatrix,group=group)yAnobjectofclassDGEList$countsCA_1CA_2CA_3CC_1CC_2CC_3gene1314000000gene1315000000gene1316000000gene1317000000gene131800000014212morerows...$samplesgrouplib.sizenorm.factorsCA_1CA_117885371CA_2CA_218255461CA_3CA_319030171CC_1CC_118260421CC_2CC_221244681CC_3CC_320250631过滤过滤掉那些count结果都为0的数据，这些没有表达的基因对结果的分析没有用，过滤又两点好处:1可以减少内存的压力2可以减少计算的压力keep-rowSums(cpm(y)1)=2y-y[keep,,keep.lib.sizes=FALSE]yAnobjectofclassDGEList$countsCA_1CA_2CA_3CC_1CC_2CC_3gene1321161138129218194220gene1322231133gene1323202733475146gene132460877986100132gene13253229215875563877morerows...$samplesgrouplib.sizenorm.factorsCA_1CA_117883621CA_2CA_218253081CA_3CA_319027961CC_1CC_118258891CC_2CC_221241551CC_3CC_320247861标准化处理edgeR采用的是TMM方法进行标准化处理,只有标准化处理后的数据才又可比性y-calcNormFactors(y)yAnobjectofclassDGEList$countsCA_1CA_2CA_3CC_1CC_2CC_3gene1321161138129218194220gene1322231133gene1323202733475146gene132460877986100132gene13253229215875563877morerows...$samplesgrouplib.sizenorm.factorsCA_1CA_117883620.9553769CA_2CA_218253080.9052539CA_3CA_319027960.9686232CC_1CC_118258890.9923455CC_2CC_221241551.1275178CC_3CC_320247861.0668754设计矩阵为什么要一个设计矩阵呢，道理很简单，有了一个设计矩阵才能够更好的分组分析subGroup-factor(substring(colnames(countMatrix),4,4))design-model.matrix(~subGroup+group)rownames(design)-colnames(y)design(Intercept)subGroup2subGroup3groupCCCA_11000CA_21100CA_31010CC_11001CC_21101CC_31011attr(,assign)[1]0112attr(,contrasts)attr(,contrasts)$subGroup[1]contr.treatmentattr(,contrasts)$group[1]contr.treatment评估离散度y-estimateDisp(y,design,robust=TRUE)y$common.dispersion[1]0.02683622#plotplotBCV(y)差异表达基因fit-glmQLFit(y,design,robust=TRUE)qlf-glmQLFTest(fit)topTags(qlf)Coefficient:groupCClogFClogCPMFPValueFDRgene7024-5.5156489.612809594.92326.431484e-442.496702e-40gene66125.1302828.451143468.20601.557517e-393.023140e-36gene27434.3774925.586773208.02683.488383e-264.513967e-23gene120324.7343835.098148192.93784.359649e-254.231040e-22gene491-2.73391010.412673190.98396.104188e-254.739291e-22gene89412.9971856.839106177.76146.332836e-244.097345e-21gene2611-2.8469247.216173174.73321.099339e-236.096619e-21gene62422.5291259.897771169.26583.022914e-231.466869e-20gene72523.7323156.137670188.20943.890569e-231.678132e-20gene61252.8754236.569935160.31891.656083e-226.428914e-20查看差异表达基因原始的CMPtop-rownames(topTags(qlf))cpm(y)[top,]CA_1CA_2CA_3CC_1CC_2CC_3gene70241711.3830021405.8618991480.12111533.1141837.1604029.62696gene661217.55864912.10384826.585753403.99298582.457961044.35046gene27434.6823061.8155775.96823062.9169487.26431114.34156gene120321.7558652.4207702.71283265.6764647.5987275.45617gene4912811.1397272059.4696692222.351938444.83381385.38258253.68087gene894123.99682024.81288824.415488131.35291244.67410225.90560gene2611245.821088310.463691225.16505243.0484326.3045539.81123gene6242231.188880299.570228298.4115151348.298991343.619882191.93237gene72529.36461313.3142325.42566492.71970108.55847181.92807gene612523.41153214.52461729.841152145.70239160.75005185.16852查看上调和下调基因的数目summary(dt-decideTestsDGE(qlf))[,1]-1536027931553挑选出差异表达基因的名字isDE-as.logical(dt)DEnames-rownames(y)[isDE]head(DEnames)[1]gene1325gene1326gene1327gene1331gene1340gene1343差异表达基因画图plotSmear(qlf,de.tags=DEnames)abline(h=c(-1,1),col=blue)DESeq2包的安装安装:##try()biocLite(DESeq2)数据导入导入count矩阵，导入数据的方式很多这里直接导入count矩阵count结果如下:library(DESeq2)sampleNames-c(CA_1,CA_2,CA_3,CC_1,CC_2,CC_3)mydata-read.table(counts.txt,header=TRUE,quote='\t',skip=1)names(mydata)[7:12]-sampleNamescountMatrix-as.matrix(mydata[7:12])rownames(countMatrix)-mydata$Geneidtable2-data.frame(name=c(CA_1,CA_2,CA_3,CC_1,CC_2,CC_3),condition=(CA,CA,CA,CC,CC,CC))rownames(table2)-sampleNameshead(countMatrix)CA_1CA_2CA_3CC_1CC_2CC_3gene1314000000gene1315000000gene1316000000gene1317000000gene1318000000gene1319000000把count矩阵转化为DESeq2的数据格式dds-DESeqDataSetFromMatrix(countMatrix,colData=ta