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

1、edgeR 包的安装edgeR 包是基于 Bioconductor 平台发布的，所以安装不能直接用install.packages()命令从 CRAN 上来下载安装：#tryhttp:/ifhttps:/URLsarenotsupportedsource(/biocLite.R)biocLite(edgeR)数据导入由于 edgeR 对测序结果的下游分析是依赖 count 计数来进行基因差异表达分析的，在这里使用的是 featureCounts 来进行统计.bam文件中 Map 的结果count 结果如下：library(edgeR)mydatas

2、ampleNamesnames(mydata)7:12head(mydata)GeneidChrStartEndStrandLengthCA1CA2CA3CC1CC2CC3在这里我们只是需要 Geneid 和后 6 列的样本的 count 信息来组成矩阵，所以要处理下countMatrixrownames(countMatrix)head(countMatrix)1gene1314NW_139421.12gene1315NW_139421.13gene1316NW_139421.14gene1317NW_139421.15gene1318NW_139421.16gene1319NW_13942

3、1.112571745+48900000021153452+133800000038564680+82500000048665435-57000000060666836-77100000072949483+2190000000CA_1CA_2CA_3CC_1CC_2CC_3gene1314000000gene1315000000gene1316000000gene1317000000gene1318000000gene1319000000*要导入的矩阵由 3V3 样本组成(三组生物学重复)创建 DEGlistgroupyyAnobjectofclassDGEList$countsCA_1CA_

4、2CA_3CC_1CC_2CC_3gene1314000000gene1315000000gene1316000000gene1317000000gene131800000014212morerows$samplesgrouplib.sizenorm.factorsCA_1CA_117885371CA_2CA_218255461CA_3CA_319030171CC_1CC_118260421CC_2CC_221244681CC_3CC_320250631过滤过滤掉那些 count 结果都为 0 的数据，这些没有表达的基因对结果的分析没有用，过滤又两点好处：1 可以减少内存的压力 2 可以减少计

5、算的压力keep1)=2yyAnobjectofclassDGEList$countsCA1CA2CA3CC1CC2CC3gene1321161138129218194220gene1322231133gene1323202733475146gene132460877986100132gene13253229215875563877morerows.$samplesgrouplib.sizenorm.factorsCA_1CA_117883621CA_2CA_218253081CA_3CA_319027961CC_1CC_118258891CC_2CC_221241551CC_3CC_3202

6、47861标准化处理edgeR 采用的是 TMM 方法进行标准化处理，只有标准化处理后的数据才又可比性yyAnobjectofclassDGEList$countsCA1CA2CA3CC1CC2CC3gene1321161138129218194220gene1322231133gene1323202733475146gene132460877986100132gene13253229215875563877morerows.$samplesgrouplib.sizenorm.factorsCA_1CA_117883620.9553769CA_2CA_218253080.9052539CA_3

7、CA_319027960.9686232CC_1CC_118258890.9923455CC_2CC_221241551.1275178CC_3CC_320247861.0668754设计矩阵为什么要一个设计矩阵呢，道理很简单，有了一个设计矩阵才能够更好的分组分析subGroupdesignrownames(design)design(Intercept)subGroup2subGroup3groupCCCA_11000CA_21100CA_31010CC_11001CC_21101CC_31011attr(,assign)10112attr(,contrasts)attr(,contrast

8、s)$subGroup1contr.treatmentattr(,contrasts)$group1contr.treatment评估离散度yy$common.dispersion10.02683622#plotplotBCV(y)QF055E2eCMO0S1015AveragelogCpM差异表达基因fitqlftopTags(qlf)Coefficient:groupCClogFClogCPMFPValueFDRgene7024-5.5156489.612809594.92326.431484e-442.496702e-40gene66125.1302828.451143468.20601

9、.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.0966

10、19e-21gene62422.5291259.897771169.26583.022914e-231.466869e-20gene72523.7323156.137670188.20943.890569e-231.678132e-20gene61252.8754236.569935160.31891.656083e-226.428914e-20查看差异表达基因原始的 CMPtopcpm(y)top,CA_1CA_2CA_3CC_1CC_2CC_3gene70241711.3830021405.8618991480.12111533.1141837.1604029.62696gene66121

11、7.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.90560gen

12、e2611245.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(dtisDEDEnameshead(DEnames)1gene1

13、325gene1326gene1327gene1331gene1340gene1343差异表达基因画图plotSmear(qlf,de.tags=DEnames)abline(h=c(-1,1),col=blue)DESeq2 包的安装安装：#tryhttp:/ifhttps:/URLsarenotsupportedsource(/biocLite.RbiocLite(DESeq2)数据导入导入 count 矩阵，导入数据的方式很多这里直接导入 count 矩阵count 结果如下：library(DESeq2)sampleNames-c(CA_1

14、,CA_2,CA_3,CC_1,CC_2,CC_3)mydata-read.table(counts.txt,header=TRUEquote=t,skip=1)names(mydata)7:12-sampleNamescountMatrix-as.matrix(mydata7: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)ddsddsclass:DESeqData

15、Setdim:142176metadata(0):assays(1):countsrownames(14217):gene1314gene1315gene6710gene6709rowRangesmetadatacolumnnames(0):colnames(6):CA_1CA_2CC_2CC_3colDatanames(2):namecondition过滤过滤掉那些 count 结果都为 0 的数据，这些没有表达的基因对结果的分析没有用dds1,ddsclass:DESeqDataSetdim:41906metadata(0):assays(1):countsrownames(4190):g

16、ene1321gene1322gene6712gene6710rowRangesmetadatacolumnnames(0):colnames(6):CA_1CA_2CC_2CC_3colDatanames(2):nameconditionPC 矽析rld-rlog(dds)plotPCA(rld,intgroup=c(name,condition)PCl:的%Marine当然也可以使用 ggplot2 来画 PCA 图library(ggplot2)rld-rlog(dds)data-plotPCA(rld,intgroup=c(condition,name),returnData=TRUE

17、)percentVar-round(100*attr(data,percentVar)p-ggplot(data,aes(PC1,PC2,color=condition,shape=name)+geom_point(size=3)+xlab(paste0(PC1:,percentVar1,%variance)+ylab(paste0(PC2:,percentVar2,%variance)p注意在进行 PCA 分析前不要library(DESeq)5-groupCA_1:CACA_2：CACA3:GAGG_1:CCGC_2:0CCC_3lOC否则无法进行 PCA 分析PC1:variance差异

18、表达基因分析分析结果输出library(DESeq)dds-DESeq(dds)res-results(dds)write.table(res,result.csv,sep=,s=TRUE)head(res)10g2foldchange(MAP):conditionCCvsCAWaldtestp-value:conditionCCvsCADataFramewith6rowsand6columnsbaseMeanlog2FoldChangelfcSEstatpvaluegene1321173.2886810.262679590.20499831.28137422.000623e-

19、01condition*CA*ecCA_I*CA_2CA_3+GG.1闰QC_2*CCJJgene13222.118367-0.052379520.4989589-0.10497769.163936e-01gene132335.9737010.500545800.30380961.64756419.944215e-02gene132488.4216610.176776050.24027270.73573094.618945e-01gene132543.0018280.811431040.29193962.77944865.445127e-03gene1326662.136259-1.05356

20、1050.1752230-6.01268801.824720e-09padjgene13213.790396e-01gene13229.559679e-01gene13232.337858e-01gene13246.565731e-01gene13252.447141e-02gene13264.520861e-08注：（1）rownames:基因 ID（2）baseMean:所有样本矫正后的平均 reads 数（3）log2FoldChange:取 log2 后的表达量差异（4）pvalue:统计学差异显著性检验指标（5）padj:校正后的 pvalue,padj 越小，表示基因表达差异越显著summary查看整体分析结果summary(res)outof4190withnonzerot