RNA 7. SCI 文章中的基因表达——主成分分析 (PCA)

2023-11-06 11:43:58

在RNA-seq中，主成分分析（PCA）是最常见的多元数据分析类型之一，这期主要介绍一下利用已有的表达差异数据如何分析，别着急，见下文。

1. 前言

1. 相关背景

在RNA-seq中，主成分分析（PCA）是最常见的多元数据分析类型之一。基因表达定量后获得了各样本中所有基因的表达值信息，随后我们通常会期望比较样本之间在基因表达值的整体相似性或者差异程度。基因数量成千上万，肯定不能对每个基因的表达都作个比较，这时候就要用到"降维"算法，PCA分析因此派上用场。PCA设法将N维（N=基因数量）的表达矩阵降维成只有少数几个主成分的形式，这几个主成分就可以代表所有基因的整体表达格局，进而据此描述样本差异。

2. 数据降维

降维就是一种对高维度特征数据预处理方法。降维是将高维度的数据保留下最重要的一些特征，去除噪声和不重要的特征，从而实现提升数据处理速度的目的。在实际的生产和应用中，降维在一定的信息损失范围内，可以为我们节省大量的时间和成本。降维也成为应用非常广泛的数据预处理方法。

降维具有如下一些优点：

使得数据集更易使用。
降低算法的计算开销。
去除噪声。
使得结果容易理解。

PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。以此类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。降维的算法有很多，比如奇异值分解(SVD)、主成分分析(PCA)、因子分析(FA)、独立成分分析(ICA)。

2. 实例解析

数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，我们使用的是edgeR 软件包计算出来的差异表达结果，合并了原始的 Count 数据，如下：

1. 数据读取

之前有保存最后的表达差异的结果，直接读取数据即可，这里面我们是比较癌和癌旁的表达差异，数据获得，如下：

group<-read.table("DEG-group.xls",sep="\t",check.names=F,header = T)
DEG=read.table("DEG-resdata.xls",sep="\t",check.names=F,header = T)

PCAmat<-DEG[,8:ncol(DEG)]
rownames(PCAmat)=DEG[,1]
PCAmat[1:5,1:3]


##                 TCGA-3L-AA1B-01A-11R-A37K-07
## ENSG00000142959                           20
## ENSG00000163815                          175
## ENSG00000107611                           49
## ENSG00000162461                           55
## ENSG00000163959                          153
##                 TCGA-4N-A93T-01A-11R-A37K-07
## ENSG00000142959                           15
## ENSG00000163815                          108
## ENSG00000107611                           13
## ENSG00000162461                           89
## ENSG00000163959                          259
##                 TCGA-4T-AA8H-01A-11R-A41B-07
## ENSG00000142959                           49
## ENSG00000163815                           59
## ENSG00000107611                            6
## ENSG00000162461                           48
## ENSG00000163959                          102

将基因表达值矩阵作个转置，使行为样本，列为基因，如下：

PCAmat_t<-t(PCAmat)
PCAmat_t[1:5,1:3]


##                              ENSG00000142959 ENSG00000163815
## TCGA-3L-AA1B-01A-11R-A37K-07              20             175
## TCGA-4N-A93T-01A-11R-A37K-07              15             108
## TCGA-4T-AA8H-01A-11R-A41B-07              49              59
## TCGA-5M-AAT4-01A-11R-A41B-07              16              60
## TCGA-5M-AAT5-01A-21R-A41B-07              24              60
##                              ENSG00000107611
## TCGA-3L-AA1B-01A-11R-A37K-07              49
## TCGA-4N-A93T-01A-11R-A37K-07              13
## TCGA-4T-AA8H-01A-11R-A41B-07               6
## TCGA-5M-AAT4-01A-11R-A41B-07              24
## TCGA-5M-AAT5-01A-21R-A41B-07              10

2. PCA {FactoMineR}

软件包 FactoMineR 安装并加载，如下：

if(!require(FactoMineR)){
  install.packages("FactoMineR")
}

library(FactoMineR)

在这个 FactoMineR 中提出的方法是探索性的多元方法，如主成分分析，对应分析或聚类。我们通过这个包，实现 PCA 分析，如下：

#############PCA {FactoMineR}
library(ggplot2)
library(ggrepel)
#样本中基因表达值的 PCA 分析
gene.pca <- PCA(PCAmat_t, ncp = 2, scale.unit = TRUE, graph = FALSE)

对结果进行数据处理，如下：

#提取样本在 PCA 前两轴中的坐标
pca_sample <- data.frame(gene.pca$ind$coord[ ,1:2])
pca_sample$Sample=row.names(pca_sample)
#提取 PCA 前两轴的贡献度
pca_eig1 <- round(gene.pca$eig[1,2], 2)
pca_eig2 <- round(gene.pca$eig[2,2],2 )
pca_sample <- merge(pca_sample, group,by="Sample")
head(pca_sample)


##                         Sample      Dim.1       Dim.2 Group
## 1 TCGA-3L-AA1B-01A-11R-A37K-07  1.1786551  -0.7589326    TP
## 2 TCGA-4N-A93T-01A-11R-A37K-07 -1.9561110  -5.4213738    TP
## 3 TCGA-4T-AA8H-01A-11R-A41B-07 -5.1936884  -7.9477300    TP
## 4 TCGA-5M-AAT4-01A-11R-A41B-07  1.7605422 -11.8875367    TP
## 5 TCGA-5M-AAT5-01A-21R-A41B-07  0.9007755 -10.0044126    TP
## 6 TCGA-5M-AAT6-01A-11R-A41B-07 -1.8586282  -3.0063580    TP

利用ggplot2 绘制二维散点图，我们看到两个组癌和癌旁基本上是能分开的，这里我们选择的是差异较大的基因有关系，如下：

library(ggplot2)
p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
  geom_point(aes(color = Group), size = 2) +  #根据样本坐标绘制二维散点图
  scale_color_manual(values = c('orange', 'purple')) +  #自定义颜色
  theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.key = element_rect(fill = 'transparent')) +  #去除背景和网格线
  labs(x =  paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '')  #将 PCA 轴贡献度添加到坐标轴标题中

p

有的时候也需要标识样本名称，使用 ggplot2 的拓展包 ggrepel 来完成，如下：

library(ggrepel)
p <- p + 
  geom_text_repel(aes(label = Sample), size = 1, show.legend = FALSE, 
                  box.padding = unit(0.2, 'lines'))

p

添加 95% 置信椭圆，可用于表示对象分类，但只能适用于各组样本数大于 5 个的情况，我们这里样本比较多，做个例子而已，如下：

p + stat_ellipse(aes(color = Group), level = 0.95, show.legend = FALSE)

添加阴影区域，如下：

p + stat_ellipse(aes(fill = Group), geom = 'polygon', level = 0.95, alpha = 0.3, show.legend = FALSE) +
  scale_fill_manual(values = c('orange', 'purple'))

3. plotPCA {DESeq2}

当我们使用的数据是表达量的时候，我们可以首选利用 DESeq2 软件包中的内置函数 plotPCA 来绘制主成分分析图，非常方便，当然我们这里使用差异表达挑选出来的差异表达基因，在做主成分分析时能够更好的区分癌和癌旁组织。

首先需要构造 dds 对象，如下：

#####plotPCA {DESeq2}
library(DESeq2)
condition<-factor(group$Group)
dds <- DESeqDataSetFromMatrix(PCAmat, 
                              DataFrame(condition), 
                              design= ~ condition )
dds


## class: DESeqDataSet 
## dim: 4128 519 
## metadata(1): version
## assays(1): counts
## rownames(4128): ENSG00000142959 ENSG00000163815 ...
##   ENSG00000230184 ENSG00000244217
## rowData names(0):
## colnames(519): TCGA-3L-AA1B-01A-11R-A37K-07
##   TCGA-4N-A93T-01A-11R-A37K-07 ... TCGA-AZ-6605-11A-01R-1839-07
##   TCGA-F4-6704-11A-01R-1839-07
## colData names(1): condition

数据进行标准化处理，如下：

rld <- vst(dds)   #DEseq2自己的方法标准化数据`
exprSet_new=assay(rld)   #提取DEseq2标准化后的数据
str(exprSet_new)


##  num [1:4128, 1:519] 3.87 6.69 4.96 5.11 6.5 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:4128] "ENSG00000142959" "ENSG00000163815" "ENSG00000107611" "ENSG00000162461" ...
##   ..$ : chr [1:519] "TCGA-3L-AA1B-01A-11R-A37K-07" "TCGA-4N-A93T-01A-11R-A37K-07" "TCGA-4T-AA8H-01A-11R-A41B-07" "TCGA-5M-AAT4-01A-11R-A41B-07" ...

主程序绘制 PCA 图形，如下：

plotPCA(rld, intgroup=c('condition'))+
  theme_bw()

怎么样，看懂了没？没看懂也没关系，下期就出视频教程了，也可以找我代劳，保您满意！

关注公众号，别走开，精彩下期继续！

References:

Le, S., Josse, J. & Husson, F. (2008). FactoMineR: An R Package for Multivariate Analysis. Journal of Statistical Software. 25(1). pp. 1-18.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139-140.
Agarwal A, Koppstein D, Rozowsky J, et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC Genomics. 2010;11:383.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

43篇原创内容

公众号：桓峰基因

码农公寓