生物本身就是一個(gè)復(fù)雜的體系，我們必須使用多維數(shù)據(jù)來對之進(jìn)行準(zhǔn)確描述冤议，而這一點(diǎn)在生物信息學(xué)領(lǐng)域體現(xiàn)的尤為突出师坎，以經(jīng)典的轉(zhuǎn)錄組數(shù)據(jù)為例：
- bulk RNA-seq作為一種經(jīng)典的基因相對表達(dá)量定量方法廣泛應(yīng)用于各種研究中胯陋，如果我們在一個(gè)課題中對4個(gè)樣品進(jìn)行了bulk RNA-seq測序，那么我們將得到一個(gè)4*N維的矩陣义矛，這里N就是全基因組范圍內(nèi)能轉(zhuǎn)錄的基因盟萨。
- scRNA-seq是近十年來蓬勃發(fā)展的一種分辨率更高的RNA-seq捻激，如果我們進(jìn)行測序的組織樣本中含有10000個(gè)細(xì)胞，那么我們在最后將得到一個(gè)10000*N維的矩陣垃杖。
數(shù)據(jù)紛繁復(fù)雜丈屹，我們不可能將所有的數(shù)據(jù)都進(jìn)行分析，無論是從時(shí)間和空間上都是不現(xiàn)實(shí)的苞冯，因此在生物信息學(xué)（尤其是NGS組學(xué)）中數(shù)據(jù)降維幾乎是一個(gè)無法逃避的主題舅锄。
那么究竟什么是數(shù)據(jù)降維呢司忱？簡單來說就是只保留那些在不同的組別（對bulk RNA-seq而言就是不同的樣本，對scRNA-seq而言就是不同的細(xì)胞）中具有較大差異的特征鳍烁，而舍棄掉那些非翅；模“均一”的特征梳玫，所謂“撥云見日”。
數(shù)據(jù)降維究竟有什么實(shí)際應(yīng)用呢姚垃？一個(gè)最簡單的應(yīng)用就是找到一個(gè)群體中相似的個(gè)體积糯，而這一點(diǎn)也被大家廣泛地應(yīng)用到NGS測序數(shù)據(jù)質(zhì)量評估中：一般而言谦纱，我們會給同一個(gè)生物學(xué)條件下準(zhǔn)備多個(gè)生物學(xué)重復(fù)，我們更期待不同樣品的組內(nèi)差異小于組間差異绍昂，只有這樣才能說明我們的數(shù)據(jù)可重復(fù)性好偿荷，實(shí)驗(yàn)條件更加穩(wěn)定唠椭。
1. 常見數(shù)據(jù)降維方法
生物學(xué)中的數(shù)據(jù)降維方法都是從機(jī)器學(xué)習(xí)領(lǐng)域借鑒過來的贪嫂，目前常用的主要包括線性降維的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維特征彼此正交话速，并且可以盡可能多的解釋數(shù)據(jù)中的variation芯侥。這樣得到新的k維特征就被稱為k個(gè)主成分柱查。k個(gè)主成分實(shí)際上對應(yīng)著我們從輸入數(shù)據(jù)估計(jì)出的協(xié)方差矩陣最大的k個(gè)特征值對應(yīng)的特征向量物赶。PCA在在生物信息學(xué)分析中，是樣本量不多的bulk RNA-seq轉(zhuǎn)錄組的可視化最常用的方法酵紫。
為了便于讓大家更好的簡單理解什么是PCA奖地，在這里我用兩張PPT做簡單說明：
- tSNE與UMAP
t-SNE與UMAP主要的優(yōu)勢就是保持非線性的局部結(jié)構(gòu)的能力参歹，比較適用于樣本量較大犬庇，樣本之間相互關(guān)系比較復(fù)雜的數(shù)據(jù)臭挽，所以在單細(xì)胞轉(zhuǎn)錄組的可視化中有大量應(yīng)用。需要注意的是葬荷，這兩種將為的方法通常只應(yīng)用于數(shù)據(jù)的可視化。
2. 代碼實(shí)現(xiàn)
在本部分我將從數(shù)據(jù)的預(yù)處理到降維分析可視化的過程用R語言進(jìn)行展示举反。
-
數(shù)據(jù)預(yù)處理
在進(jìn)行降維分析之前扒吁，首先必須要做的一個(gè)工作就是數(shù)據(jù)的預(yù)處理瘦陈，其中非常關(guān)鍵的一步就是scale（有時(shí)可以翻譯為“歸一化”）晨逝。簡單來說，未經(jīng)處理的數(shù)據(jù)可能具有非常大的變異度支鸡，這可能會對后續(xù)的分析產(chǎn)生影響牧挣，具體來說瀑构，以RNA-seq為例：不同的基因具有不同的表達(dá)水平寺晌，有些基因表達(dá)水平很高澡刹，有的很低；同一個(gè)基因在不同的生物學(xué)條件下表達(dá)的差異也不同陆赋，有的“巋然不動”攒岛，有的“聞風(fēng)而動”阵子，甚至具有非常大的差異思杯。用魯志老師的tutorial來說挠进，對涉及距離計(jì)算的模型而言领突，在實(shí)際應(yīng)用中君旦，不同特征的數(shù)值大小不同嘲碱，會在求距離時(shí)造成很大的影響。例如恕稠，一個(gè)基因的表達(dá)量很高鹅巍，另一個(gè)基因的表達(dá)量很低骆捧，如果不進(jìn)行數(shù)據(jù)的歸一化敛苇，表達(dá)量高的基因就會主導(dǎo)距離的計(jì)算接谨。
在R語言中我們常用的就是使用z-score里進(jìn)行歸一化：
這里我們使用鳶尾花數(shù)據(jù)集來進(jìn)行演示脓豪，該數(shù)據(jù)集為R包datasets的一個(gè)內(nèi)置數(shù)據(jù)集忌卤，記錄了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()函數(shù)進(jìn)行歸一化：
iris.scale <- scale(as.matrix(iris[, 1:4]),
scale = TRUE,
center = TRUE)
-
數(shù)據(jù)降維
1. PCA
首先是PCA，這里介紹兩種方法驰徊，分別基于stats包的prcomp()函數(shù)與FactoMineR包的PCA()函數(shù)笤闯。
#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()
我們可以看到同一個(gè)物種的不同個(gè)體基本都被聚在了一起，而不同物種之間“涇渭分明”棍厂。（versicolor和virginica未被很好分開是因?yàn)檫@兩個(gè)物種本身就具有相似的特征）
#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()
結(jié)果類似颗味，不過看起來好像PC2反過來了∥可以看到這里有兩個(gè)數(shù)字：73%和22.9%浦马，你可以簡單理解為每個(gè)主成分能夠解釋的數(shù)據(jù)變異度比例时呀，根據(jù)前面的PPT磺陡，理論上來講PC1解釋度最高坞靶，PC2次之……所以我們嘗試看一下PC3和PC4對鳶尾花個(gè)體的區(qū)分度如何：
fviz_pca_ind(pca.res,
axes = c(3, 4),
geom = 'point',
habillage = iris$Species,
addEllipses = TRUE) +
theme_test()
可以說基本沒有分開辫封，原因就在于這兩個(gè)主成分沒有很好的解釋整個(gè)數(shù)據(jù)集的變異度（3.7%，0.5%）欣福。
看吧，所以我們可以用這種PCA的方法來檢驗(yàn)bulk RNA-seq不同生物學(xué)條件下的不同生物學(xué)重復(fù)的相似性和差異性栖博，當(dāng)然DESeq2這個(gè)包本身也可以幫你來繪制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. 寫在最后
推薦一個(gè)學(xué)習(xí)網(wǎng)站：清華大學(xué)魯志實(shí)驗(yàn)室的Bioinformatics Tutorial，點(diǎn)擊鏈接Bioinformatics Tutorial - Bioinformatics Tutorial (ncrnalab.org)
即可跳轉(zhuǎn)哦粗俱！
