可视化

rm(list = ls())
load("TCGA-CHOL_gdc.Rdata")
load("TCGA-CHOL_DEG.Rdata")
if(!require(tinyarray))devtools::install_local("tinyarray-master.zip",upgrade = F)  #画图代码包裹成函数放在了tinyarray包里面
library(ggplot2)
library(tinyarray)
exp[1:4,1:4]  #count矩阵仅用于做差异分析，画图需要对其标准化

cpm 去除文库大小的影响

dat = log2(cpm(exp)+1)

pca图

pca.plot = draw_pca(dat,Group);pca.plot
save(pca.plot,file = paste0(cancer_type,"_pcaplot.Rdata"))   

画热图

 #取出差异基因
cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change !="NOT"] 
cg2 = rownames(edgeR_DEG)[edgeR_DEG$change !="NOT"]
cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change !="NOT"]
![截屏2021-09-19 下午7.30.25.png](https://img.haomeiwen.com/i25840690/841c14763730bd77.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

h1 = draw_heatmap(dat[cg1,],Group,n_cutoff = 2)
 #提供数据，分组，n_cutoff设置颜色跨度，可根据基因表达量主要分布范围设置-2-2
h2 = draw_heatmap(dat[cg2,],Group,n_cutoff = 2)
h3 = draw_heatmap(dat[cg3,],Group,n_cutoff = 2)

火山图

#计算阈值的函数，也可自己设定
m2d = function(x){
  mean(abs(x))+2*sd(abs(x))  
}  
#每个R包算出来的基因表达不一样，算出来的域值也不一样

v1 = draw_volcano(DESeq2_DEG,pkg = 1,logFC_cutoff = m2d(DESeq2_DEG$log2FoldChange),symmetry = T)  

 #pkg =  1,2,3,4 means用"DESeq2","edgeR","limma(voom)","limma"不同包处理得到的数据，p value默认设置=o.o5 ；
#因为取名有重叠，后面再分别带入m2d这个函数计算不同包的阈值;前后画图阈值要一致
#symmetry = T调节0位置居中对称
v2 = draw_volcano(edgeR_DEG,pkg = 2,logFC_cutoff = m2d(edgeR_DEG$logFC))
v3 = draw_volcano(limma_voom_DEG,pkg = 3,logFC_cutoff = m2d(limma_voom_DEG$logFC))

截屏2021-09-19 下午7.30.25.png
截屏2021-09-19 下午7.32.49.png

拼图

patchwork包只能用来拼ggplot2画的图

library(patchwork)  
(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect') &theme(legend.position = "none")

ggsave(paste0(cancer_type,"_heat_vo.png"),width = 15,height = 10)

三大R包差异基因对比（同up同down）

三个R包得到的差异基因都出自同一表达矩阵，差别在于算法不同

#用函数挑选出表达矩阵上下调基因行
rm(list = ls())
load("TCGA-CHOL_gdc.Rdata")
load("TCGA-CHOL_DEG.Rdata")
load("TCGA-CHOL_pcaplot.Rdata")
UP=function(df){
  rownames(df)[df$change=="UP"]
}   
DOWN=function(df){
  rownames(df)[df$change=="DOWN"]
}
#对3个R包的上下调基因取交集
up = intersect(intersect(UP(DESeq2_DEG),UP(edgeR_DEG)),UP(limma_voom_DEG))  
down = intersect(intersect(DOWN(DESeq2_DEG),DOWN(edgeR_DEG)),DOWN(limma_voom_DEG))
dat = log2(cpm(exp)+1)

根据公认数据做热图

hp = draw_heatmap(dat[c(up,down),],Group,n_cutoff = 2) 

上调、下调基因分别画维恩图

#要对哪3个向量画韦恩图就把它组织成三个有名字的列表
up_genes = list(Deseq2 = UP(DESeq2_DEG),
          edgeR = UP(edgeR_DEG),
          limma = UP(limma_voom_DEG))

down_genes = list(Deseq2 = DOWN(DESeq2_DEG),
          edgeR = DOWN(edgeR_DEG),
          limma = DOWN(limma_voom_DEG))

up.plot <- draw_venn(up_genes,"UPgene")
down.plot <- draw_venn(down_genes,"DOWNgene")

维恩图拼图

library(patchwork)
#up.plot + down.plot
# 拼图
pca.plot + hp+up.plot +down.plot+ plot_layout(guides = "collect")
ggsave(paste0(cancer_type,"_heat_ve_pca.png"),width = 15,height = 10)

image.png