一区二区三区日韩精品-日韩经典一区二区三区-五月激情综合丁香婷婷-欧美精品中文字幕专区

分享

小澤看文獻 | 誠意滿滿的綜述之單細胞轉錄組分析最佳思路

 漠藩 2020-08-06

劉小澤寫于2020.5.6-5.8

歷時三天,終于理解完!不僅僅是綜述,更不是純翻譯文,而是將其中重要的知識點和之前個人所學結合起來
”一陽指和獅吼功合并為一整招“

1 文章信息

題目:Current best practices in single-cell RNA-seq analysis: a tutorial

發(fā)表日期:2019年6月19日

雜志:Mol Syst Biol

文章在:https://www./doi/10.15252/msb.20188746

DOI:https:///10.15252/msb.20188746

圖1

2 摘要

單細胞領域日新月異,大量的工具被開發(fā)出來,但很難去判斷是否好用,而且如何組建一個分析流程是一個難點。本文將詳細介紹單細胞轉錄組數(shù)據(jù)分析的步驟,包括預處理(質(zhì)控、歸一化標準化、數(shù)據(jù)矯正、挑選基因、降維)以及細胞和基因?qū)用娴南掠畏治?。并且作者將整個流程應用在了一個公共數(shù)據(jù)集作為展示(詳細說明在:https://www.github.com/theislab/single-cell-tutorial),目的是幫助新入坑用戶建立一個知識體系,已入坑用戶更新知識體系。

3 前言

需要注意,雖然在原文鏈接中這些文獻鏈接可以打開,但會鏈接到原文的該文獻位置,而不是直接打開該文獻

現(xiàn)在已經(jīng)可以利用scRNA研究斑馬魚、青蛙、渦蟲的細胞異質(zhì)性(Briggs et al, 2018; Plass et al, 2018; Wagner et al, 2018) ,重新理解以前的細胞群體,但這個領域面臨的一個問題就是沒有成熟的標準化流程。標準化之路的困難有:大量分析方法和工具的誕生(截止2019.3.7 已經(jīng)有385種工具)、爆炸式增長的數(shù)據(jù)量(Angerer et al, 2017; Zappia et al, 2018)。另外根據(jù)不同研究目的,各種分支也突顯,例如在細胞分化過程中預測細胞命運(La Manno et al, 2018)。在我們眼界大開的同時,分析流程標準化就變得更加困難。

在未來分析流程標準化之路上,困難還會存在于技術整合層面。比如現(xiàn)在大量的scRNA工具都是用R和Python寫的,跨平臺分析需求在增長,而對編程語言的喜好也決定了工具的選擇。很多好用的分析工具將自己限制在用各自的編程語言開發(fā)的環(huán)境中,例如Seurat、Scater、Scanpy。

接下來,就一起看看作者列出了哪些他認為比較好的軟件和流程吧

先上一個scRNA分析總體流程圖:

圖2

4 預處理和可視化

4.1 首先看一下實驗過程

比較詳細的介紹可以看:Ziegenhain et al (2017); Macosko et al (2015); Svensson et al (2017).

原文描述的關鍵點是:

  • 4步走:Typical workflows incorporate single‐cell dissociation, single‐cell isolation, library construction, and sequencing.
    組織裂解=》細胞分離=》文庫構建=》測序

  • 第一步:As a first step, a single‐cell suspension is generated in a process called single‐cell dissociation in which the tissue is digested.

  • 第二步:To profile the mRNA in each cell separately, cells must be isolated. 寫了主要的2種方法:plate‐based、droplet‐based,當然也都存在一些問題: In both cases, errors can occur that lead to multiple cells being captured together (doublets or multiplets), non‐viable cells being captured, or no cell being captured at all (empty droplets/wells)

  • 第三步:Each well or droplet contains the necessary chemicals to break down the cell membranes and perform library construction. Furthermore, many experimental protocols also label captured molecules with a unique molecular identifier (UMI).

    UMI的作用主要是區(qū)分:UMIs allow us to distinguish between amplified copies of the same mRNA molecule and reads from separate mRNA molecules transcribed from the same gene.

  • 第四步:Libraries are labelled with cellular barcodes and pooled together (multiplexed) for sequencing.

感覺原文描述的還沒有illumina給出的詳細,那么就看看illumina的圖文并茂版:

illumina單細胞測序工作流程:關鍵步驟和注意事項:http://web./landing/products_view.asp?newsid=324

圖3
圖4

原始測序數(shù)據(jù)要經(jīng)過處理得到表達矩陣,注意這里有兩種表述方式:molecular counts (count matrices) 【也即是使用UMI的】和 read counts (read matrices),取決于是否使用UMI。而作者介紹的流程中,默認使用 count matrices,除非readmatrices和 count matrices得到的結果存在差異,才會特別介紹read matrices

關于比較read and molecule counts,有人寫了一個R流程:https://jdblischak./singleCellSeq/analysis/compare-reads-v-molecules.html

原始數(shù)據(jù)處理工具主要有:CellRanger、indrops、SEQC、zUMIs

它們主要做了這么幾件事:

  • read quality control (QC)
  • assigning reads to their cellular barcodes and mRNA molecules of origin (also called “demultiplexing”)
  • genome alignment
  • quantification

得到的矩陣行是轉錄本,列是barcodes【這里用barcodes而不是直接叫細胞,是因為不同細胞的reads也可能屬于同一個barcode =》如果出現(xiàn)一孔/液滴多細胞(doublet情況),那么barcode在多個細胞都是一樣的】當然也會出現(xiàn)有barcode但實際沒有細胞的情況(一個孔/液滴沒有細胞即droplet,但這個孔/液滴也會賦予barcode)

反正記?。篵arcode和孔/液滴是對應的,但一個孔/液滴中有一個細胞還是多個細胞或者沒有細胞,都會存在barcode。只不過最后可能會看到:多個細胞對應一個barcode、即使沒有細胞也會有barcode這樣的情況

關于10X實驗環(huán)節(jié),可以看我之前寫的:https://mp.weixin.qq.com/s/0DEybX7GnuDFhfY1uj9t9A

圖5

4.2 質(zhì)控

在正式分析之前,先要確定barcode是不是對應真正的細胞(上面已經(jīng)了解了barcode和細胞的關系),也就是進行Cell QC,主要考慮三個因素(這幾個因素也就是現(xiàn)在流程中常用的過濾指標):

  • the number of counts per barcode (count depth)
  • the number of genes per barcode
  • the fraction of counts from mitochondrial genes per barcode

從下面??圖中,感覺過濾的一個方向就是:保留大山頭,去掉小山頭(把略有增長但不礙大局的小山頭炸掉)

先看圖A:其中這個小的直方圖就是把count depth小于4000的放大,這里設定了一個閾值1500,也就是一個barcode中至少有1500的表達量

圖B:每個細胞中包含的基因數(shù)直方圖??梢钥吹綑M坐標有一個小的峰在400附近,這里設定的閾值是700

圖C:依舊是看count depth。從高到低排列count depth值,可以過濾一些空的液滴(empty droplets),看到從”肘部“也就是縱坐標1500左右開始迅速下降

圖D:看線粒體比例。如果占比很高并且細胞類型不是線粒體特別豐富的那種(如心肌細胞),可能說明這個細胞本身的基因數(shù)不多并且總體表達量也不高

圖6

以上三個指標固然重要,但如果只關注其中某一個,也會產(chǎn)生誤導作用,所以作者建議看問題一定要全面,并且要把數(shù)據(jù)和生物學知識結合起來。作者舉了個例子:比如線粒體表達量相對較高的細胞也可能參與了呼吸過程。細胞總體表達量低或者基因數(shù)量少,也可能是因為當時取的細胞處于靜止;細胞表達量很高,也可能因為本身細胞體積就比較大。的確,細胞與細胞之間的總表達量還是存在較大差異的。未來也許QC會提供更多的選擇。

除了檢查細胞完整度,QC還要進行轉錄本層面上的檢查。原始的count矩陣一般包含超過20000個基因。這里一般要根據(jù)在細胞中有表達的數(shù)量進行過濾,但這個閾值要根據(jù)總體細胞數(shù)和預計的分群情況來靈活調(diào)整。比如有的細胞類型本身就數(shù)量比較少(也許就50個),那么如果我們要設定”在少于50個細胞中有表達的基因“這種條件,那么可能會丟失那些總共就50個細胞中的marker基因,最終導致鑒定的細胞亞群會缺失。

質(zhì)控的目的就是給下游提供更高質(zhì)量的數(shù)據(jù),但一開始誰也不知道這個質(zhì)量高不高,只能先進行下游分析,看看結果(比如細胞分群結果)再判斷。尤其是針對異質(zhì)性高的細胞群體

文章又額外介紹了一些QC指標:
https://www./action/downloadSupplement?doi=10.15252%2Fmsb.20188746&file=msb188746-sup-0001-Appendix.pdf
比如在CellRanger的結果報告中就會有:Q30指標(一般要高于60-70%)、比對到外顯子的比例(一般認為比對到非外顯子區(qū)域超過40%就造成了測序的浪費)、看看實驗中用了多少個細胞,以及結果表達矩陣得到多個barcode

小結
  • 三種QC指標(the number of genes、the count depth 、the fraction of mitochondrial reads)要放在一起思考,而不是單獨看某一個
  • 先盡可能地設定寬泛的QC閾值,如果下游聚類無法解釋再回過頭來反思QC
  • 如果看到每個樣本的QC指標分布不同,那么就要對每個樣本分別設定閾值,而不是一刀切

4.3 歸一化/標準化

作為背景知識,首先來看:歸一化和標準化的區(qū)別(https://cloud.tencent.com/developer/article/1486102)但二者的界限也沒有特別明顯,也沒有必要把這兩個概念分的特別清楚。只要清楚它們大概的使用范圍就可以了:

  • 常用的歸一化是log處理,之前離散程度很大的數(shù)據(jù)就被集中了;
  • 常用的標準化是z-score:考慮到了不同樣本對表達量的影響,消除到了表達的平均水平和偏離度的影響

它們的使用范圍:

  • 如果對表達量的范圍有要求,用log歸一化
  • 如果表達量較為穩(wěn)定,不存在極端最大最小值,使用歸一化
  • 如果表達量離散程度很大,存在異常值和較多噪音,用標準化可以避免異常值和極端值的影響
  • 在分類、聚類、PCA算法中,使用z-score值的結果更好
  • 數(shù)據(jù)不太符合正態(tài)分布時,可以使用歸一化
  • 機器學習的算法(SVM、KNN、神經(jīng)網(wǎng)絡等)要求歸一化/標準化

繪制熱圖會經(jīng)常用到z-score去除極端值

pheatmap(dat) # scale之前n=t(scale(t(dat)))n[n>2]=2 # 限定上限n[n< -2]= -2 # 限定下限pheatmap(n,show_colnames =F,show_rownames = F) # scale之后

接著:在單細胞分析中,也會同時用到Normalize和Scale(可以看:單細胞Seurat包升級之2,700 PBMCs分析

  • 歸一化 Normalize做的就是將數(shù)據(jù)進行一個轉換,可以讓同一基因在不同樣本中具有可比性(例如RPKM、TPM等);另外降低離散程度??词褂玫暮瘮?shù)LogNormalize背后的計算方法就是:log1p(value/colSums[cell-idx] *scale_factor) ,它同時考慮到了這兩點
  • 標準化Scale就是基于之前歸一化的結果(也就是log后的結果),再添z-score計算

最后,在對細胞文庫差異進行normalization 這一篇中也提到了:

  • Normalization 'normalizes' within the cell for the difference in sequenicng depth / mRNA thruput
  • Scaling 'normalizes' across the sample for differences in range of variation of expression of genes
  • normalization一般是對文庫處理,目的消除一些技術差異;scale一般對基因表達量處理(典型的z-score:表達量減均值再除以標準差),目的是后續(xù)分析不受極值影響

表達矩陣中的每個count值都表示成功的細胞捕獲、成功的反轉錄、成功的測序。但即使是相同類型的細胞,它們的count depth(也就是每個細胞的全部表達量)也會有變化,變化的來源就在于上面說的那三步。因此在比較兩個細胞時,任何差異都可能由于實驗測序誤差產(chǎn)生,而不是真的生物學差異。歸一化就是解決這個問題,它把要比較的兩個count值根據(jù)各自身處的環(huán)境求出一個相對豐度,也就是放在了一個水平上考慮,減少實驗測序誤差,突出更多的生物學差異。

最常用的歸一化方法就是:count depth scaling,也稱為counts per million(CPM),這個方法常用于bulk轉錄組,它會根據(jù)每個細胞的總表達量計算一個 size factor ,然后對其中各個基因表達量進行normalize。

這里再回顧下其他一些方法:跟著豆豆一起回顧標準化方法
另外來自:https://cloud.tencent.com/developer/article/1484078

  • RPM沒有考慮轉錄本的長度的影響。適合于產(chǎn)生的read讀數(shù)不受基因長度影響的測序方法,比如miRNA-seq測序,miRNA的長度一般在20-24個堿基之間
  • RPKM/FPKM考慮了轉錄本的長度的影響。適用于基因長度波動較大的測序方法,如lncRNA-seq測序,lncRNA的長度在200-100000堿基不等
  • TPM是先去除了基因長度的影響,而RPKM/FPKM是先去除測序深度的影響。TPM實際上改進了RPKM/FPKM方法在跨樣品間定量的不準確性。

單細胞測序中使用的歸一化方法由于細胞種類和基因錯綜復雜,有人就在bulk的基礎上進行了改動。例如:Weinreb et al (2018) 先排除了表達量超過總體5%的基因,然后再計算size factor,主要是預防少量極高表達量基因的存在;Scran包有個pooling‐based size factor estimation方法,允許更高的細胞異質(zhì)性存在;另外Scran包在批次矯正和差異分析環(huán)節(jié)也比其他歸一化方法表現(xiàn)更好(Buttner et al, 2019)。

在單細胞RNA測序領域,目前有三種常用方法:其一是以10x Genomics為代表的微滴(droplet-based)測序;其二是以Namocell為代表的PCR板(plate-based)測序;其三是以BD Rhapsody為代表的微孔(micro-well-based)測序。就測序長度來說,Smart-seq/C1和Smart-seq2基于full length的測序方案,CEL-seq2, Drop-seq, MARS-seq, SCRBseq是基于UMI的測序方案。

不能指望某一種方法適用于所有類型的scRNA數(shù)據(jù),(Cole et al, 2019)就發(fā)現(xiàn)不同的歸一化方法對于不同類型數(shù)據(jù)集表現(xiàn)不同,使用scone工具可以幫助選擇合適的方法。

一般在歸一化后,數(shù)據(jù)都會變成log(x+1)的樣子,但之后是否對基因進行z-score的標準化上,沒有一個共識。Seurat的教程基本都使用了scale這一步,但Slingshot的作者就反對對基因進行scale (Street et al, 2018)。在本文中,作者傾向于避免對基因進行scale。

使用log轉換的一個好處就是:讓數(shù)據(jù)更加集中,減少數(shù)據(jù)的偏斜度,從而近似于許多下游分析工具對數(shù)據(jù)為正態(tài)分布的假設(盡管scRNA數(shù)據(jù)并不是真正的符合正態(tài)分布),比如在差異表達分析和批次矯正環(huán)節(jié)

小結
  • 對于非全長scRNA數(shù)據(jù)(如10X),推薦使用scran的歸一化方法;
  • 對于plate‐based 的數(shù)據(jù),可以用scone工具來進行評價,進而可以更好地處理plate之間的批次效應(Cole et al, 2019);
  • 對于全長scRNA數(shù)據(jù)(如smart-seq)可以借用bulk的方法(如TPM),來矯正基因長度【這個問題在10X中不存在:in 10x single cell 3' or 5' gene expression assay, this gene-length bias does not exist. https://kb./hc/en-us/articles/115003684783-How-to-calculate-TPM-RPKM-or-FPKM-instead-of-counts-
  • 歸一化的數(shù)據(jù)應該是log(x+1)這種形式的
  • 是否進行scale沒有共識,這里作者不推薦scale

4.4 數(shù)據(jù)矯正與整合

數(shù)據(jù)矯正的對象種技術和生物因素都有,例如:不同批次、捕獲失?。╠ropout)、不同細胞周期。這些在之前的歸一化中沒有被矯正,但這些差異因素都可能會后面的分析產(chǎn)生影響,它們現(xiàn)在都是導致差異的”嫌疑人“之一。這里要做的就是把這些差異來源去掉(Regressing out 《=》【專門查的詞典】 同義詞partialling out :剔除)

4.4.1 首先是生物因素

最常見的生物矯正因素就是:轉錄組中的細胞周期信息。簡單一點的方式就像Scanpy和Seurat對細胞周期評分進行簡單線性回歸;復雜點的方式就像scLVM和f‐scLVM。用來計算細胞周期分數(shù)的marker基因可以從文獻中獲得 (Macosko et al, 2015)。另外,這些方法還能用來去除其他已知的生物因素,例如線粒體基因表達量(可以作為細胞應激的標記)。

需要注意的是:

  • 細胞周期因素并非一無是處,例如在一個增殖的細胞群中,所有細胞不是同步增殖的,那么就可以根據(jù)細胞周期評分來識別。所以需不需矯正還要根據(jù)研究目的判斷
  • 需要結合具體分析的生物問題來判斷是否去除。生物體的多個生物過程往往存在依賴性,因此矯正其中一個過程,可能無意間掩蓋了另一個過程
  • 有人認為細胞大小的變化和細胞周期有關 (McDavid et al, 2016),因此在歸一化過程中對細胞大小進行矯正,或者使用專用的工具如cgCorrect,也可以部分修正細胞周期的影響
4.4.2 然后是技術因素

最常見的技術矯正因素就是:樣本測序深度、批次、噪音。

去除測序深度的影響,可以促進軌跡推斷算法的表現(xiàn),因為它需要在細胞之間找變化的路徑,只要放在同一水平才能看到更準確的總體表達高低。

批次的來源可能是:細胞捕獲的時期不同、文庫制備使用的芯片不同、測序使用的lane不同。由此產(chǎn)生的效應存在于多個層面:一次實驗中各個細胞群之間、同一實驗室中進行的不同實驗之間、或來自不同實驗室的數(shù)據(jù)集之間。這里主要介紹第一種和最后一種情況:

  • 第一種:一次實驗中各個細胞群之間是最經(jīng)典的情形,在bulk轉錄組也是常見的。使用線性方法進行矯正。例如ComBat工具就是利用線性矯正
  • 最后一種:來自不同實驗室的數(shù)據(jù)集之間的”數(shù)據(jù)整合“。使用非線性方法進行矯正。例如:CCA(Canonical Correlation Analysis)、MNN(Mutual Nearest Neighbours )、Scanorama、RISC、scGen、LIGER、BBKNN、Harmony。雖然這些非線性矯正方法也能用于第一種經(jīng)典的批次情形,但可能會由于自由度增加而導致矯正過度。例如,在第一種經(jīng)典模式下,Combat表現(xiàn)就比MNN好 (Buttner et al, 2019)

看一下Combat矯正前后的差別:其中顏色表示不同樣本

圖7

去噪也是矯正的一種類型。單細胞數(shù)據(jù)的一個特點就是含有許多噪音來源,其中一個就是dropout。一些工具就用來推斷dropout,用適當?shù)谋磉_量來替代0,例如:MAGIC、DCA、scVI、SAVER、scImpute。去噪可以提高基因間相關性的估計。這一步可以和歸一化、批次矯正及其他下游分析整合起來,例如基于Python的scVI工具。但任何方法都可能導致矯正過度或不足。

4.4.3 小結
  • 判斷要不要進行矯正生物因素:主要看后續(xù)分析是不是用于研究發(fā)育軌跡等特定生物過程
  • 技術因素和生物因素需要放在一起矯正,而不是先矯正這個,后矯正那個
  • plate-based數(shù)據(jù)的預處理一般需要利用非線性歸一化方法
  • 需要關注降噪前表達量為0,而降噪后才有表達的基因

4.5 挑選基因、降維、可視化

人類的scRNA數(shù)據(jù)中可能會包含25000個基因,但其中許多基因并非能提供有用信息,還有很多基因表達量直接為0。即使在QC階段去掉這些表達量為0的基因,一個單細胞數(shù)據(jù)的基因空間依然會有超過15000個維度(一個基因表示一個維度),因此需要降低維度

4.5.1 首先挑選基因

就是挑那些真正”具有情報價值“的基因,也就是會數(shù)據(jù)變化起作用的基因。因此我們這里會挑選名為HVG的基因,也就是highly variable genes。根據(jù)數(shù)據(jù)集的復雜程度不同,HVGs一般會有1000-5000個(如下圖就對不同數(shù)據(jù)集的HVGs做了個統(tǒng)計)

圖8

之前有研究表明,HVGs數(shù)量從200到2400,它們降維后的表現(xiàn)差不多(Klein et al (2015),作者建議先盡量多選一些HVGs。

比較流行的挑選HVGs的方法有Scanpy和Seurat,而且最好是在去除技術因素后挑選,避免因為批次、測序等因素導致錯誤挑選HVG。當然還有其他挑選的方法,看Yip et al (2018).

4.5.2 接著降維

挑出來HVGs后,就是降維了,力求在最少的維度中捕捉到最多的數(shù)據(jù)特征。

常用的降維方法:A-F分別是:PCA、t-SNE、diffusion maps、UMAP、ForceAtlas2(force‐directed graph)、Variance explained by the first 31 principal components (PCs)。關于單細胞數(shù)據(jù)的降維方法,詳細可以看:Moon et al (2018)

圖9

其中兩個應用比較廣的方法是:PCA(Pearson, 1901)和diffusion maps (Coifman et al, 2005) 【diffusion maps 于2015年在單細胞領域走紅 Haghverdi et al (2015) 】

  • 主成分分析PCA是一種線性方法,通過最大化每個其他維度中捕獲的殘差來生成縮減的維度。而且,PCA常作為非線性降維方法的預處理手段。PCA一般通過前N個主成分來表示整個數(shù)據(jù)集,其中N可以用F中的”肘elbow“部判斷,或者用基于置換檢驗的jackstraw方法確定
  • diffusion maps 是非線性的方法,它強調(diào)數(shù)據(jù)之前的轉換。當研究連續(xù)型數(shù)據(jù)例如感興趣的分化過程時會使用。它的每個成分(component)都能突出不同類型細胞間的異質(zhì)性
4.5.3 最后可視化

可視化一般使用非線性降維的方法。最常用的就是2008年提出的t-SNE( t‐distributed stochastic neighbour embedding)。t-SNE的一個特性就是關注局部而忽視整體,因此帶來的一個影響就是:可視化結果可能夸大了細胞群之間的差異,忽略了這些細胞群之間的潛在聯(lián)系

另外,使用t-SNE的一大難點就是perplexity參數(shù)的設定,因為這個數(shù)不同,結果顯著的cluster數(shù)也會不同 (Wattenberg et al, 2016)。

除了t-SNE,還有2018年推出的UMAP和SPRING可以用,在缺乏明確的生物學問題時,可以用UMAP作為不錯的數(shù)據(jù)探索。

小結
  • 根據(jù)數(shù)據(jù)集的復雜性,推薦選擇1000-5000個HVGs
  • 推薦UMAP進行數(shù)據(jù)探索;PCA獲得一般性數(shù)據(jù)總結; diffusion maps作為PCA的替代,可用于軌跡推斷
  • PAGA方法與UMAP連用適用于特別復雜的數(shù)據(jù)集

4.6 「總結」 預處理的各個階段

作者貼心將預處理比作5種類型數(shù)據(jù)的處理:

原始數(shù)據(jù)(raw data)、歸一化數(shù)據(jù)(normalized data)、矯正后的數(shù)據(jù)(corrected data)、挑選后的數(shù)據(jù)(feature‐selected data)、降維后的數(shù)據(jù)(dimensionality‐reduced data)

這5個階段又分成3個層次:

  • measured data:用于統(tǒng)計檢驗
  • corrected data:用于數(shù)據(jù)比較可視化
  • reduced data:用于下游分析

其中每個步驟適時調(diào)整,例如單一批次的數(shù)據(jù)集,就可以跳過矯正批次這一步

圖10

5 下游分析之細胞層面

下游分析的目的是解釋生物問題,例如根據(jù)表達量將細胞劃分成不同的類型;相似細胞間表達量的微小變化也會體現(xiàn)連續(xù)的分化路徑;基因表達量之間的相關性可能與基因共表達有關...

下游分析也是有細胞層面和基因?qū)用妫?/p>

  • 細胞層面主要關注:分出多少群細胞、細胞的軌跡。細胞類型為了解釋異質(zhì)性的問題;軌跡作為一個動態(tài)發(fā)育過程中的一個”快照“,可以幫助理解某個動態(tài)發(fā)育過程
  • 基因?qū)用婢褪牵翰町惙治觥⒏患治?、互作網(wǎng)絡
圖11

下面??先看看看細胞層面的分析之分群和軌跡

5.1 細胞分群

5.1.1 先是:分群方法

這里主要都是算法相關,簡單了解即可

將細胞分群基本就是任何單細胞分析的必經(jīng)之路。群的劃分就是根據(jù)細胞中基因表達譜的相似性,表達譜的相似性是由于歐幾里得距離量度決定的,而距離量度又是利用的降維的數(shù)據(jù)。一般有兩種方法計算:clustering algorithms、community detection methods

  • clustering algorithms是直接基于距離的經(jīng)典無監(jiān)督機器學習方法。最常見的是k-means。k-means使用的空間距離量度也有不同(默認是歐氏距離),比如cosine similarity、 correlation‐based distance metrics、the SIMLR method。最近的研究表示correlation‐based distances優(yōu)于其他距離 (Kim et al, 2018)。

  • community detection methods是graph‐partitioning algorithms,利用了K近鄰法 K‐Nearest Neighbour approach (KNN graph)作圖。細胞就是圖上的一個節(jié)點,每個細胞都和K個最相似的細胞連接,這個相似性也是根據(jù)降維后空間中的歐氏距離計算的。根據(jù)數(shù)據(jù)集的大小,K一般設為5-100

    community detection methods一般比clustering algorithms速度快,因為只有相鄰的細胞對被認為屬于相同的群,大大減少了可能的細胞群的搜索范圍

5.1.2 然后是:分群后的注釋

這個過程主要是基因?qū)用娴牟僮?,為每個cluster找marker gene(也就是能代表這個cluster的基因,而這個基因又和已知的細胞類型有關)。任何的分群算法和參數(shù)設置都會將一整團細胞分成多個群,但這些群是否真的有意義,就要靠這一步來和生物背景結合起來。

我們希望看到的是存在很多類型的細胞,來說明細胞異質(zhì)性的問題,但這里關于細胞類型這個定義還是存在爭議。首先,細胞類型的劃分怎樣算是清楚,對于一些人來說,”T cells“這個名稱可以叫一個細胞類型,但還有人認為,必須繼續(xù)深入,像”CD4+ T cells“、”CD8+ T cells“才叫細胞類型;另外,即使是同一種細胞類型的細胞也會有不同的發(fā)育狀態(tài),因此它們也會顯示不同的分群結果。但不管如何,它們都是當時細胞的一種身份(identity)

這個很好理解,就像人一樣,人生階段不同身份也不同,但不能簡單說它的類型發(fā)生了變化。

因此,我們將分群的結果稱為不同身份的細胞(cell identities)會比不同類型的細胞(cell types)要好一些【即每個亞群可能并不是真的不同類型細胞,只是顯示了此時此刻的細胞身份

對于不同細胞身份的注釋,近年來也隨之細胞圖譜的研究而加速,例如小鼠腦細胞圖譜 (Zeisel et al, 2018) 、人類細胞圖譜 (Regev et al, 2017)的發(fā)現(xiàn),產(chǎn)生了許多參考數(shù)據(jù)庫。在缺乏相關背景的情況下,我們可以借用數(shù)據(jù)庫中已發(fā)現(xiàn)的細胞marker 基因套入我們的細胞,幫助判斷細胞身份。需要注意:通常使用的細胞表面marker基因在細胞身份鑒定方面存在局限性(Tabula Muris Consortium et al, 2018)

看這個注釋結果:
圖A是利用Louvain方法分群+UMAP可視化;
圖B是細胞身份的鑒定:stem cells (Slc12a2), enterocytes (Arg2), goblet cells (Tff3) and Paneth cells (Defa24)。但要注意marker基因可能也會在其他身份的細胞中表達,例如很多marker都在右上角(goblet and Paneth 細胞)中有表達,但最后還是根據(jù)表達量來指定特定的細胞身份(比如Slc12a2基因雖然在很多細胞都表達,但就是在中部偏右這一坨細胞中表達量相對高,所以把它當做stem cell)
圖C是近端(上圖)和遠端(下圖)腸上皮區(qū)域的細胞身份組成圖(顏色越深細胞密度越大)

圖12

上面提到根據(jù)marker基因進行細胞分群注釋。那么marker基因怎么獲得?

利用差異分析,分成兩組:某個cluster中的細胞、數(shù)據(jù)集中其余全部的細胞。然后重點關注這個cluster中上調(diào)的基因,因為marker基因一般具有更強的表達作用。差異分析也會使用簡單的統(tǒng)計檢驗,例如Wilcoxon rank‐sum test、t-test,將基因的差異大小排個序,選出排名靠前的基因來作為marker基因

有了marker基因,再進行注釋

將數(shù)據(jù)集中選出的marker基因和參考數(shù)據(jù)集進行比對,統(tǒng)計方法可以是:enrichment tests、the Jaccard index、other overlap statistics

參考數(shù)據(jù)集可以是網(wǎng)頁工具: www.mousebrain.org、 http:///,可以將選出的marker基因在參考數(shù)據(jù)集中進行可視化,幫助判斷這個marker基因是什么細胞身份

注釋并非一蹴而就,這個很麻煩...

細胞分群、分群注釋、重分群、重注釋...這個循環(huán)很耗費時間。自動化注釋方法加快了這個過程,例如scmap (Kiselev et al, 2018b) 、Garnett (preprint: Pliner et al, 2019) ,但這樣的方法有利有弊。自動化提高了速度,但相比手動注釋也降低了靈活性。畢竟自動化工具使用的參考數(shù)據(jù)集中可能并不包含我們數(shù)據(jù)中的這樣細胞。因此,有自動化工具也不能完全拋棄手動挑選,尤其針對大型數(shù)據(jù)集中多種多樣的細胞。自動化的過程可以先幫我們粗略地給細胞加個標記,如果有需要,我們可以繼續(xù)手動對這種細胞繼續(xù)劃分子細胞。對于小型數(shù)據(jù)集或者缺乏參考基因集的,手動注釋就足夠了。

5.1.3 注意
  • 同一細胞身份的marker基因在不同數(shù)據(jù)集之間可能由于數(shù)據(jù)集細胞類型和狀態(tài)組成而不同【選出的marker基因并不是說以后遇到它,就一定等同于這種類型的細胞。只是說在某種細胞的某個狀態(tài)下,這個marker基因更符合
  • 如果存在參考數(shù)據(jù)集(例如 www.mousebrain.org、 http:///),建議輔助自動化工具進行注釋,減少手動查基因的時間
5.1.4 細胞分群衍生——細胞組成分析(Compositional analysis)

就像上面的圖12中的C圖,顯示的是近端(上圖)和遠端(下圖)腸上皮區(qū)域的細胞身份組成圖(顏色越深細胞密度越大)。研究細胞組成的變化也是一個新方向,例如沙門氏菌感染已被證明會增加小鼠腸上皮細胞的比例 (Haber et al, 2017)。

這個分析既需要足夠多的細胞數(shù)量來推斷各個cluser的占比,又需要足夠的樣本數(shù)量來證明是單純一個樣本得cluster數(shù)量這樣變還是總體都會這樣變。相關的分析工具還沒有太多,未來的開發(fā)可能會借鑒單細胞質(zhì)譜流式(mass cytometry)或者是宏基因組分析【單細胞與宏基因組的結合...】

5.2 軌跡分析

5.2.1 軌跡推斷Trajectory inference

軌跡推斷就是為了找到不同細胞身份、分化或者生物過程中漸進式非同步的變化,構建出的一個動態(tài)模型。它認為單細胞數(shù)據(jù)實際上就是一個連續(xù)過程中的快照(snapshot),這個過程可以通過在細胞空間中尋找最小化相鄰細胞間轉錄變化的路徑來重建

例如:
圖A就是利用Slingshot推斷近端(proximal)和遠端(distal)腸上皮細胞的分化軌跡;
圖B就是在PCA空間中進行的Slingshot推斷。圖中細胞的路徑就叫做”擬時序“(pseudotime)

圖13

2014年Monocle和Wanderlust先推出了軌跡推斷,之后誕生的分析方法更加豐富,它們在建模路徑的復雜性上有所不同,從簡單的linear or bifurcating(分叉) trajectories,到復雜的graphs, trees, or multifurcating(多叉) trajectories。Saelens et al, 2018)進行過軌跡推斷方法的比較,結論是沒有一種方法對所有類型的軌跡推斷有效,應該根據(jù)預期軌跡的復雜度來選擇。不過,Slingshot在簡單軌跡推斷中優(yōu)于其他方法(Street et al, 2018) 。如果期望得到更復雜的軌跡,PAGA值得推薦。軌跡推斷是一個不確定的過程,可以用多種方法來進行佐證。

  • 細胞內(nèi)通常會同時發(fā)生多個生物學過程,因此在進行發(fā)育軌跡推斷時,可以將其他生物因素去掉,例如T細胞在逐漸成熟的過程中就可能會經(jīng)歷細胞周期轉變(Buettner et al, 2015)。

  • 另外軌跡推斷最好是在細胞分群之后進行,因為一個cluster的形成可能意味著這一坨細胞處于比較穩(wěn)定的狀態(tài)了。

  • 此外,RNA速率(RNA velocities)可以添加發(fā)育軌跡的方向,例如:scVelo

    RNA速率
  • 當然,推斷的軌跡不一定就代表一個生物學過程,因為畢竟是根據(jù)”快照“數(shù)據(jù)中的轉錄狀態(tài)推測的。后續(xù)可以借鑒:perturbation experiments、inferred regulatory gene dynamics、support from RNA velocity

5.2.2 基因表達量的動態(tài)變化

在擬時序(pseudotime)中變化的基因描述了軌跡,這組與軌跡相關的基因有望包含調(diào)控建模過程的基因,可以用來識別潛在的生物過程。

目前很少有專門分析基因表達動態(tài)變化的工具。BEAM將Monocle的軌跡推斷整合進來,允許檢測在軌跡分支過程中相關基因的動態(tài)變化。另外還有LineagePulse (https://github.com/YosefLab/LineagePulse)考慮了dropout技術噪音但還在開發(fā)中。

下面這樣的圖在Slingshot的幫助文檔就有提及:https:///packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html 【4.1:Identifying temporally expressed genes】

require(gam)t <- sim$slingPseudotime_1# for time, only look at the 100 most variable genesY <- log1p(assays(sim)$norm)var100 <- names(sort(apply(Y,1,var),decreasing = TRUE))[1:100]Y <- Y[var100,]# fit a GAM with a loess term for pseudotimegam.pval <- apply(Y,1,function(z){    d <- data.frame(z=z, t=t)    suppressWarnings({      tmp <- suppressWarnings(gam(z ~ lo(t), data=d))    })    p <- summary(tmp)[3][[1]][2,3]    p})topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:100]heatdata <- assays(sim)$norm[topgenes, order(t, na.last = NA)]heatclus <- sim$GMM[order(t, na.last = NA)]heatmap(log1p(heatdata), Colv = NA,        ColSideColors = brewer.pal(9,'Set1')[heatclus])
Slingshot基因表達量的動態(tài)變化
5.2.3 細胞亞穩(wěn)態(tài)分析 Metastable states

亞穩(wěn)態(tài)常見于物理化學。在物理學中,亞穩(wěn)性(Metastable)是動力系統(tǒng)的一種穩(wěn)定狀態(tài),而不是系統(tǒng)能量最低的狀態(tài)
Metastable:stable provided it is subjected to no more than small disturbances.
另外這個狀態(tài)可以用這個圖幫助理解:大概就是「一個相對穩(wěn)定但又會變化的一個狀態(tài)」


擬時序分析會展示出不同階段細胞數(shù)量的多少。假設細胞以無偏的方式采樣,其中軌跡中的稠密區(qū)域就表示轉錄時首選的方案。當把軌跡理解為一條時間線時(例如在發(fā)育這個時間線),這些密集的區(qū)域可能代表細胞的亞穩(wěn)態(tài),可以結合擬時間坐標來繪制直方圖,找到這些亞穩(wěn)態(tài)【因此看到B圖中很多種狀態(tài),但C中直方圖認為這幾個密集的區(qū)域才屬于亞穩(wěn)態(tài)】

Metastable states
5.2.4 整合分群與軌跡分析

分群是由整體到部分,是靜態(tài)的;而軌跡又是由部分推斷整體,是動態(tài)的。二者結合起來又產(chǎn)生了一種新的分析模式

將分群的結果當成節(jié)點(node),將軌跡當成節(jié)點之間的橋梁(edge),所以將動靜數(shù)據(jù)結合在了一起。利用partition‐based graph abstraction(PAGA)這個工具就能得到類似下面這個圖。

It was the only reviewed method able to cope with disconnected topologies and complex graphs containing cycles.

整合分群與軌跡分析

6 下游分析之基因?qū)用?/h3>

之前都是對細胞進行分析,但細胞中的基因分析會提供更多的信息。例如差異表達分析、基因集分析和基因調(diào)控網(wǎng)絡推斷,不是表面上研究細胞異質(zhì)性,而是基于異質(zhì)性探索基因表達相關的原因

6.1 差異表達分析

基因?qū)用娴臄?shù)據(jù)探索,一個經(jīng)常遇到的問題就是:兩個組之間有沒有表達量的差異?

這個方法也是常規(guī)bulk轉錄組中經(jīng)常做的。不過單細胞相比于bulk轉錄組的一個優(yōu)勢就是:可以深入一個層次,原來bulk只是看一塊組織的平均表達量,但現(xiàn)在經(jīng)過分群后,能得到一塊組織中各種各樣的亞群,再結合差異分析,對理解異質(zhì)性問題更有幫助。

雖然都是朝著一個方向前進,但單細胞和bulk轉錄組的差異分析方法還是不同的。

  • bulk轉錄組存在樣本數(shù)量的限制,因此算法需要對少量樣本進行準確估計,而單細胞則不同,一個細胞作為一個樣本,成百上千不在話下;
  • 單細胞數(shù)據(jù)又有自己的特點:特異的人為技術噪音(dropout、high cell‐to‐cell variability ),因此單細胞分析方法需要額外考慮這些因素

但最近(Soneson & Robinson, 2018)研究表明,基于大批量的差異分析,bulk分析方法的性能與最好的單細胞分析方法相當。當bulk方法進行改進,加入基因權重分析后,表現(xiàn)要好于單細胞原有工具。例如:bulk差異分析工具DESeq2/EdgeR + ZINB‐wave工具估算的權重。

不過,bulk差異分析工具的性能雖然好,但是計算的效率很難提升。畢竟單細胞數(shù)據(jù)樣本數(shù)量越來越多,程序跑的時間長短也成了衡量工具優(yōu)劣的重要因素。單細胞工具MAST脫穎而出。在單個數(shù)據(jù)集的小范圍比較中,完勝bulk和其他單細胞方法(Vieth et al, 2017)。而且MAST比bulk方法快了10到100倍 (Van den Berge et al, 2018) 。

小結
  • 差異分析使用MAST或limma
  • 差異分析不能使用矯正后的數(shù)據(jù)(denoised, batch corrected, etc.),而是應該在計算過程中去指定需要矯正的技術因素
  • 我們給差異分析算法提供的矯正的因素(稱之為協(xié)變量covariates)不能太混亂,因為工具不會去智能識別,必須要清楚需要矯正什么

6.2 基因集分析

基因?qū)用娴姆治觯鶗a(chǎn)生大量的基因,但很難去解釋。

例如差異分析我們往往能得到上千基因,為了比較方便解讀,一般會把有共同特性的基因歸為一組,然后檢查我們歸類的可靠性 【grouping the genes into sets based on shared characteristics and testing whether these characteristics are overrepresented in the candidate gene list.】

我們一般關注基因在生物過程(biological processes, BP)中的富集,可以使用MSigDB、GO、KEGG pathway、Reactome數(shù)據(jù)庫

另外,單細胞中的一個新進展就是利用成對基因標簽進行配體受體分析( ligand–receptor analysis)

來自:https://www./showarticle.asp?id=453107210
腫瘤內(nèi)細胞-細胞相互作用的研究將通過對配體和受體的表達分析來探究細胞間的相互交流。運用配體-受體復合物的數(shù)據(jù)庫,通過scRNA-seq數(shù)據(jù)與腫瘤細胞亞群的定義相結合,來推斷潛在的細胞-細胞相互作用,可以理解為其中一個群體產(chǎn)生配體,向另一個表達相應受體的群體發(fā)信號

配體-受體成對標簽可以從:CellPhoneDB數(shù)據(jù)庫獲得,然后用來解釋cluster之間高表達基因的聯(lián)系

例如,利用celltalker 就可以做

Celltalker分析

6.3 基因調(diào)控網(wǎng)絡 gene regulatory network (GRN)

基因并非獨立發(fā)揮作用的。相反,基因的表達水平取決于與其他基因和小分子之間的相互調(diào)控

方法例如:SCONE、PIDCSCENIC (Single-Cell rEgulatory Network Inference and Clustering),但發(fā)展還不是很完善,推斷的調(diào)控關系不是很穩(wěn)定【謹慎使用】

7 分析平臺

現(xiàn)在開發(fā)了很多平臺,整合了一套分析流程,有基于R的(McCarthy et al, 2017; Butler et al, 2018) ,python的 (Wolf et al, 2018),本地的(Patel, 2018; preprint: Scholz et al, 2018) ,網(wǎng)頁版帶可視化的(Gardeux et al, 2017; Zhu et al, 2017)

Zhu et al (2017) and Zappia et al (2018).列出了各種平臺

Seurat是使用最廣泛的,Scater在QC和預處理中表現(xiàn)優(yōu)異;除此以外,基于Python的scanpy也逐漸發(fā)展起來,它對于大量細胞的標準化方面表現(xiàn)不錯

如果不使用命令行,可視化界面也有,只不過用戶只能跑人家已經(jīng)寫好的腳本,操作靈活性不足。這樣的平臺更多的用處是在可視化探索上,例如Granatum、ASAP。未來 Human Cell Atlas(HCA)會在數(shù)據(jù)可視化探索上迅速發(fā)展: https://www./data-sharing

8 結語

8.1 作者的結語

作者把流程測試和說明都放在了:https://github.com/theislab/single-cell-tutorial

感興趣的可以跟著走一遍,比較一下不同的工具。作者希望這一篇能代表單細胞領域目前發(fā)展的一個最新動向。他也提到,新方法層出不窮,本文介紹的大量的方法是經(jīng)過實踐比較、驗證過的。目前可用的方法不管是運行效率還是易用性可能都不如最新開發(fā)的方法,但要注意:新方法在未被大量驗證之前都需小心使用。而且新方法一般都是針對單個層面(比如降維、分群、軌跡推斷等),大體的分析流程基本固定了。

未來整合深度學習和單細胞多組學是兩個重要的發(fā)展方向,流程化運行更是趨勢。

隨著文庫制備和測序技術的進步,未來的單細胞平臺必將可以處理多種類型數(shù)據(jù):DNA甲基化、蛋白豐度等等。

8.2 劉小澤的結語

截止到2020年5月8日下午15.35,打卡看完!

三天的時間,基本每天都會花半天時間在閱讀這篇綜述上。從第一眼看到它的文章邏輯,就感覺:嗯是它,沒錯了!連午覺都不想睡了。

一開始想強迫自己看下去,沒想到,越看越精彩。尤其是將整個流程和自己的知識結合起來,就看得比較順暢。為了更加易讀,我在其中加了很多注釋,包括之前自己寫的一些推文和網(wǎng)上一些好的資源,可以幫助梳理知識點。

最后,希望看完本文對你有幫助??!


歡迎關注我們的公眾號~_~  
我們是兩個農(nóng)轉生信的小碩,打造生信星球,想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!

    本站是提供個人知識管理的網(wǎng)絡存儲空間,所有內(nèi)容均由用戶發(fā)布,不代表本站觀點。請注意甄別內(nèi)容中的聯(lián)系方式、誘導購買等信息,謹防詐騙。如發(fā)現(xiàn)有害或侵權內(nèi)容,請點擊一鍵舉報。
    轉藏 分享 獻花(0

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多

    中文字幕在线区中文色| 国产精品成人免费精品自在线观看| 黄色片一区二区在线观看| 91人妻人人澡人人人人精品| 国产精品福利精品福利| 亚洲综合色婷婷七月丁香| 蜜桃传媒视频麻豆第一区| 免费啪视频免费欧美亚洲| 亚洲欧美日韩另类第一页| 日韩偷拍精品一区二区三区| 99久热只有精品视频免费看| 国产日韩欧美在线播放| 91精品国产品国语在线不卡| 免费在线观看激情小视频| 免费高清欧美一区二区视频 | 1024你懂的在线视频| 国产一区二区三区香蕉av| 成人精品视频在线观看不卡| 亚洲av日韩av高潮无打码| 国产熟女一区二区精品视频| 亚洲国产成人一区二区在线观看| 东京热男人的天堂久久综合| 亚洲国产另类久久精品| 99久久国产精品免费| 亚洲综合一区二区三区在线| 好骚国产99在线中文| 欧美小黄片在线一级观看| 精品日韩中文字幕视频在线| 国产日产欧美精品视频| 在线一区二区免费的视频| 亚洲成人免费天堂诱惑| 蜜桃av人妻精品一区二区三区 | 亚洲性生活一区二区三区| 欧美日韩高清不卡在线播放| 中文字幕乱码一区二区三区四区| 久久精品中文字幕人妻中文| 亚洲日本久久国产精品久久| 熟女体下毛荫荫黑森林自拍| 免费在线观看激情小视频| 成年女人午夜在线视频| 91精品蜜臀一区二区三区|