miRNA seq差异表达分析练习（三）——差异结果可视化

上一篇已经完成了差异分析的整个过程，先看一下结果

结果表明，在p<0.001条件下，有105个上调表达的miRNA，26个下调表达的miRNA

此外，还有996个低表达的miRNA。

绘制MA-plot

plotMA(res,ylim=c(-3,3))

 

图中绘制了log2foldchange在(-3,+3)以内的基因，即改变程度在(0.125,8)倍的基因，对于不在此范围内的基因，用三角形的标志画在边界线上。从图像可以看出，DESeq2的负二项分布模型找出来的差异表达基因，大部分都是reads数目多（测序深度高），且表达量差异很大的基因。

绘制火山图

# plot作图
plot(res$log2FoldChange,res$padj)

说实话这图太丑了，无法忍受！而且也很难解释图像的含义，于是用ggplot2重新做了一个。

library(ggplot2)
data_res <- as.data.frame(res)
data_res$color <- ifelse(data_res$padj<0.05&abs(data_res$log2FoldChange)>2,'LFC>2,p<0.05',ifelse(abs(data_res$log2FoldChange)>2,'LFC>2,p>0.05','LFC<2,p>0.05'))
ggplot(data_res,aes(log2FoldChange,padj,col=factor(color))) +
  geom_point(alpha=0.8, size = 1) +
  theme_bw(base_size = 15) +
  theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())+
  geom_vline(xintercept=c(-2,2) ,linetype=4 ) +
  scale_color_manual(name = "", values = c("red", "green", "black"), limits = c("LFC>2,p<0.05", "LFC>2,p>0.05", "LFC<2,p>0.05")) 

这下好看多了！也可以考虑将图例移到图像内的右上角

解释一下，红点表示差异倍数大于2且矫正p值小于0.05，蓝点表示差异倍数大于2但矫正p值没有小于0.05，黑点就是剩下的那些差异倍数不大，矫正p值却较大的miRNA。图中可以明显的看出上调的miRNA（右边）要多余下调的miRNA（左边）。

提取差异miRNA并保存

diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#write.csv(diff_gene_deseq2,file= "D:\\用户\\桌面\\DEG_treat_vs_control.csv")

这里，我仅保存了差异倍数大于2且矫正p值小于0.05，即color列为“red"的miRNA。

绘制热图

选择差异最显著的前20个miRNA绘制热图

# 选择差异最显著的前20个miRNA绘制热图
library('pheatmap')
diff_miRNA <- row.names(data_res)[1:20] # 提取前20的差异miRNA
# 从mydata中提取差异最显著的前20个miRNA的表达数据
myvars <- row.names(mydata) %in% diff_miRNA
heatmap_data <- mydata[myvars,]
mymatrix <- as.matrix(heatmap_data)
annotation_col = data.frame(
  CellType = factor(c(rep("N",20),rep("T",20),rep("P",20))), 
  condition
)
rownames(annotation_col) <- c(paste('N-',1:20,sep = ''),paste('T-',1:20,sep = ''),paste('P-',1:20,sep = ''))
pheatmap(heatmap_data,scale = 'row', annotation_col = annotation_col, show_rownames = FALSE)

看最上面的聚类轴,正常组织和癌变组织的差异还是比较明显的.

主成分分析

# 所有差异miRNA主成分分析
miRNA <- miRNA <- diff_gene_deseq2@rownames
vars <- row.names(mydata) %in% miRNA
pca_data <- mydata[vars,]
pca_matrix <- as.matrix(pca_data)
dds1 <- DESeqDataSetFromMatrix(pca_matrix,DataFrame(coldata), design = ~org)
dds1 = DESeq(dds1)
rld <- rlog(dds1)
plotPCA(rld,intgroup=c('org'))

散点图

选择差异最显著的miRNA绘制散点图观察在三个组织的表达差异

# 选择差异最显著的miRNA绘制散点图观察在三个组织的表达差异
point_data <- data.frame(counts=as.numeric(mydata[diff_miRNA[1],]))
point_data$org <- c(rep('N',20),rep('P',20),rep('T',20))
ggplot(point_data,aes(org,counts))+geom_point()

 

不愧是差异最显著的miRNA,在正常组织中几乎不表达，而在癌变组织中表达丰富。

关于miRNA差异表达分析的练习就到这了，有什么问题留言交流吧！