对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择
原文鏈接：
- https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6072887
摘要
單細胞RNA測序技術的發展加深了我們對于細胞作為功能單元的理解，不僅能基于成百到成千上萬個單細胞的基因表達譜得到新的結論，還能發現新的具有特異基因表達譜的細胞群（這在傳統轉錄組測序數據中是很難發現的）。
但是合適的分析和使用單細胞RNA測序后的大量數據并不容易，需要了解從準備細胞到獲得合理的結果之間采用的實驗和分析方式。
在本篇綜述中，作者討論了這些新技術的基本原理，重點強調單細胞轉錄組分析中的重要概念。具體地說就是總結了單細胞過濾的質控方法、由于mRNA捕獲率低需要進行的標準化和歸一化方法和用于降維數據繪制二維圖的聚類和可視化的算法。
-
Hemberg-lab單細胞轉錄組數據分析（九）- Scater包單細胞過濾
-
Hemberg-lab單細胞轉錄組數據分析（十）- Scater基因評估和過濾
分析流程
圖示：對得到的DGE數據進行質控→依據轉錄本和細胞數、轉錄本比對到參考基因組的比例和線粒體mRNA的比例對數據進行質控→用標準化和歸一化來除去批次效應和誤差→降維和可視化。
正文
早期，單細胞基因表達圖譜用于研究少數單個細胞里面特定選擇的一些轉錄本。后來，高通量測序技術和高效細胞分離方法的發展促進了現代單細胞測序平臺（例如Fluidigm C1，DropSeq，Chromium 10X，SCI-Seq以及過去十年中開發的許多其他單細胞測序技術）的發展。(Hemberg-lab單細胞轉錄組數據分析（二）- 實驗平臺)
表格展示了現在流行的單細胞測序方法，這些方法都依賴于細胞解離技術。價格很大程度上取決于測序細胞的數量、測序深度以及測序平臺。表中的價格是單細胞文庫制備價格區間的最低價。
| Fluidigm C12 | ～1,000 cells | cells capture insize-specific chambers | must know the size of cells of interest; allows forstaining and imaging prior to cell rupture | 96- or 800-chamber units areavailable | an average of 6,606 genes/cell (no data on percentage) | $1.70 |
| DropSeq | ～150,000 cells/run | droplet-based separation | remains the mostcost-effective and most customizable | ～5% of cells per run (approximately 7,000 cells) | ～10.7% of the cell’s transcripts | $0.06 |
| Chromium 10X | ～1,700 cells/run | droplet-based separation | the most commerciallysuccessful method; almost fully automated | ～65% of cells per run (approximately 1,000 cells) | ～14% of the cell’s transcripts | $0.10 |
| SCI-Seq | ～500,000 cells (depends on experimental design) | FACS sorter; cellsare never singly isolated | combinatorial indexing of individual methanol-fixedpermeable cells | 5%–10% of cells | ～10%–15% of the cell’s transcripts | $0.05–$0.143 |
表1
這些技術能夠一次檢測到數百至成千上萬個單細胞的轉錄組數據。在逆轉錄或后續步驟中，這些方法都是使用DNA barcode標記mRNA，來判斷轉錄本來自于哪個細胞。盡管每種技術在分離細胞和標記mRNA方面處理方式都不一樣，但是它們都使用相似的計算流程處理單細胞轉錄組數據。
在本篇綜述中，作者討論了流程中最常使用的一些算法，并以DropSeq為主要示例，因為它是性價比最高并且使用最廣泛的單細胞基因表達平臺（表1）。不過這些方法也同樣適用于絕大多數使用DNA barcode 來標記mRNA來源的單細胞測序技術。
單細胞轉錄組主要應用于：
-
深入了解組織中的細胞異質性。
-
鑒定未知的細胞類型。
-
鑒定已知細胞類型的亞型，原理是在感興趣的細胞群中尋找差異基因表達模式。
-
從稀有細胞群中分離出信號，這些信號在普通轉錄組中很難被分離出來。
-
給未知Maker的細胞類型推斷可能的Maker，如細胞表面蛋白等。
這個原理是由于單細胞轉錄組分析中會根據細胞間差異表達的基因進行聚類，這樣就可以把對聚類影響最大的基因認為是感興趣細胞群可能的Maker。
-
細胞譜系和分化調控的研究，例如可以誘導一組干細胞分化，并在不同時間點進行單細胞測序可得到分化各個階段的“snapshot”。這些snapshot可用于推斷細胞到達終末分化狀態所遵循的軌跡和在每個分支點受到差異調控的關鍵基因。
不過這些應用幾乎都依賴于由頂尖生物信息學實驗室開發的并釋放的一些特定的算法。而在本篇綜述中，作者重點介紹了用這些特定算法之前必須做的——對數據進行質量控制和標準化，并討論了簡單的細胞聚類和可視化的算法。
基于Droplet方法生成單細胞基因表達數據集
DropSeq和商用的10X都是基于液滴（Droplet）的方法來產生單細胞基因表達數據集。基于液滴型方法是使用微流體芯片將單個細胞與單個凝珠(beads)包裹進油囊化液滴中，理想情況下每個液滴最多包含一個細胞。(Hemberg-lab單細胞轉錄組數據分析（七）- 導入10X和SmartSeq2數據Tabula Muris)
凝珠 (beads)上附有數目極多的DNA寡核苷酸探針（DNA oligos）。DNA寡核苷酸的3’末端是一個poly（T）尾巴，可用于捕獲細胞中的mRNA （生信寶典注：更準確的說是捕獲細胞中有poly-A尾巴的RNA，既有mRNA，也有ncRNA）；5’末端是用于標記細胞的cell barcode序列，一個凝珠上結合的所有寡核苷酸的cell barcode都一樣；中間還有一個具有高度多樣性的唯一分子標識符（unique molecular identifier，UMI），磁珠上的每個寡核苷酸的UMI都不一樣 (生信寶典注：這個不能同意，UMI的種類是少于一個凝珠上所有寡核苷酸的數目的)。（見下圖1）
在液滴中，細胞破裂，磁珠上的DNA寡核苷酸捕獲并標記釋放的轉錄本；隨后液滴破裂，所有細胞一起進行、逆轉錄、PCR擴增并通過高通量平臺測序。測序得到序列與參考基因組進行比對，對應到注釋的基因；再根據比對序列上的cell barcode區分來自同一個細胞序列。最后使用UMI計算每個細胞中表達的單個基因的轉錄本的拷貝數，從而可以生成基因表達矩陣（DGE, digital gene expression），這個矩陣就是包含了細胞barcode和基因counts的表格文件。
圖一：Drop-Seq Bead的結構和所得序列文庫
（A）DropSeq單細胞測序珠的結構。
DNA寡核苷酸的3’末端是一個poly（T）尾巴，可用于捕獲細胞中的mRNA ；5’末端是用于標記細胞的cell barcode序列，一個凝珠上結合的所有寡核苷酸的cell barcode都一樣；中間還有一個具有高度多樣性的唯一分子標識符（unique molecular identifier，UMI），磁珠上的每個寡核苷酸的UMI都不一樣。
（B）測序文庫的結構。
-
紅色：PCR引物，也可用作測序引物；
-
綠色和藍色：來自珠子的細胞和分子barcode；
-
橙色：捕獲的帶有poly（A/T）尾巴的轉錄本。
從如此復雜的測序數據中得到的可靠結論取決于后續的計算分析。
大多數可用的單細胞測序算法沒有圖形用戶界面。因此分析單細胞數據需要一些編程的基礎，這樣才能更好地去執行比對、聚類和可視化數據的算法。
除此之外，對目標細胞有深刻的生物學認識是很有必要的，這樣才能正確地解釋數據并選擇合適的質控參數等。
在單細胞測序方面的生信專家需要在分析過程中對應用的算法能選擇適當的閾值，并避免產生誤導性的結果，才可以對其分析結果進行有意義的生物學推斷。
質控指標
基于液滴的實驗可以視為對單個液滴內的單個細胞進行的成千上萬次的獨立實驗，所以必須要對數據進行質量控制（QC）去除低質量數據。而QC是通過使用不同的指標來判斷并過濾掉不合格（如技術問題或細胞質量問題等導致的）的數據。
QC指標—每個細胞檢測到的轉錄本數量或測序序列比對到參考基因組的比例
QC參數的閾值在不同分析中不一定相同，閾值的設置取決于測序的細胞或組織。
常見的QC指標是每個細胞的轉錄本數量或每個細胞能比對到參考基因組的測序序列的百分比。
若細胞的轉錄本數量低于或高于定義好的閾值，該細胞會被標記為異常細胞并從分析除去；閾值既可以由分析者自定義（例如，細胞的轉錄本少于20個或者超過5,000），也可以由程序自動判斷（例如，轉錄本總數大于所有細胞平均轉錄本數目2倍標準差的細胞需要被移除，cells with a sum of transcripts larger than 2 SDs from the mean are removed）。
因為如果一個細胞包含大量的轉錄本，可能是由于doublets（即兩個或兩個以上的細胞懸浮在一個液滴中）造成，這種數據要從分析中除去；如果一個細胞的檢測到的轉錄本數量很少，意味著捕獲質量較差，這可能是因為細胞死亡、細胞過早破裂或者是捕獲了從細胞中逸出并漂浮在細胞懸液中的隨機mRNA。（單細胞預測Doublets軟件包匯總-過渡態細胞是真的嗎？）
也可以應用其他QC指標，例如，直接刪除表達某個特定基因的所有細胞，這樣可以刪除不感興趣的污染細胞，或者更復雜一些，只包含兩個或更多特定基因表達達到一定比例的細胞。
在確定QC閾值時，必須考慮所分析的組織的多樣性。例如，在設計實驗研究血液中轉移的癌細胞時，癌細胞的數量相較于正常血細胞的數量而言非常低，因此必須調整QC指標中的轉錄本數量（counts of transcripts）。在該組織中血細胞是優勢細胞，但與活躍的癌細胞相比，它們的表達卻被認為處于相對靜止狀態，具有相對較低的RNA量。故而如果設置閾值為刪除那些轉錄本數量高于平均值2倍標準差的細胞，癌細胞因為轉錄活性比較高，就可能會被誤認為是doublets，并被全部移除。（生信寶典注：相比于很多人生搬硬套Seurat示例數據中的200，2500的篩選標準，采用n倍標準差是適應性更廣的方式，尤其是不關注稀有類型時。如果自己比較了解，還是需要好好看下數據分布再定標準。如果不了解，可以先松后緊，根據最后結果再回來看轉錄本數目異常的細胞聚類在什么地方再做評判。）
QC指標—線粒體基因的數量
另一個常見的QC參數是線粒體基因的數量。高比例的線粒體基因表達細胞處于應激狀態的指標之一，因此分析中通常需要移除線粒體基因表達占比較高的細胞，因為大多數實驗不研究這一類特殊狀態的細胞。
但是，與轉錄本數量一樣，此參數高度依賴于組織類型和所研究的問題。例如，由于心肌細胞的高能量需求，心臟中總mRNA的30％是線粒體，而低能量需求的組織中占比則為5％或更少。故而線粒體mRNA占30％在心肌細胞表示健康，但在淋巴細胞表示不正常。(生信寶典注：這一步篩選也不要受Seurat文檔影響太深，參數都是可以改的，只要有合適的原因。最近一期的單細胞培訓，這個也是討論的重點，國內外學者濟濟一堂討論這個參數選擇。)
QC指標—篩選基因
根據實驗的目的，也可以添加基因特異性的QC指標。在所有細胞內表達量都很低并在細胞類型之間無統計意義的基因，可以考慮設置閾值過濾掉，減少后期的計算量：設置每個細胞內的基因count閾值（例如，基因在在每個測序的細胞中，count值都小于5）或設置所有或一個細胞子集中該基因count總和的閾值（例如，所有測序細胞中的基因∑count<=300）。
雖然排除此類基因將加快計算過程，不過可能會丟失一些表達差異很小但對數據差異有貢獻的基因。（生信寶典注：不排除有一些基因表達量比較低，并且較小的變化幅度就可以帶來有意義的生物效果。但表達低的基因本身檢測的噪音也大，比較難區分哪些是生物差異，哪些是技術差異。私以為，原文這句描述有誤。）
數據標準化和歸一化
在分析測序數據時，如果要對多批測序數據進行相互比較，需要消除批次效應。這些批次效應可能是由不可避免的技術差異引起的，例如將樣本冷凍存放時間、反復凍融的次數、提取RNA的方法、測序深度等。
研究人員應努力保持這些實驗和測序過程中的變量恒定。但是基于液滴的測序還包含數千個單獨的細胞實驗，因此在標準化時還必須考慮細胞特異性偏差，以便能夠將一個細胞與另一個細胞進行比較。
特異性偏差是由mRNA捕獲效率引起的，在所有液滴中mRNA分子沒有以相同比例被磁珠捕獲，這被稱為“dropout events”，它也是數據稀疏的主要原因，數據稀疏將在下一段深入討論。
此外在bulk RNA測序中，需要被標準化的多批數據幾乎來自相似的生物材料（例如將血細胞與血細胞進行比較），但是在單細胞測序中，單個細胞并不屬于同一類型，這就需要調整標準化的參數以保留細胞間差異，同時還要消除技術差異帶來的批次效應和細胞特異性偏差。（DESeq2差異基因分析和批次效應移除）
mRNA捕獲效率很低（例如，DropSeq被認為最多能捕獲每個細胞**10％**左右的mRNA），這是液滴型單細胞測序數據的分析面臨的最大挑戰。由于這些“dropout events”，DGE矩陣大部分數據都會是0，這就是數據稀疏了。因此在解釋數據之前，標準化和歸一化至關重要。不過，這需要假設細胞在生物學上不需要嚴格準確（Unfortunately, this requires making assumptions about the cells that can be biologically inaccurate）。
一種可接受的標準化測序數據的方法是利用管家基因進行比較。
首先基于文獻資料和對測序的生物樣品的了解，選擇一個管家基因用于后續標準化。假定所選的管家基因在所有細胞中均以相同的水平表達，然后對測序數據進行歸一化使所選的管家基因的表達水平在所有細胞中均相等。（什么？你做的差異基因方法不合適？）
但是這個方法也可能不準確，因為這些持家基因在不同細胞中表達量并不總是一致的。另一個思路是基于在所有或一部分細胞中所有表達無差異的基因進行標準化。這一方法基于所有細胞或部分細胞之間表達無差異的所有基因均在所有或部分細胞中均等表達的假設進行歸一化，并推斷出每個細胞的歸一化因子來標準化轉錄本的計數。
降維和可視化
PCA
對基因表達譜標準化后，應用無偏聚類的算法可以確定哪些細胞更為相似。
Principal component analysis（PCA）通常是首選的聚類算法，因為它是一種相對簡單的線性降維算法，可以預測多維數據的相關性，具體的在單細胞分析中指只需要依賴高可變基因的表達譜就可以預測細胞間的相似關系。
PCA把相關的基因合并到“metagene”或主成分（PC）中。PC1解釋最大的數據差異，具有最大的標準差（例如對于一個實驗，細胞之間30％的差異由定義了PC1的基因解釋），PC2解釋了數據的第二大部分差異（例如，細胞之間20％的差異可歸因于PC2中的基因，而8％則歸因于PC3中的基因），然后依此類推，簡單來說PC的排名就是解釋數據差異貢獻的順序，其中PC1是排名最高的PC，同時也說明PC排名越低，對解釋數據差異的貢獻就越小。
關于PCA的解釋，還是推薦我們的文章：
-
PCA主成分分析實戰和可視化 | 附R代碼和測試數據
-
用了這么多年的PCA可視化竟然是錯的！！！
使用排名較低的PC一般都沒什么好處，因為它既增加了計算量，又幾乎沒有將任何信息添加到細胞間差異的展示中。因此，決定用于可視化的PC數非常重要。常用的判斷方式就是繪制knee圖或elbow圖，如下圖所示。
圖中展示了每個主成分的標準差，代表每個PC對數據差異解釋的貢獻度。PC4、PC5和PC6都在拐點附近，說明推薦使用前四個、前五個和前六個主成分用于后續分析。
-
Hemberg-lab單細胞轉錄組數據分析（十一）- Scater單細胞表達譜PCA可視化
-
Hemberg-lab單細胞轉錄組數據分析（十二）- Scater單細胞表達譜tSNE可視化
t-SNE
t-Distributed stochastic neighbor embedding（t-SNE）是一種常見的可視化方法。
它使用機器學習的算法來降低維數，非常適合將高維數據放到二維或三維空間中可視化展示，并且不會丟失細胞之間的相對距離的信息。
例如，如果發現用七個PC可以很好地表示細胞的多樣性，就得需要七個軸或維度來展示細胞的空間分布。t-SNE能維持細胞在七維空間的關系并在二維圖上展示細胞，所以在七維圖上相鄰的細胞在二維圖上仍然相鄰。同時PCA分析是線性的，t-SNE是非線性降維方法。（還在用PCA降維？快學學大牛最愛的t-SNE算法吧（附Python/R代碼））
注意事項：有關數據生成效率和可替代的單細胞平臺
在本篇綜述中討論的計算方法主要是用于基于液滴的分離方法，例如DropSeq和Chromium 10X。
不過大多數單細胞測序平臺都是用特異的DNA barcode標記每個細胞的mRNA，從而得到每個細胞的基因表達信息的，所以上述介紹的類似的原理和算法也可用于其他方法的數據集。但是要注意不同平臺之間始終存在著技術的區別或儀器的差異。
Valihrach等人發表了一篇綜述(https://www.ncbi.nlm.nih.gov/pubmed/29534489)，詳細地描述了細胞分離、標記和DNA擴增的平臺和方法，并且討論了每種平臺的基本原理以及不同方法的優點和缺陷。
例如，在SCI-Seq中，先使用酒精固定細胞，使其具有滲透性，再使用流式細胞儀對固定的細胞進行分選，最后將特定數量的細胞分配到多孔板的每個孔中（見表1）。
每個孔中細胞的mRNA均通過反轉錄結合了該孔特有的寡核苷酸，然后合并所有孔中的細胞，并以較低的密度對細胞進行另一輪熒光激活細胞分選（FACS），然后添加第二個獨特的孔特異性barcode，為每個細胞創建唯一的barcode組合。此過程可以再次重復，幫助降低同一barcode的組合標記兩個細胞的可能性。
這種barcode組合標記單個細胞的方法需要專門的算法來得到DGE矩陣，因為與基于液滴的方法相比，單個細胞不是由單個barcode確定，而是由barcode的特定組合確定的。值得注意的是，由于這種方法至少需要兩輪細胞分選，可能會對細胞產生更大的影響并且影響基因的表達。
另一個示例是設置轉錄本/細胞參數的數目以排除doublets，因為每種方法都會有不同的doublets比例。
在Fluidigm C1系統中，單個細胞被隔離在特定大小的區室捕獲，在隔離的中等大小的96-區室中對細胞進行顯微鏡檢查后，doublets的比例從7％下降至3％。這個比例沒有為0是因為細胞有時會在隔離室中相互堆疊，使它們看起來像單個細胞，顯微鏡檢查就無法發現這種堆疊的細胞。
如果經顯微鏡檢查后還有3％以上的細胞或未經檢查的數據中有7％的以上的細胞的轉錄本的數量顯著更高（例如比轉錄本的平均值高2倍的標準差以上），這表明可能這批細胞是由少數的有轉錄活性細胞和大多數的無轉錄活性的細胞組成，或者可能是由于doublets的比例很高，如果是在這種情況，可能需要換更小的隔離室來選擇細胞了。
當前主流的單細胞分析流程是以細胞之間表達差異最大的基因為基礎。這有助于發現未知細胞類型的基因marker。但是如果研究人員打算研究非常相似的細胞類型，或想從一種主要細胞類型中找到亞型，則可以在分析之前對這些細胞進行分選和增加感興趣的細胞數量，從而提高分析精度。（單細胞分群后，怎么找到Marker基因定義每一類群？）
即使熒光激活細胞分選（FACS）被證明對基因表達的影響很小，分選仍延長了細胞不在最佳培養條件下而在單細胞懸浮液中的時間，這可能會對細胞造成影響并可能改變mRNA和線粒體mRNA的表達。此外，使細胞通過小區室、微流控分選儀或細胞分選儀會導致應激反應和影響對應激更敏感和更容易死亡的一些細胞類型。因此，在基于液滴的單細胞測序實驗中，自身比較脆弱的細胞亞型等可能很難被發現，特別是如果在單細胞分離之前對細胞進行了分選。
結論
在本篇綜述中討論了一些重要概念，這些概念在單細胞基因表達數據分析和根據細胞類型或條件來選擇參數非常重要。另外還提供了一些其他類型技術的示例，使得分析可以用于基于非液滴的單細胞測序數據。
分析流程首先從原始測序文件生成包含每個細胞的基因計數DGE矩陣開始；接下來通過QC除去可能由于doublets和細胞應激等產生的錯誤細胞；再進行標準化和歸一化解決不可比的問題（由mRNA捕獲率低導致）；然后基于細胞間高可變的基因進行降維和聚類；最后在二維或三維空間上展示數據中每個細胞與其他細胞的相關性。
通用算法，通常這些算法包含在易于使用的程序包中：
-
Seurat，這是一個基于R的程序包，可創建與其他下游算法兼容的R對象(https://satijalab.org/seurat/)；
-
scran，還包括用于細胞周期分配的算法(http://bioconductor.org/packages/release/bioc/html/scran.html)；
-
ascend，其中包括完善的新算法，提供了靈活的分析框架(https://github.com/IMB-Computational-Genomics-Lab/ascend)。
評估以上和其他一些軟件包檢測到的高變基因的準確性和精確性：https://www.ncbi.nlm.nih.gov/pubmed/29481632
后續分析算法，基于實驗目的，可選的更有針對性地下一步分析的算法：
-
Monocle*，*該算法旨在分析單個細胞的分化軌跡(http://cole-trapnell-lab.github.io/monocle-release/)；
-
SingleSplice，用于研究單細胞群體中的可變剪接(https://github.com/jw156605/SingleSplice)；
-
OncoNEM，這是一種推斷腫瘤細胞之間基于層次進化關系的工具(https://bitbucket.org/edith_ross/onconem/src)。
單細胞測序工具及其應用的大集合可在下列網站上找到，這些網站也會更新有最新推出的工具：
-
https://github.com/seandavi/awesome-single-cell
-
https://www.scrna-tools.org/
推薦閱讀
-
Seurat亮點之細胞周期評分和回歸
-
Celaref | 單細胞測序細胞類型注釋工具
-
cellassign：用于腫瘤微環境分析的單細胞注釋工具（9月Nature）
-
Nature重磅綜述 |關于RNA-seq，你想知道的都在這
-
NBT|45種單細胞軌跡推斷方法比較，110個實際數據集和229個合成數據集
-
如何火眼金睛鑒定那些單細胞轉錄組中的混雜因素
-
讓你的單細胞數據動起來！|iCellR(一)
-
Cell子刊 | 植物單細胞轉錄組綜述·植物功能基因組學的高分辨率研究方法
總結
以上是生活随笔為你收集整理的对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 这个保守的 RNA motif是病毒侵染
- 下一篇： 液滴型单细胞测序技术比较（二）