转录组-样品表达总体分布及质控可视化

转录组-样品表达总体分布及质控可视化

在拿到表达矩阵时我们常常会对其基因表达的总体分布（可选），以及质量控制进行可视化（必须）。这里总结记录相关代码。

1 总体分布

表达矩阵总体分布常常用3种图可视化：

代码语言：r复制

rm(list = ls())
options(stringsAsFactors = F)

# 加载包，设置绘图参数
library(ggplot2)
library(ggsci)
library(tidyverse)

mythe <- theme_bw()   theme(panel.grid.major=element_blank(),
                            panel.grid.minor=element_blank())


# 加载原始表达的数据
load(file = 'symbol_matrix.Rdata')
dat = log2(edgeR::cpm(symbol_matrix) 1)

1.1 箱式图

代码语言：r复制

data <- dat %>% 
  as.data.frame() %>% 
  pivot_longer(cols = everything(), names_to = "sample",values_to = "expression")

head(data)

mythe <- theme_bw()   theme(panel.grid.major=element_blank(),
                            panel.grid.minor=element_blank())
p <- ggplot(data = data, aes(x=sample,y=expression,fill=sample))
p1 <- p   geom_boxplot()   
  mythe   theme(axis.text.x = element_text(angle = 90))   
  xlab(NULL)   ylab("log2(CPM 1)")   scale_fill_lancet()

p1
# 保存图片
png(file = "1.sample_boxplot.png",width = 800, height = 900,res=150)
print(p1)
dev.off()

1.2 小提琴图

代码语言：r复制

## 2.样本表达总体分布-小提琴图
p2 <- p   geom_violin()    mythe  
  theme(axis.text = element_text(size = 12),
        axis.text.x = element_text(angle = 90))   
  xlab(NULL)   ylab("log2(CPM 1)") scale_fill_lancet()
p2

# 保存图片
png(file = "1.sample_violin.png",width = 800, height = 900,res=150)
print(p2)
dev.off()

1.3 概率密度分布图

代码语言：r复制

## 3.样本表达总体分布-概率密度分布图
m <- ggplot(data=data, aes(x=expression)) 
p3 <- m    geom_density(aes(fill=sample, colour=sample),alpha = 0.1)   
  xlab("log2(CPM 1)")   mythe   scale_fill_npg()
p3

# 保存图片
png(file = "1.sample_density.png",width = 800, height = 700, res=150)
print(p3)
dev.off()

2 质量控制

对表达矩阵质量控制可视化是转录组标准分析流程中必备的一步，通常包含3张图：

代码语言：r复制

rm(list = ls()) 
options(stringsAsFactors = F)
library(ggplot2)
library(ggstatsplot)
library(ggsci)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
load(file = 'symbol_matrix.Rdata') 
colnames(symbol_matrix)
symbol_matrix[1:4,1:4] ## 基因名字的样品，矩阵
dat = log2(edgeR::cpm(symbol_matrix) 1)
pro = 'test'

2.1 单独的基因可视化

这里可挑选感兴趣的基因进行可视化，这里的target_gene以表达矩阵中的第一个基因为例。

代码语言：r复制

bp=function(g){         #定义一个函数g，函数为{}里的内容
  library(ggpubr)
  df=data.frame(expression = g,group = group_list)
  p <- ggboxplot(df, x = "group", y = "expression",
                 color = "group", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p   stat_compare_means()   ggtitle(target_gene) 
    theme_bw()
}
target_gene=head(rownames(dat),1)
p1 <-  bp(dat[target_gene,])
p1

2.2 样本的PCA图

代码语言：r复制

exp=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame 

library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
#～～～主成分分析图p2～～～
dat.pca <- PCA(exp , graph = FALSE)
this_title <- paste0(pro,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
                   geom.ind = "point",  #c( "point", "text" ), # show points only (nbut not "text")
                   col.ind = group_list, # color by groups
                   palette = "Dark2",
                   addEllipses = TRUE, # Concentration ellipses
                   legend.title = "Groups") 
  ggtitle(this_title)  
  theme(plot.title = element_text(size=12,hjust = 0.5))

p2
ggsave('qc_pca.pdf',width = 5,height = 5)

2.3 方差最大的前1000个基因热图

代码语言：r复制

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2 
n[n< -2]= -2
head(n)
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(Group=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
         annotation_col=ac)

#～～～top1000热图p3～～～
p3 <- pheatmap::pheatmap(n,
                         show_colnames =F,
                         show_rownames = F,
                         main = pro,
                         annotation_col=ac,
                         breaks = seq(-3,3,length.out = 100))
p3
ggsave('qc_top1000_pheatmap.pdf',plot = p3,width = 5,height = 5)

2.4 组间差异与组内差异的比较

组内的样本的相似性应该是要高于组间的

代码语言：r复制

#~~~~~corplot~~~~~
cg=names(tail(sort(apply(dat,1,sd)),1000))
exprSet=dat[cg,] # dat
pheatmap::pheatmap(cor(exprSet)) 
# 组内的样本的相似性应该是要高于组间的！
colD=data.frame(Group=group_list)
rownames(colD)=colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),
                   annotation_col = colD,
                   show_rownames = F)
#～～～样品相关性热图p4～～～
p4 <- pheatmap::pheatmap(cor(exprSet),
                         annotation_col = colD,
                         display_numbers = T,
                         show_rownames = F,
                         show_colnames =F,
                         main = pro,
)
p4
ggsave('qc_top1000_corplot.pdf',plot = p4,width = 5,height = 5)

可以拼图总体展示

代码语言：r复制

p_check <- (p1 p2 as.ggplot(p3) as.ggplot(p4) 
              plot_layout(widths = c(2.5,1.5,4.5,4.5)))
p_check
path='./'
ggsave(plot = p_check,# path = file.path(path,'plot'),
       'check_plot.pdf',width = 18,height = 4)

数据挖掘

0 人点赞