劉小澤寫于19.10.30 更新于2020-06-24
為何取名叫“交響樂”剖煌？因為單細胞分析就像一個大樂團拥诡，需要各個流程的協(xié)同配合
交響樂1-理解單細胞常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
這次的內(nèi)容是：實驗到后續(xù)分析做一個簡單的總結(jié)
一張非常重要的圖
這張圖和上一章的那張可以說是貫穿始終的兩張圖
實驗設(shè)計環(huán)節(jié)
在正式分析之前，關(guān)于實驗問題的探討是很有必要的乎完，最重要的一個就是技術(shù)的選擇：
- Droplet-based: 10X Genomics, inDrop, Drop-seq
- Plate-based with unique molecular identifiers (UMIs): CEL-seq, MARS-seq
- Plate-based with reads: Smart-seq2
- Other: sci-RNA-seq, Seq-Well
每種方法都有優(yōu)劣（Mereu et al. 2019; Ziegenhain et al. 2017)辱揭，目前以10X為代表的droplet-based方法由于高通量和低細胞成本成為了約定俗成的技術(shù)；Plate-based方法可以捕獲其他的一些表型信息（如細胞形態(tài)）忍疾，另外可以根據(jù)實驗?zāi)康倪M行調(diào)整，靈活性比較好谨朝；Read-based方法可以覆蓋全轉(zhuǎn)錄本卤妒，在分析可變剪切、外顯子突變等方面很有用字币；UMI-based方法可以減輕PCR擴增偏差则披。
下圖來自文章的評測：Benchmarking single-cell RNA-sequencing protocols for cell atlas projects
下一個問題就是：到底要捕獲多少細胞？測序要測多深洗出？
答案言簡意賅：As much as you can afford to spend.
如果再補充一下這個答案就是：想要發(fā)現(xiàn)罕見細胞類群士复，就要多獲得細胞；想要探索潛在的微小差異翩活，就要加大測序深度阱洪。目前常用的droplet-based儀器可以捕獲1萬到10萬細胞，測序深度是每個細胞1000到10000 UMIs隅茎，在經(jīng)濟條件一定的前提下澄峰，它們之間一般是成反比。另外它還要權(quán)衡高細胞捕獲通量和影響捕獲效率的“雙細胞比例”辟犀。
實驗設(shè)計和常規(guī)轉(zhuǎn)錄組類似
也是要考慮一個實驗條件下多個生物重復(fù)俏竞，而且實驗條件最好不要混雜批次。需要注意的是：生物重復(fù)不是指的單個細胞堂竟，而是指的提供細胞的供體（donors）或者細胞培養(yǎng)體系（cultures）
獲得表達矩陣（count matrix）
和常規(guī)轉(zhuǎn)錄組一樣魂毁，單細胞轉(zhuǎn)錄組也是需要得到表達矩陣，才能進行下游分析出嘹。表達矩陣包含的信息就是：每個細胞中比對到每個基因的UMIs或者reads數(shù)席楚。有一點需要注意：它的定量方法和具體的實驗技術(shù)相關(guān)
10X的數(shù)據(jù)：使用
CellRanger軟件，基于STAR比對到參考基因組税稼，然后統(tǒng)計每個基因的UMIs數(shù)量Pseudo-alignment方法（如
alevin）：就像之前用的salmon烦秩、kallisto意思一樣，不需要比對參考基因組郎仆，節(jié)省時間只祠、內(nèi)存-
對于一些高度multiplexed的方法：可以使用scPipe 包：提供了一套綜合的分析流程，利用
Rsubread比對扰肌，然后統(tǒng)計每個基因的UMIs數(shù)量multiplexed：翻譯叫做”多路復(fù)用“抛寝，即：large numbers of libraries to be pooled and sequenced simultaneously during a single run，可以節(jié)省成本和時間
CEL-seq、CEL-seq2數(shù)據(jù)：scruff 包可以專門分析
read-based方法：可以使用常規(guī)bulk 轉(zhuǎn)錄組定量的流程（比如smartseq2就可以用hisat2+featureCounts）
任何包含spike-in轉(zhuǎn)錄本的數(shù)據(jù)：spike-in序列都要在比對盗舰、定量之前加到參考基因組中
定量結(jié)束后晶府，一般是先導(dǎo)入表達矩陣然后創(chuàng)建一個SingleCellExperiment對象（例如：read.table() + SingleCellExperiment()）。除此以外钻趋，還有一些特定的文件格式需要用特定的包川陆，比如DropletUtils可以分析10X數(shù)據(jù)，tximport/tximeta 可以分析pseudo-alignment數(shù)據(jù)
需要注意
- 如果分析的是人類數(shù)據(jù)并且加入了ERCC爷绘，我們很多時候直接用
^ERCC在行名中進行正則匹配书劝，但是這時要小心进倍，因為ERCC基因家族在人類基因組注釋中確實存在土至，很有可能將真的基因作為外源轉(zhuǎn)錄本進行分析。這個問題可以通過將表達矩陣的行名設(shè)置為Ensembl,或Entrez來解決 - 一些定量工具會統(tǒng)計表達矩陣中的reads比對率猾昆，會存在一些未必對的情況陶因。盡管這些信息可以用作質(zhì)控，但這些數(shù)值如果被誤認為是表達量信息垂蜗，那么就會干擾下游分析楷扬。因此在進行下游分析之前，這部分信息可以去掉或者保存在
colData中
數(shù)據(jù)處理與下游分析
- 首先進行質(zhì)控：去掉低質(zhì)量細胞贴见。這些細胞可能在建庫環(huán)節(jié)被破壞烘苹，可能沒有被有效捕獲（這就是所謂的“dropout”）。一般會統(tǒng)計：每個細胞的全部count數(shù)片部、spike-in或線粒體reads比例镣衡、檢測到基因的數(shù)量
- 表達矩陣歸一化：為了減小細胞文庫的偏差（可能由于細胞捕獲效率不同、測序深度的差異而造成文庫大小差異）档悠，把細胞們放在同一起跑線上廊鸥，才能進行下面的細胞相似性比較，后面再根據(jù)相似性進行細胞分群辖所。一般是基于log轉(zhuǎn)換（當然有的函數(shù)也涉及了一些size factor的計算）惰说，從而對均值-方差進行校正
- 挑選一些特征基因（一般是高變化基因HVGs，Highly Variable Genes）進行下游分析缘回。原理是根據(jù)每個基因在細胞之間的差異構(gòu)建變化模型吆视，然后找那些變化差異大的基因。使用HVGs不用全部基因的原因一是為了減少計算量酥宴，二是減少不感興趣基因（比如在細胞之間沒什么差異）對分析產(chǎn)生的噪音
- 降維處理：讓數(shù)據(jù)更“緊湊”啦吧，一般是線性降維PCA+非線性降維tSNE/umap。PCA一般是先獲得初步的低維數(shù)據(jù)（可能會挑出幾十個主成分）幅虑，然后傳給t-SNE進一步壓縮丰滑，進行可視化
- 細胞聚類：根據(jù)細胞歸一化后的表達量相似性分成組，然后根據(jù)每個組marker基因（可理解為這一群細胞的標志性基因）的差異表達對分群進行生物學(xué)定義
比如用來自scRNAseq的一個droplet-based的視網(wǎng)膜數(shù)據(jù)【Macosko et al. (2015)】，就從原始矩陣得到了分群結(jié)果褒墨，可以看到這里不使用Seurat也能做質(zhì)控炫刷、挑高變化基因等等
library(scRNAseq)
sce <- MacoskoRetinaData()
# 質(zhì)控
library(scater)
is.mito <- grepl("^MT-", rownames(sce))
qcstats <- perCellQCMetrics(sce, subsets=list(Mito=is.mito))
filtered <- quickPerCellQC(qcstats, percent_subsets="subsets_Mito_percent")
sce <- sce[, !filtered$discard]
# 歸一化
sce <- logNormCounts(sce)
# 挑高變化基因
library(scran)
dec <- modelGeneVar(sce)
hvg <- getTopHVGs(dec, prop=0.1)
# 降維
set.seed(1234)
sce <- runPCA(sce, ncomponents=25, subset_row=hvg)
sce <- runUMAP(sce, dimred = 'PCA', external_neighbors=TRUE)
# 聚類
g <- buildSNNGraph(sce, use.dimred = 'PCA')
sce$clusters <- factor(igraph::cluster_louvain(g)$membership)
# 可視化
plotUMAP(sce, colour_by="clusters")
最后注意這里的分群并不一定是真正有生物學(xué)意義的，根據(jù)不同的參數(shù)可以得到不同的分群結(jié)果郁妈，而且這里看到的多個小群也有可能是同屬一個大群浑玛。最后的分群需要計算+生物知識共同實現(xiàn)。
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩噩咪，打造生信星球顾彰，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺胃碾。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
