本节概览：
- 数据预处理, 进行样本间具有可比性
- boxplot查看样本的基因整体表达情况
- 查看不同分组的聚类情况：样本hclust 图、距离热图、PCA图、差异基因热图、相关性热图
承接上节RNA-seq入门实战（三）：从featureCounts与Salmon输出文件获取counts矩阵
在进行差异分析前需要进行数据检查，保证我们的下游分析是有意义的。
以下展示了样本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)
出图如下：
从以上各图可以看出，我们进行归一化后的数据在各样本间分布一致，因此各样本间是可比较的。各种聚类可视化图也可以明显看出我们的两个分组之间确实存在有很大的差异，组间样品是分开的，组内是聚在一起的，因此我们就可以自信地进行下一步的差异分析啦
参考资料
Analyzing RNA-seq data with DESeq2 (bioconductor.org)
GitHub - jmzeng1314/GEO （强烈推荐学习）
本实战教程基于以下生信技能树分享的视频
【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili
RNA-seq实战系列文章：
RNA-seq入门实战（零）：RNA-seq流程前的准备——Linux与R的环境创建
RNA-seq入门实战（一）：上游数据下载、格式转化和质控清洗
RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon
RNA-seq入门实战（三）：从featureCounts与Salmon输出文件获取counts矩阵
RNA-seq入门实战（四）：差异分析前的准备——数据检查
RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较
RNA-seq入门实战（六）：GO、KEGG富集分析与enrichplot超全可视化攻略
RNA-seq入门实战（七）：GSEA——基因集富集分析
RNA-seq入门实战（八）：GSVA——基因集变异分析
RNA-seq入门实战（九）：PPI蛋白互作网络构建（上）——STRING数据库的使用
RNA-seq入门实战（十）：PPI蛋白互作网络构建（下）——Cytoscape软件的使用
RNA-seq入门实战（十一）：WGCNA加权基因共表达网络分析——关联基因模块与表型
