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

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)