生物本身就是一个复杂的体系，我们必须使用多维数据来对之进行准确描述，而这一点在生物信息学领域体现的尤为突出，以经典的转录组数据为例：
- bulk RNA-seq作为一种经典的基因相对表达量定量方法广泛应用于各种研究中，如果我们在一个课题中对4个样品进行了bulk RNA-seq测序，那么我们将得到一个4*N维的矩阵，这里N就是全基因组范围内能转录的基因。
- scRNA-seq是近十年来蓬勃发展的一种分辨率更高的RNA-seq，如果我们进行测序的组织样本中含有10000个细胞，那么我们在最后将得到一个10000*N维的矩阵。
数据纷繁复杂，我们不可能将所有的数据都进行分析，无论是从时间和空间上都是不现实的，因此在生物信息学（尤其是NGS组学）中数据降维几乎是一个无法逃避的主题。
那么究竟什么是数据降维呢？简单来说就是只保留那些在不同的组别（对bulk RNA-seq而言就是不同的样本，对scRNA-seq而言就是不同的细胞）中具有较大差异的特征，而舍弃掉那些非常“均一”的特征，所谓“拨云见日”。
数据降维究竟有什么实际应用呢？一个最简单的应用就是找到一个群体中相似的个体，而这一点也被大家广泛地应用到NGS测序数据质量评估中：一般而言，我们会给同一个生物学条件下准备多个生物学重复，我们更期待不同样品的组内差异小于组间差异，只有这样才能说明我们的数据可重复性好，实验条件更加稳定。
1. 常见数据降维方法
生物学中的数据降维方法都是从机器学习领域借鉴过来的，目前常用的主要包括线性降维的PCA（principal component analysis，主成分分析）、非线性降维的tSNE（t-Distribution Stochastic Neighbour Embedding）与UMAP（Uniform Manifold Approximation and Projection for Dimension Reduction）等。
- PCA：
PCA的主要思想是找到k组线性组合，将n维特征投影到k维上，使得新的k维特征彼此正交，并且可以尽可能多的解释数据中的variation。这样得到新的k维特征就被称为k个主成分。k个主成分实际上对应着我们从输入数据估计出的协方差矩阵最大的k个特征值对应的特征向量。PCA在在生物信息学分析中，是样本量不多的bulk RNA-seq转录组的可视化最常用的方法。
为了便于让大家更好的简单理解什么是PCA，在这里我用两张PPT做简单说明：
- tSNE与UMAP
t-SNE与UMAP主要的优势就是保持非线性的局部结构的能力，比较适用于样本量较大，样本之间相互关系比较复杂的数据，所以在单细胞转录组的可视化中有大量应用。需要注意的是，这两种将为的方法通常只应用于数据的可视化。
2. 代码实现
在本部分我将从数据的预处理到降维分析可视化的过程用R语言进行展示。
-
数据预处理
在进行降维分析之前，首先必须要做的一个工作就是数据的预处理，其中非常关键的一步就是scale（有时可以翻译为“归一化”）。简单来说，未经处理的数据可能具有非常大的变异度，这可能会对后续的分析产生影响，具体来说，以RNA-seq为例：不同的基因具有不同的表达水平，有些基因表达水平很高，有的很低；同一个基因在不同的生物学条件下表达的差异也不同，有的“岿然不动”，有的“闻风而动”，甚至具有非常大的差异。用鲁志老师的tutorial来说，对涉及距离计算的模型而言，在实际应用中，不同特征的数值大小不同，会在求距离时造成很大的影响。例如，一个基因的表达量很高，另一个基因的表达量很低，如果不进行数据的归一化，表达量高的基因就会主导距离的计算。
在R语言中我们常用的就是使用z-score里进行归一化：
这里我们使用鸢尾花数据集来进行演示，该数据集为R包datasets的一个内置数据集，记录了150朵鸢尾花不同的特征及其种属，一般可以直接加载使用即可：
head(iris)
# Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#1 5.1 3.5 1.4 0.2 setosa
#2 4.9 3.0 1.4 0.2 setosa
#3 4.7 3.2 1.3 0.2 setosa
#4 4.6 3.1 1.5 0.2 setosa
#5 5.0 3.6 1.4 0.2 setosa
#6 5.4 3.9 1.7 0.4 setosa
利用R语言的scale()函数进行归一化：
iris.scale <- scale(as.matrix(iris[, 1:4]),
scale = TRUE,
center = TRUE)
-
数据降维
1. PCA
首先是PCA，这里介绍两种方法，分别基于stats包的prcomp()函数与FactoMineR包的PCA()函数。
#use prcomp()
pca.res <- prcomp(iris.scale, scale. = FALSE, center = FALSE)
#set FALSE, because we have already done
plot.data <- cbind(as.data.frame(pca.res$x[, 1:2]),
Species = iris$Species)
library(ggplot2)
ggplot(data = plot.data, aes(PC1, PC2)) +
geom_point(aes(color = Species)) +
theme_test()
我们可以看到同一个物种的不同个体基本都被聚在了一起，而不同物种之间“泾渭分明”。（versicolor和virginica未被很好分开是因为这两个物种本身就具有相似的特征）
#use PCA()
#install.packages('factoextra')
#install.packages('FactoMineR')
library(factoextra)
library(FactoMineR)
pca.res <- PCA(iris.scale,
scale.unit = TRUE,
graph = FALSE)
fviz_pca_ind(pca.res,
geom = 'point',
habillage = iris$Species,
addEllipses = TRUE) +
theme_test()
结果类似，不过看起来好像PC2反过来了。可以看到这里有两个数字：73%和22.9%，你可以简单理解为每个主成分能够解释的数据变异度比例，根据前面的PPT，理论上来讲PC1解释度最高，PC2次之……所以我们尝试看一下PC3和PC4对鸢尾花个体的区分度如何：
fviz_pca_ind(pca.res,
axes = c(3, 4),
geom = 'point',
habillage = iris$Species,
addEllipses = TRUE) +
theme_test()
可以说基本没有分开，原因就在于这两个主成分没有很好的解释整个数据集的变异度（3.7%，0.5%）。
看吧，所以我们可以用这种PCA的方法来检验bulk RNA-seq不同生物学条件下的不同生物学重复的相似性和差异性，当然DESeq2这个包本身也可以帮你来绘制PCA图，但我们还是要理解背后的原理。
2. tSNE
#install.packages('Rtsne')
library(Rtsne)
tsne.res <- Rtsne(iris.scale, check_duplicates = FALSE)
plot.data <- cbind(as.data.frame(tsne.res$Y[, 1:2]),
Species = iris$Species)
colnames(plot.data)[1:2] <- c('tSNE_1', 'tSNE_2')
ggplot(data = plot.data, aes(tSNE_1, tSNE_2)) +
geom_point(aes(color = Species)) +
theme_test()
可以看到相同物种的不同样本被很好的聚在了一起。
3. UMAP
#install.packages('umap')
library(umap)
umap.res <- umap(iris.scale)
plot.data <- cbind(as.data.frame(umap.res$layout[, 1:2]),
Species = iris$Species)
colnames(plot.data)[1:2] <- c('UMAP_1', 'UMAP_2')
ggplot(data = plot.data, aes(UMAP_1, UMAP_2)) +
geom_point(aes(color = Species)) +
theme_test()
3. 写在最后
推荐一个学习网站：清华大学鲁志实验室的Bioinformatics Tutorial，点击链接Bioinformatics Tutorial - Bioinformatics Tutorial (ncrnalab.org)
即可跳转哦！
