RNA-seq入门实战（四）：差异分析前的准备——数据检查

连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中几个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！

他前面的分享是：

下面是他对我们b站转录组视频课程的详细笔记

本节概览：

1. 数据预处理, 进行样本间具有可比性

2. boxplot查看样本的基因整体表达情况

3. 查看不同分组的聚类情况：样本hclust 图、距离热图、PCA图、差异基因热图、相关性热图

承接上节 RNA-seq入门实战（三）：在R里面整理表达量counts矩阵 和  RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon  

在进行差异分析前需要进行数据检查，保证我们的下游分析是有意义的。

以下展示了样本hclust 图、距离热图、PCA图、前500差异性大的基因热图、相关性热图（选取了500高表达基因，防止低表达基因造成的干扰），确定我们不同样本间确实是有差异的。这些图并不是全都是必须的，它们全都是为了说明一个问题：我们的不同分组间确实存在差异。

rm(list = ls())  
options(stringsAsFactors = F)
library(FactoMineR)
library(factoextra)  
library(tidyverse) # ggplot2 stringer dplyr tidyr readr purrr  tibble forcats
library(pheatmap)
library(DESeq2)
library(RColorBrewer)

#### 载入数据 设置目录
setwd("C:/Users/Lenovo/Desktop/test")
load(file = '1.counts.Rdata')
dir.create("2.check")
setwd("2.check")

#### 数据预处理  # (任选以下一种作为dat即可，主要是进行样本间归一化，使得样本具有可比性)
#dat <- as.data.frame(log2(edgeR::cpm(counts)+1))  #简单归一化 CPM:Counts per million
dat <- log2(tpm+1)
#DESeq2_normalize   rld 
if (F) { 
  dds <- DESeqDataSetFromMatrix(countData = counts,
                                colData = gl,
                                design = ~ group_list)
  rld <- rlog(dds, blind=FALSE)
  write.table(assay(rld),  file="Deseq2_rld.txt", sep="\t", quote=F, col.names=NA)
  dat <- as.data.frame(assay(rld))     
}

###boxplot 查看样本的基因整体表达情况
boxplot(dat,col=color, ylab="dat", main=" normalized data ",
        outline = F, notch = F)
dev.off()

###################### hclust and Heatmap of the sample-to-sample distances ###########################
sampleDists <- dist(t(dat))   #dist默认计算矩阵行与行的距离， 因此需要转置
sampleDistMatrix <- as.matrix(sampleDists)  
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)  #选取热图的颜色
p0 <- pheatmap::pheatmap(sampleDistMatrix,
                         fontsize=7,
                         clustering_distance_rows=sampleDists,
                         clustering_distance_cols=sampleDists,
                         angle_col=45,
                         col=colors)
ggsave(p0,filename = 'check_dist.pdf',width = 7.5,height =6)
dev.off()

pdf("check_hclust.pdf")
plot(hclust(sampleDists))
dev.off()

################################# PCA检测 #####################################
#PCA查看实验和对照组情况
dat.pca <- PCA(t(dat) , graph = F) 
pca <- fviz_pca_ind(dat.pca,
                    title = "Principal Component Analysis",
                    legend.title = "Groups",
                    geom.ind = c("point", "text"), 
                    pointsize = 1.5,
                    labelsize = 4,
                    col.ind = group_list, # 分组上色
                    axes.linetype=NA,  # remove axeslines
                    mean.point=F#去除分组中心点
                    ) + 
  coord_fixed(ratio = 1)+  #坐标轴的纵横比
  xlab(paste0("PC1 (",round(percentVar[1,'variance.percent'],1),"%)")) +
  ylab(paste0("PC2 (",round(percentVar[2,'variance.percent'],1),"%)"))
#若用 rld 数据，还可使用DESeq2自带函数 
#pca <- plotPCA(rld, ntop = 500, intgroup=c("group_list"))
ggsave(pca, filename = 'check_PCA.pdf',width = 7.5,height =6)
dev.off()

####################### heatmap检测——取500差异大的基因 ##########################################
cg <- names(tail(sort(apply(dat,1,sd)),500)) #取每一行的方差，从小到大排序，取最大的500个
p1 <- pheatmap::pheatmap(n,show_colnames =T,show_rownames = F,
                         fontsize=7,
                         legend_breaks = -3:3,
                         #scale = "row",
                         angle_col=45,
                         annotation_col=gl) 
}
ggsave(p1,filename = 'check_heatmap_top500_sd.pdf',width = 7.5,height =6)
dev.off()

#######################样本相关性检测————取500高表达基因##################################
dat_500 <- dat[names(sort(apply(dat,1,mad),decreasing = T)[1:500]),]#取高表达量前500基因
M <- cor(dat_500)

p2 <-pheatmap::pheatmap(M,
                   show_rownames = T,
                   angle_col=45,
                   fontsize=7,
                   annotation_col = gl ) }
ggsave(p2,filename = 'check_cor_top500.pdf',width = 7.5,height =6)

出图如下：

image.png

image.png

image.png

---

从以上各图可以看出，我们进行归一化后的数据在各样本间分布一致，因此各样本间是可比较的。各种聚类可视化图也可以明显看出我们的两个分组之间确实存在有很大的差异，组间样品是分开的，组内是聚在一起的，因此我们就可以自信地进行下一步的差异分析啦。

这就是生信技能树一直提到的三张图，生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

• 左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
• 中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
• 右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。

参考资料

Analyzing RNA-seq data with DESeq2 (bioconductor.org)

GitHub - jmzeng1314/GEO （强烈推荐学习）
本实战教程基于以下生信技能树分享的视频

【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili

【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年5~6月场，你的生物信息学第一课