edgeR_差异表达分析

差异表达分析

用到的R包是edgeR

在进行差异表达分析之前，会先查看一下生物学重复的数据可用性，主要会看一下MDS样本无监督聚类图和dispression离散分析图

MDS样本无监督聚类

MDS展示的是样本的无监督聚类，可以展示出样品间的相似性和不相似性，可以对于能检测到多少的差异表达基因有个大致的概念，可以查看样本的分组情况

stage <- factor(MD.targets$stage)

View(stage)

colnames(MD$counts)

colnames(MD$counts) <- c("MD1_1","MD1_2","MD1_3","MD1_4","MD2_1","MD2_2","MD2_3","MD2_4","MD3_1","MD3_2","MD3_3","MD3_4")

colnames(MD$counts)

MD.design <- model.matrix(~stage)

MD.design

## MDS画图

> plotMDS(MD)
> png("./output/mdsPlot.maleflowerRNAseq.png")
> plotMDS(MD, main="MDS plot of Male flower RNA-Seq data", labels=colnames(MD$counts))
> dev.off()
RStudioGD 
        2 

dispersion 离散分析图

MD <- estimateGLMCommonDisp(MD, MD.design, verbose=T)       

#为所有的基因都计算同样的dispersion
#Disp = 0.02097 , BCV = 0.1448 

MD <- estimateGLMTrendedDisp(MD, MD.design)     #根据不同基因的均值-方差之间的关系来计算dispersion，并拟合一个trended model

MD <- estimateGLMTagwiseDisp(MD, MD.design)     #计算每个基因的dispersion，并利用经验性贝叶斯方法压缩至Trend dispersion

summary(MD$tagwise.dispersion)
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.004798 0.008932 0.014360 0.026628 0.030889 0.910934

summary(MD$common.dispersion)
 #  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02097 0.02097 0.02097 0.02097 0.02097 0.02097 

plotBCV(MD, main="BCV plot of Tag wise dispersion estimation (prior.df=20)\nMaleflower RNA-Seq data", xlab="logCPM", ylab="Biological coefficient of variation", pch=21, cex=0.5)
#默认df=20，使用默认值会好一点，当然其他的值也可以用来尝试

plotBCV 反映不同表达量的基因与模型拟合的情况，如果模型拟合得好则tagwise点的分布会拟合到trend这条曲线上，低表达量的图就会出现离散的情况

BCV（Biological coefficient of variance ）生物学差异，是样本之间真实存在的差异，之后通过差异表达分析得到的数据其实就是在这一步中不在离散曲线上面的点

edge R是通过估计dispersion来估计BCV，进而拟合出线性回归模型的参数

差异表达分析

先对MD2与MD1进行差异分析

MD.fit <- glmQLFit(MD, MD.design)
 
colnames(MD.fit)        #查看列名，然后选择我们需要做比较的列，因为默认时拿MD1去进行比较的

MD.lrt_MD2vsMD1 <- glmQLFTest(MD.fit, coef=2)

topTags(MD.lrt_MD2vsMD1)

查看上调和下调基因的数量

>summary(decideTestsDGE(MD.lrt_MD2vsMD1, adjust.method="BH", p.value=0.05))
       stageMD2
Down       2305
NotSig    14409
Up         2364

#疑惑会不会差异表达基因有点少

>MD.lrt_MD2vsMD1cpm <- cpm(MD)[rownames(topTags(MD.lrt_MD2vsMD1, n=4669)$table)[1:4669],]

>head(MD.lrt_MD2vsMD1cpm)

#分别导出所有的基因

> write.table(topTags(MD.lrt_MD2vsMD1, n=4669), file="./output/MD.lrt_MD2vsMD1.DElist.txt", sep="\t", quote=F, row.names=T)

> write.table(topTags(MD.lrt_MD2vsMD1, n=NULL), file="./output/MD.lrt_MD2vsMD1.noDE_and_DElist.txt", sep="\t", quote=F, row.names=T)

> MD.lrt_MD2vsMD1.topTags <- topTags(MD.lrt_MD2vsMD1, n=4669)$table

> nrow(subset(MD.lrt_MD2vsMD1.topTags, logFC > 0))
[1] 2364

> MD.lrt_MD2vsMD1.UP <- subset(MD.lrt_MD2vsMD1.topTags, logFC > 0)

> write.table(MD.lrt_MD2vsMD1.UP, file="./output/MD.lrt_MD2vsMD1_UP.DElist.txt", sep="\t", quote=F, row.names=T )

> nrow(subset(MD.lrt_MD2vsMD1.topTags, logFC < 0))
[1] 2305

> MD.lrt_MD2vsMD1.DW <- subset(MD.lrt_MD2vsMD1.topTags, logFC < 0)

> write.table(MD.lrt_MD2vsMD1.DW, file="./output/MD.lrt_MD2vsMD1_DW.DElist.txt", sep="\t", quote=F, row.names=T )

根据得到差异表达基因的数量指定颜色然后进行可视化

> MD.DEtags.lrt_MD2vsMD1 <- rownames(MD)[as.logical(decideTestsDGE(MD.lrt_MD2vsMD1, adjust.method="BH", p.value=0.05))]

> MD.col_MD2vsMD1 <- rep("grey55", nrow(MD))

> MD.col_MD2vsMD1[MD.lrt_MD2vsMD1$table$PValue < 0.05 & MD.lrt_MD2vsMD1$table$logFC >0 ] <- "red"

> MD.col_MD2vsMD1[MD.lrt_MD2vsMD1$table$PValue < 0.05 & MD.lrt_MD2vsMD1$table$logFC <0 ] <- "blue"

> plotSmear(MD.lrt_MD2vsMD1, de.tags=MD.DEtags.lrt_MD2vsMD1, xlab="log-counts per million (logCPM)", ylab="log-fold change (logFC)", main="Stage 2 compare with Stage 1", pch=19, cex=0.4, smearWidth=0.5, panel.first=grid(), smooth.scatter=FALSE, ylim=c(-7,7), yaxs="i")

> abline(h=c(-1,1),col="dodgerblue")

> plotSmear(MD.lrt_MD2vsMD1, xlab="log-counts per million (logCPM)", ylab="log-fold change (logFC)", main="Male flower stage 2 compare with Stage 1", smearWidth=0.5, pch=21, cex=0.4, deCol="red", col=MD.col_MD2vsMD1, ylim=c(-7,7), yaxs="i")

> abline(h=c(-1,1),col="dodgerblue")

导出需要的图像

> png("./output/smearPlot.maleflowerRNAseq_MD2vsMD1.png")

> plotSmear(MD.lrt_MD2vsMD1, xlab="log-counts per million (logCPM)", ylab="log-fold change (logFC)", main="Male flower Stage 2 compare with Stage 1", smearWidth=0.5, pch=21, cex=0.4, deCol="red", col=MD.col_MD2vsMD1, ylim=c(-7,7), yaxs="i")

> abline(h=c(-1,1),col="dodgerblue")

> dev.off()