RNA-Seq(5):构建表达矩阵，归一化样本数据以及检测数据可信度并可视化

通过featurecounts计数后最后得到一个merge.txt文件，如下图：

image.png

但是我们并不能直接拿着数据去做后面的差异分析，得将数据归一化后才能进行下一步操作。那么问题来了，为什么得归一化后才能进行下一步操作，如何归一化。

1.为什么要归一化：RNA-Seq归一化问题 - 简书 (jianshu.com)

2.如何归一化（归一化步骤）

这里我采用的是使用DESeq2包归一化的两种方法获取归一化数值（当然你也可以采用其它办法归一化），有对归一化过程感兴趣的可以看归一化的算法，excel演示DESeq2归一化原理 - 简书 (jianshu.com)

1.归一化前数据准备，通过featurecounts计数后最后得到一个merge.txt文件，导入r,并查看数据

options(stringsAsFactors = FALSE)#stringsAsFactors = F意味着，“在读入数据时，遇到字符串之后，不将其转换为factors，仍然保留为字符串格式”。
raw_count <- read.table("E:/ZYH/R.project/rna-seq/lianxi1/merge.txt",sep="/t", header = T)

image.png

2.对数据进行处理

#先把第一列当作行名来处理
row.names(raw_count)<-make.names(raw_count[,1],TRUE)
#删除第一列
raw_count<-raw_count[,-1]
删除表达矩阵中和为0的行
raw_count= raw_count[which(rowSums(raw_count==0)==0),]
#筛选所有样本中表达量的和大于10的基因
raw_count = raw_count[rowSums(raw_count)>10,]

3.制作colData

condition <- factor(c(rep("control", 4), rep("treat", 4)), levels = c("control","treat"))
colData <- data.frame(row.names=colnames(raw_count), condition)
colData
             condition
X1_untreated   control
X2_untreated   control
X3_untreated   control
X4_untreated   control
X1_Dex           treat
X2_Dex           treat
X3_Dex           treat
X4_Dex           treat

4.构建dds对象，dds=DESeqDataSet Object

载入r包
install.packages("pacman")# pacman包如果已经安装则不需要这句
r包需要许多依赖包，可以一次性载入
pacman::p_load(GenomicRanges, matrixStats, S4Vectors, MatrixGenerics, Biobase, IRanges, BiocGenerics, GenomeInfoDb, SummarizedExperiment)#如果其中有包未安装，则该条命令会安装并加载它
> library(DESeq2)
构建dds对象
dds <- DESeqDataSetFromMatrix(raw_count, colData, design= ~ condition)
#查看dds对象  
dim(dds)  查看行列
assay(dds)  
assayNames(dds)  
colSums(assay(dds))  查看列名
rowRanges(dds)  
colData(dds)  

5.在正式分析前，我们需要对数据预筛，虽然在运行DESeq2前没必要对低reads基因进行预过滤;但为了减少dds对象占内存的空间，提高计算速度，这里做了简单过滤，保留至少有10个reads的行。

dds = dds[rowSums(counts(dds)) > 10 ,]
查看dds
dim(dds)第一次dim(dds)有57905行，这次有21032行，筛掉了三万多

6.使用DESeq2包中的DESeq（）函数获取归一化数值

看一下原来矩阵长啥样
> head(assay(dds),3)
                X1_untreated X2_untreated X3_untreated X4_untreated X1_Dex X2_Dex
ENSG00000223972            0            0            0            0      1      0
ENSG00000227232          117          159           81          104     70     77
ENSG00000243485            0            0            0            0      0      0
                X3_Dex X4_Dex
ENSG00000223972      0      0
ENSG00000227232    137     82
ENSG00000243485      0      0
=======================进行归一化==============
方法一：分步进行均一化，均一化后取对数
dds <- estimateSizeFactors(dds) 计算归一化系数sizeFactor，colData(dds)多了sizeFactor这一列，对测序深度和文库组成进行校正。

norm_dds <- counts(dds,normalized=T) 
均一化之后的数据
> head(norm_dds,3)
                X1_untreated X2_untreated X3_untreated X4_untreated    X1_Dex
ENSG00000223972       0.0000       0.0000      0.00000       0.0000  1.119579
ENSG00000227232     115.5627     134.6122     68.60638     112.8743 78.370552
ENSG00000243485       0.0000       0.0000      0.00000       0.0000  0.000000
                  X2_Dex   X3_Dex   X4_Dex
ENSG00000223972   0.0000  0.00000  0.00000
ENSG00000227232 114.5897 97.91823 86.42603
ENSG00000243485   0.0000  0.00000  0.00000
对归一化后的数据取log值,这样能降低数据的离散程度，方便后面的PCA和聚类分析
r.norm_dds = log2(norm_dds)

===============
方法二：使用DESeq2包中的rlog函数，它可以直接对数据归一化并自动进行对数转换，rlog转换可以得到与普通log2转换相似的结果，推荐使用这个方法,因为上一种方法得手动取对数，这种办法不用，哈哈，想偷懒的可以用这种
rld <- rlogTransformation(dds)

> head(assay(rld),3)
                X1_untreated X2_untreated X3_untreated X4_untreated    X1_Dex    X2_Dex
ENSG00000223972    -2.089725    -2.090478    -2.090475    -2.089265 -2.079351 -2.087724
ENSG00000227232     6.760110     6.892673     6.346646     6.739761  6.450597  6.749787
ENSG00000243485     0.000000     0.000000     0.000000     0.000000  0.000000  0.000000
                   X3_Dex    X4_Dex
ENSG00000223972 -2.091304 -2.089408
ENSG00000227232  6.623112  6.524621
ENSG00000243485  0.000000  0.000000
# 如果有需要，可以导出保存
> write.csv(assay(rld), file="rld.csv")

2.对归一化后的表达矩阵进行质控，样本间相关性分析

1.为什么要质控：拿到归一化数据后，我们不能直接去做差异分析，得先考虑这些问题：哪些样品相似，哪些不同?数据集变异的主要来源是什么?这符合实验设计的期望吗?所以得对计数数据进行QC检查，以帮助我们确保样本/重复看起来符合我们的实验设计。

为了探究样本之间的相关性，将使用两个方法主成分分析(PCA)和相关性分析并进行层次聚类分析方法执行样本级质量控制。这些方法允许我们检查重复样本彼此之间的相似程度(聚类)，并确保实验条件是数据变异的主要来源。
1.什么是PCA：笔记正在写

2.什么是相关性分析并进行层次聚类分析方法：
1.协方差与相关系数 - 简书 (jianshu.com)
2.层次聚类是将研究对象按照它们的相似性关系用树形图进行呈现。

2.R语言如何质控以及对质控结果可视化代码：1.主成分分析(PCA) 法 2.相关性分析并进行层次聚类分析方法

也可以根据这篇文章画PCA分析图R语言主成分分析(PCA)加“置信椭圆” - 简书 (jianshu.com)

library(ggplot2)
library("pheatmap")  
library("RColorBrewer")  

## 方法1--PCA-----------
作图方法一：DESeq自带的函数plotPCA，可以一步完成PAC分析并画图
plotPCA(rld, intgroup = "condition") 
作图方法二：使用factoextra包分析
###PCA分析，PCA分析很简单，只需要一个函数就可以搞定。prcomp（）就可以，然后导出分析结果。
library(factoextra)
data <- t(assay(rld))#PCA分析需要将表达矩阵转置 
PCA降维计算
> data.pca <- prcomp(data[ , which(apply(data, 2, var) != 0)], scale=TRUE)   #这是后续作图的文件
#输出特征向量
write.table(data.pca$rotation, file="PC.xls", quote=F, sep = "/t")
#输出新表
write.table(predict(data.pca), file="newTab.xls", quote=F, sep = "/t")
#输出PC比重
pca.sum=summary(data.pca)
write.table(pca.sum$importance, file="importance.xls", quote=F, sep = "/t")
设置分组
group=c(rep("control",4),rep("treat",4))
fviz_pca_ind(data.pca,
             col.ind=condition,
             mean.point=F,
             addEllipses = T, 
             legend.title="Groups",
             ellipse.type="confidence",
             ellipse.level=0.9,
             palette = c("#CC3333", "#339999"))+ Cell配色theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))


#方法2--相关性分析--------------------------  

# 利用绝对中位数mad或标准差sd统计学方法进行数据异常值检测
# 将表达量的绝对中位数mad从大到小排列取前500的结果
library(corrplot)
dat <- assay(rld)
tmp <- sort(apply(dat,1,mad),decreasing = T)[1:500]
exprSet <- dat[names(tmp),]
# 使用500个基因的表达量来做相关性图
dim(exprSet)
# 计算相关性,计算所有样本的两两相关值。
M <- cor(exprSet)
#查看相关性系数，看一下相关矩阵的列名和行名
> head(M,3)
             X1_untreated X2_untreated X3_untreated X4_untreated    X1_Dex    X2_Dex
X1_untreated    1.0000000    0.9200403    0.9451242    0.9455734 0.9452363 0.8451718
X2_untreated    0.9200403    1.0000000    0.9114139    0.9767463 0.9133338 0.9413762
X3_untreated    0.9451242    0.9114139    1.0000000    0.9110983 0.8850036 0.8175655
                X3_Dex    X4_Dex
X1_untreated 0.8604055 0.8520192
X2_untreated 0.8873660 0.9339959
X3_untreated 0.8849059 0.8083923
# 相关性热图，绘图 
corrplot(M,order = "AOE", addCoef.col = "white")
pheatmap(M, annotation =colData,display_numbers = T)

PCA图

image.png

层次聚类相关性热图，相关性还不错

image.png

横纵坐标为各个样本，将这些样本两两进行比较，得出皮尔森相关系数。图中可以看出，样本自己与自己比较，皮尔森相关系数是1，与其他样本进行比较，皮尔森相关系数均大于0.9表明各样本之间的正相关性比较好，这也说明送样样本的生物学重复很好，数据是可信的。