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

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)

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()

## 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()

## 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()

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'

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

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)

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)

#~~~~~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)

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)