單細胞轉(zhuǎn)錄組測序--閑言碎語
時代發(fā)展的步伐總是毫不留情的將你甩在身后访递，連車尾燈都看不見汹买。當你還在沉迷于普通轉(zhuǎn)錄組數(shù)據(jù)挖掘時存崖，已經(jīng)有人悄悄的搞上單細胞了汽煮。單細胞轉(zhuǎn)錄組測序搏熄，顧名思義就是在單個細胞的分辨率基礎(chǔ)上去研究細胞內(nèi)的基因表達等茅诱，其主要目的是為了研究不同細胞類型的基因表達異質(zhì)性，從而解決相關(guān)生物學問題搬卒。談到單細胞就不得不提一下當下火爆的10x Genomics服務(wù)商了，具體參見10x Genomics翎卓。本篇文章暫時不介紹10x契邀，主要介紹單細胞轉(zhuǎn)錄組數(shù)據(jù)分析軟件Seurat。
Seurat軟件是一個R包失暴，可以說是單細胞轉(zhuǎn)錄組測序分析的明星軟件坯门，很多單細胞測序文章都會引用該軟件，引用次數(shù)也是杠杠的逗扒，而且也有詳細的在線教程古戴。本文也主要是根據(jù)其教程介紹一下使用Seurat軟件分析一個樣本的單細胞轉(zhuǎn)錄組數(shù)據(jù)的步驟及注意事項，供大家討論矩肩。
Seurat軟件介紹
導(dǎo)入分析需要的包
library(Seurat)
library(dplyr)
library(magrittr)
Seurat軟件提供了很友好的函數(shù)可以直接讀取10x Genomics的輸出結(jié)果
matrix <- Read10X_h5("filtered_gene_bc_matrices_h5.h5")#該函數(shù)直接讀取10x的h5文件现恼，
#并返回一個稀疏矩陣
#也可以使用下面這個函數(shù)導(dǎo)入數(shù)據(jù)
matrix <- Read10X("data/matrix_dir")#data/matrix_dir目錄下包含文件
#matrix.mtx, genes.tsv, and barcodes.tsv
導(dǎo)入文件后便可以創(chuàng)建Seurat對象
filter_10x_object <-CreateSeuratObject(raw.data = matrix,min.cells = 3,
min.genes = 200, project = "d615")#該函數(shù)還具有數(shù)據(jù)過濾的功能，
#基因至少在min.cells個細胞中表達黍檩，每個細胞中至少表達min.genes個基因叉袍，
#可以根據(jù)具體情況決定這些閾值
創(chuàng)建完Seurat對象后，Seurat將數(shù)據(jù)保存在不同的slot中刽酱，如filter_10x_object@raw.data, filter_10x_objectt@data, filter_10x_object@meta.data, filter_10x_object@ident喳逛，其中raw.data存放的是每個細胞中每個gene的原始UMI數(shù)據(jù)，data存放的是gene的表達量棵里，meta.data存放的是每個細胞的統(tǒng)計數(shù)據(jù)如UMI數(shù)目润文，gene數(shù)目等，ident此時存放的是project信息殿怜。
由于技術(shù)原因典蝌，一個GEM中可能會包含2個或多個細胞，也可能不包含細胞稳捆，這時候可以通過觀察每個barcode中的基因數(shù)目或UMI數(shù)目來判斷赠法。
GenePlot(object = filter_10x_object, gene1 = "nUMI", gene2 = "nGene")
上圖展示的是每個barcode中的基因數(shù)目和UMI數(shù)目的關(guān)系，一般二者都成正相關(guān)關(guān)系乔夯，有個別barcode的基因數(shù)目和UMI數(shù)目過高砖织，有可能就是包含2個細胞的GEM，可以考慮在后續(xù)分析中將其過濾掉末荐。
我們不僅僅可以觀察每個barcode的基因數(shù)目侧纯，還可以計算每個barcode中的線粒體基因含量等，從而更加仔細的觀察數(shù)據(jù)的質(zhì)量甲脏。
mito.genes <- grep(pattern = "^MT-", x = rownames(x = filter_10x_object@data),
value = TRUE)#統(tǒng)計線粒體基因
percent.mito <- Matrix::colSums(filter_10x_object@raw.data[mito.genes, ])/
Matrix::colSums(filter_10x_object@raw.data)#計算線粒體基因含量
filter_10x_object <- AddMetaData(object = filter_10x_object, metadata = percent.mito,
col.name = "percent.mito")#將線粒體基因含量添加到meta.data中
VlnPlot(object = filter_10x_object, features.plot = c("nGene", "nUMI",
"percent.mito"), nCol = 3)#繪制小提琴圖
這張圖片展示了每個barcode中基因數(shù)目眶熬、UMI數(shù)目以及線粒體基因含量的分布情況妹笆，根據(jù)上述2張圖片就可以大致確定是否需要過濾哪些數(shù)據(jù)進行后續(xù)分析。
Seurat提供了一個很好用的數(shù)據(jù)過濾函數(shù)：
filter_10x_object_filter <- FilterCells(object = filter_10x_object, subset.names =
c("nGene", "percent.mito"), low.thresholds = c(200, -Inf),
high.thresholds = c(2500, 0.05))#過濾掉小于low.thresholds或大于high.thresholds的數(shù)據(jù)娜氏，
#這個參數(shù)值對應(yīng)subset.names的參數(shù)
以上就是數(shù)據(jù)的預(yù)處理過程了拳缠，接下來就進入正式的分析階段，包括數(shù)據(jù)的標準化贸弥、歸一化窟坐、數(shù)據(jù)降維以及聚類分析等。
#數(shù)據(jù)標準化绵疲，即表達量的計算
filter_10x_object_norm <- NormalizeData(object = filter_10x_object_filter,
normalization.method = "LogNormalize", scale.factor = 10000)#在每個細胞內(nèi)哲鸳，
#對每個基因進行標準化：該基因的UMI除以該細胞內(nèi)轉(zhuǎn)錄物的總UMI并乘以10000，
#然后在使用log1p對這些值進行自然對數(shù)轉(zhuǎn)換（標準化的數(shù)值存放在data slot中）
#計算高度變化的基因
#（并不是所有基因都有有效的信息盔憨，尋找高度變化的基因既可以節(jié)省計算資源徙菠，又能概括樣本信息）
filter_10x_object_norm_var <- FindVariableGenes(object = filter_10x_object_norm,
mean.function = ExpMean, dispersion.function = LogVMR,
x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)#具體 cutoff值可以根據(jù)下圖調(diào)整
FindVariableGenes算法：首先計算基因的平均表達量，然后計算基因的離散度郁岩；接下來根據(jù)平均表達值將基因分成20塊并計算每塊的離散度的Z值婿奔。
如上圖：橫坐標代表基因的平均表達量，縱坐標代表基因的離散度的Z值问慎，標有基因名的點就是由函數(shù)中的cutoff值決定的脸秽，改變cutoff值，這些標記也會隨之改變蝴乔。
數(shù)據(jù)的線性回歸记餐、中心化和比例化：對數(shù)據(jù)進行線性回歸分析，去除不想要的變異源薇正。
中心化：首先計算基因A在所有細胞中的平均表達量片酝，然后分別將每個細胞中基因A的表達值減去平均值。
比例化：在中心化的基礎(chǔ)上挖腰，首先計算基因A在所有細胞中的中心化值后的標準差雕沿，然后分別將每個細胞中基因A的中心化值除以標準差。這些步驟都在一個函數(shù)中完成猴仑。
filter_10x_object_norm_var_scal <- ScaleData(object = filter_10x_object_norm_var,
vars.to.regress = c("nUMI", "percent.mito"))
#vars.to.regress中的參數(shù)就是需要去除的變異來源
#該函數(shù)會先進行先行回歸分析审轮，然后再進行數(shù)據(jù)的中心化和比例化，并將值存在scale.data中
單細胞轉(zhuǎn)錄組測序產(chǎn)生的數(shù)據(jù)是數(shù)萬個基因在數(shù)萬個細胞中的表達情況辽俗，屬于典型的高維數(shù)據(jù)疾渣。如果把1個基因視為1個坐標軸的話，那么一個細胞的空間位置就是在數(shù)萬個坐標軸中的定位崖飘，這樣的話相同細胞類型的細胞就應(yīng)該挨在一起榴捡，我們就可以根據(jù)細胞的空間位置判斷細胞亞群了≈煸。可是我最多也就認識三維坐標啊吊圾，咋辦达椰，能不能把這些高維數(shù)據(jù)投影到二維坐標呢，那就交給PCA和t-SNE吧项乒。PCA和tSNE都是數(shù)據(jù)降維分析方法啰劲，PCA屬于線性降維，tSNE屬于非線性降維檀何。我們先執(zhí)行PCA分析呈枉，使高維數(shù)據(jù)的信息最大程度保留在低維數(shù)據(jù)中，PCA分析利用的是保存在scale.data的值埃碱。
filter_10x_object_norm_var_scal_pca<- RunPCA(object= filter_10x_object_norm_var_scal,
pc.genes = filter_10x_object_norm_var_scal@var.genes, pcs.compute=20)#利用高可變基因
#計算20個主成分，并將數(shù)據(jù)保存在@dr$pca中
執(zhí)行完P(guān)CA分析后酥泞，就要根據(jù)PCA得分來進行聚類分析了砚殿，但是在進行聚類分析之前，需要選擇使用對少個主成分進行計算芝囤。每個主成分實際上代表的是相關(guān)基因集的信息似炎，因此確定多少個主成分是一個重要的步驟，我們可以根據(jù)PCElbowPlot函數(shù)來判斷悯姊。
PCElbowPlot(object = filter_10x_object_norm_var_scal_pca)
從上圖可以看到羡藐，拐點出現(xiàn)在10-15之間，我們可以選擇15來進行聚類分析悯许。Seurat采用的是基于圖形的聚類方法仆嗦，即利用PCA空間中的歐幾里德距離構(gòu)造一個KNN圖（數(shù)學好的可以留下來幫忙講講）。
filter_10x_object_norm_var_scal_pca_cluster <- FindClusters(object =
filter_10x_object_norm_var_scal_pca, reduction.type = "pca",
dims.use = 1:15, resolution = 1.0, print.output = 0, save.SNN = TRUE)
#resolution會影響聚類的個數(shù)
好了先壕，到此我們就知道了我們的數(shù)據(jù)中有多少種細胞亞群了瘩扼，怎么可以少得了圖片展示呢。超棒的可視化方法tSNE要上場了垃僚。tSNE的目標是將在高維空間中具有相似局部鄰域的細胞集绰，在低維空間中放在一起。
filter_10x_object_norm_var_scal_pca_cluster_tsne <- RunTSNE(object =
filter_10x_object_norm_var_scal_pca_cluster, dims.use = 1:15, do.fast = TRUE)
TSNEPlot(object = filter_10x_object_norm_var_scal_pca_cluster_tsne,do.label = TRUE,
pt.size = 0.2)#繪圖
既然我們知道了有多少種細胞亞群谆棺，那么是不是就要分析一下這些亞群間的差異性呢栽燕，交給FindAllMarkers吧。FindAllMarkers能夠同時計算所有亞群的差異性（分別計算每個亞群與剩下的所有細胞的差異性）改淑。
pbmc.markers <- FindAllMarkers(object =
filter_10x_object_norm_var_scal_pca_cluster_tsne,
only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)#根據(jù)data slot的值進行計算
得到差異表達基因后碍岔，當然要進行展示了。
好了朵夏，剩下的就是進行生物學知識挖掘了付秕，例如根據(jù)這些差異基因推斷細胞類型啊之類的。
關(guān)于單個樣本的單細胞轉(zhuǎn)錄組數(shù)據(jù)分析就介紹到這兒了侍郭，那多個樣本的分析會有什么不同呢询吴，我們下次再說吧掠河。
單細胞轉(zhuǎn)錄組數(shù)據(jù)分析的小伙伴一起來討論吧
