作者，Evil Genius
標(biāo)準(zhǔn)的單細(xì)胞rna測序分析(scRNA-seq)工作流程軟件括通過序列排列將原始讀取數(shù)據(jù)轉(zhuǎn)換為細(xì)胞基因計數(shù)矩陣，然后進(jìn)行過濾、高變量基因選擇、降維、聚類和差異表達(dá)分析等分析。Seurat和Scanpy是實(shí)現(xiàn)這種工作流的最廣泛使用的軟件，通常被認(rèn)為是實(shí)現(xiàn)類似的單個步驟。
下面我們就需要比較一下軟件之間、以及不同版本之間的數(shù)據(jù)分析差異。
單細(xì)胞rna測序(scRNA-seq)是一種強(qiáng)大的實(shí)驗(yàn)方法，為基因表達(dá)分析提供細(xì)胞分辨率。隨著scRNA-seq技術(shù)的廣泛應(yīng)用，分析scRNA-seq數(shù)據(jù)的方法也越來越多。然而，盡管已經(jīng)開發(fā)了大量的工具，但大多數(shù)scRNA-seq分析都是在兩種分析平臺之一進(jìn)行的:Seurat或Scanpy。表面上，這些程序被認(rèn)為實(shí)現(xiàn)了分析相同或非常相似的工作流程:scRNA-seq結(jié)果計算分析的第一步是將原始讀取數(shù)據(jù)轉(zhuǎn)換為細(xì)胞基因計數(shù)矩陣X，其中輸入Xig是細(xì)胞i表達(dá)的基因g的RNA轉(zhuǎn)錄本的數(shù)量。通常，細(xì)胞和基因被過濾以去除質(zhì)量差的細(xì)胞和最低表達(dá)的基因。然后，將數(shù)據(jù)歸一化以控制無意義的可變性來源，如測序深度、技術(shù)噪聲、庫大小和批處理效果。然后從歸一化數(shù)據(jù)中選擇高度可變基因(hvg)來識別感興趣的潛在基因并降低數(shù)據(jù)的維數(shù)。隨后，基因表達(dá)值被縮放到跨細(xì)胞的平均值為0，方差為1**。這種縮放主要是為了能夠應(yīng)用主成分分析(PCA)來進(jìn)一步降低維數(shù)，并提供有意義的嵌入來描述細(xì)胞之間的可變性來源。然后通過k近鄰(KNN)算法傳遞細(xì)胞的PCA嵌入，以便根據(jù)細(xì)胞的基因表達(dá)描述細(xì)胞之間的關(guān)系。KNN圖用于生成無向共享最近鄰(SNN)圖以供進(jìn)一步分析，最近鄰圖被傳遞到聚類算法中，將相似的單元分組在一起。圖(s)也用于進(jìn)一步的非線性降維，使用t-SNE或UMAP在二維中圖形化地描繪這些數(shù)據(jù)結(jié)構(gòu)。最后，通過差異表達(dá)(DE)分析鑒定cluster特異性marker基因，其中每個基因的表達(dá)在每個cluster與所有其他cluster之間進(jìn)行比較，并通過倍比變化和p值進(jìn)行量化。
Seurat是2015年用R語言編寫的，在生物信息學(xué)領(lǐng)域特別受歡迎;它是第一個全面的scRNA-seq分析平臺之一。Scanpy是2017年繼Seurat之后開發(fā)的一個基于python的工具，提供了一組類似的特性和功能。這兩個工具都有廣泛的運(yùn)用。Seurat和Scanpy之間的選擇通常歸結(jié)為用戶的編程偏好和他們的scRNA-seq數(shù)據(jù)分析項(xiàng)目的具體要求。
Seurat和Scanpy的輸入由一個基因計數(shù)矩陣組成，通常是cellranger生成的矩陣。
一個“標(biāo)準(zhǔn)的”scRNA-seq實(shí)驗(yàn)需要花費(fèi)數(shù)千美元，具體價格在很大程度上受數(shù)據(jù)大小的影響。雖然由于不同方法之間的差異，很難提供確切的成本，但據(jù)估計，一個典型的測序試劑盒的成本大約在數(shù)百到數(shù)千美元之間，測序成本每百萬次讀取5美元。對于高質(zhì)量數(shù)據(jù)，每個cell的必要讀取次數(shù)取決于實(shí)驗(yàn)的背景，但作為一個例子，cell Ranger通常建議其v3技術(shù)每個cell 20,000對reads，其v2技術(shù)每個cell 50,000對reads。樣品制備也有很大的成本，通常需要寶貴的患者樣本，或維持細(xì)胞或動物系數(shù)月至數(shù)年，為實(shí)驗(yàn)分析做準(zhǔn)備。一個標(biāo)準(zhǔn)的10x Genomics scRNA-seq實(shí)驗(yàn)序列數(shù)千萬到數(shù)十億的reads，根據(jù)環(huán)境的不同，推薦的細(xì)胞計數(shù)范圍為500-10k+。這些估計沒有考慮額外的成本，包括人工、實(shí)驗(yàn)設(shè)置和后續(xù)分析。
生物信息學(xué)數(shù)據(jù)分析中一個典型的隱含假設(shè)是，軟件和版本之間的選擇應(yīng)該對結(jié)果的解釋幾乎沒有影響。然而，在軟件或版本之間觀察到相當(dāng)大的可變性，即使在執(zhí)行其他類似或看似相同的分析時也是如此。這里研究的目的是量化標(biāo)準(zhǔn)scRNA-seq pipeline在軟件之間(即，Seurat與Scanpy)和同一軟件的多個版本之間(即，Seurat v5與v4, Scanpy v1.9與v1.4, Cell Ranger v7與v6)的變異性。
Seurat和Scanpy在默認(rèn)的scnaseq工作流中顯示出相當(dāng)大的差異
下圖顯示了使用PBMC 10k數(shù)據(jù)集與默認(rèn)設(shè)置比較Seurat v5.0.2和Scanpy v1.9.5的結(jié)果，展示了“標(biāo)準(zhǔn)”單細(xì)胞RNA-seq工作流的兩種實(shí)現(xiàn)之間的典型可變性。
在篩選UMIs、細(xì)胞最小基因數(shù)、基因最小細(xì)胞數(shù)和最大線粒體基因含量后，不同軟件之間的細(xì)胞或基因過濾沒有差異。此外，給定相同的矩陣作為輸入，Seurat和Scanpy也以相同的方式處理日志規(guī)范化，產(chǎn)生等效的輸出。然而，HVG選擇的默認(rèn)算法產(chǎn)生了差異，Jaccard index(兩組之間差異基因的交集/并集)為0.22。This difference could be resolved either by selecting the “seurat v3” flavor for Scanpy or the “mean.var.plot” algorithm for Seurat。
PCA分析開始觀察到更多的差異，使用默認(rèn)參數(shù)運(yùn)行時也會產(chǎn)生不同的結(jié)果。PCA圖顯示PC1-2空間中每個細(xì)胞的繪制位置存在明顯差異，盡管圖的大致形狀保持不變。Scree圖也顯示出差異，最明顯的是第一個PC解釋的方差比例相差0.1。PCA的變化都可以通過HVG設(shè)置標(biāo)準(zhǔn)化來解決，并相應(yīng)地調(diào)整PCA。
接下來，這些軟件在SNN圖的生成上有很大的不同。每個細(xì)胞的每個鄰域的內(nèi)容和大小都有很大差異。Seurat和Scanpy的每個細(xì)胞鄰域間的Jaccard指數(shù)中位數(shù)為0.11，median degree ratio(Seurat/Scanpy)量級為2.05。每個函數(shù)的degree ratio幾乎總是大于1，這表明Seurat在默認(rèn)情況下比Scanpy給出更多的高連接SNN圖?？紤]到SNN圖的點(diǎn)在所有degree ratio中相對均勻地分布在0和最大潛在Jaccard指數(shù)之間，似乎不是簡單的度差驅(qū)動低平均Jaccard指數(shù)。在生成SNN圖時對函數(shù)參數(shù)進(jìn)行對齊，degree ratio沒有質(zhì)的提高，但中位Jaccard指數(shù)有輕微的提高。
使用默認(rèn)設(shè)置的聚類也會導(dǎo)致輸出的差異，即使在調(diào)整函數(shù)參數(shù)和輸入SNN圖時，Seurat和Scanpy也證明了Louvain聚類的差異，但在Leiden算法的實(shí)現(xiàn)中是相同的。
UMAP圖在視覺上顯示了局部和鄰近c(diǎn)luster形狀的一些差異，即使在控制全局移動或旋轉(zhuǎn)的情況下。比較由這些UMAP數(shù)據(jù)構(gòu)建的KNN圖的鄰域相似性，發(fā)現(xiàn)鄰域重疊較差，隨著函數(shù)參數(shù)和先前輸入之間的相似性對齊，鄰域重疊會適度改善。對這些由UMAP導(dǎo)出的KNN圖進(jìn)行Leiden聚類和隨后的UMAP繪圖，發(fā)現(xiàn)軟件之間的UMAP繪圖的一般特征保持不變，但仍然存在一些相當(dāng)大的不可調(diào)和的差異。
在DE分析中，Seurat和Scanpy的顯著性差異基因的Jaccard指數(shù)重合為0.62(即，所有聚類的基因總數(shù)調(diào)整后p值< 0.05)，但Seurat的顯著性marker基因比Scanpy多約50%。顯著marker基因的差異是軟件間默認(rèn)設(shè)置的一些差異的結(jié)果。首先，每個軟件分別實(shí)現(xiàn)Wilcoxon功能，Seurat需要tie校正，而Scanpy默認(rèn)情況下省略tie校正。此外，每個軟件在默認(rèn)情況下調(diào)整p值不同- Seurat與Bonferroni multiple testing correction，而Scanpy與Benjamini-Hochberg多重測試校正。最后，Seurat在默認(rèn)情況下，在執(zhí)行Wilcoxon秩和檢驗(yàn)之前，通過p值、每組擁有該基因的細(xì)胞百分比和對數(shù)倍變化(logFC)過濾marker;Scanpy在不調(diào)用其他函數(shù)的情況下不會執(zhí)行這種類型的過濾。將Scanpy的過濾參數(shù)和聚類設(shè)置為與Seurat相同(filter：tie-correction, Bonferroni correction)進(jìn)行DE分析，將marker基因重疊的Jaccard指數(shù)提高到0.73，提供相同的聚類分配進(jìn)一步將Jaccard指數(shù)提高到0.99。其余1%的基因由于logFC計算的差異而存在差異。將方法設(shè)置為類似于Scanpy(沒有過濾，Benjamini-Hochberg)，使Jaccard指數(shù)惡化到0.38，因?yàn)闊o法去掉Seurat中的tie校正。
當(dāng)對齊cluster信息時，可以執(zhí)行進(jìn)一步的DE分析，比較每個cluster中每個基因的表達(dá)水平差異。除了比較所有聚類中顯著marker基因的外，還可以比較marker之間的相似性(即DE分析后每個聚類的基因)。Scanpy缺乏過濾意味著Scanpy包括所有cluster中的所有基因，即使該基因在該cluster中表達(dá)最少或無差異;而Seurat，默認(rèn)情況下，根據(jù)logFC, p值和reference中表達(dá)基因的細(xì)胞數(shù)量，每個cluster只包含很小比例的基因。對Scanpy應(yīng)用類似的閾值處理大大減少了這個問題，由于Scanpy缺少過濾，它將Jaccard指數(shù)從0.22提高到0.92。然而，這仍然不能完全調(diào)整差異。
Seurat和Scanpy計算logFC的方式也不同。比較各組間相似基因的一致性相關(guān)系數(shù)(CCC)為0.98,PCA擬合線斜率為1，表明各組間具有較強(qiáng)的相關(guān)性。簡而言之，CCC衡量兩個變量在相關(guān)性和方差方面的一致性。然而，通過觀察logFC值的散點(diǎn)圖，可以發(fā)現(xiàn)大量值之間存在顯著差異。具體來說，在少數(shù)情況下(135185個marker中的4109個)，Scanpy預(yù)測cluster中某個基因的logFC接近±30，而Seurat預(yù)測的logFC接近0。
在調(diào)整后的p值方面，Seurat和Scanpy之間也存在差異。對于默認(rèn)的函數(shù)參數(shù)，Seurat預(yù)測的p值要么小于或類似于Scanpy，但不會大得多。大多數(shù)p值接近最大值1，但存在很大程度的變異性。相當(dāng)多的p值遠(yuǎn)離y=x線，包括Seurat低于1e-50但Scanpy接近1的差異基因。20%的差異基因在軟件之間的p值在p=0.05閾值上翻轉(zhuǎn)，并且在兩個方向上翻轉(zhuǎn)相當(dāng)均勻(即僅在Seurat中顯著，或僅在Scanpy中顯著)。當(dāng)函數(shù)參數(shù)像Seurat一樣對齊時，幾乎所有調(diào)整后p值的差異都消失了。然而，由于在Seurat /presto的Wilcoxon秩和計算中缺乏切換校正的能力，這些差異無法與類似scanpy的函數(shù)參數(shù)相協(xié)調(diào)。
Seurat和Scanpy的比較表明，在某些情況下，程序結(jié)果是可以調(diào)和的，但并非總是如此。函數(shù)之間有三種可能的對齊方式:默認(rèn)對齊、匹配函數(shù)參數(shù)時對齊、不兼容對齊。下表顯示了每個函數(shù)在這些類中的分類。表1 (Seurat)和2 (Scanpy)詳細(xì)介紹了分析函數(shù)名的每個步驟、默認(rèn)參數(shù)、盡可能與其他軟件匹配所需的參數(shù)以及該軟件的唯一參數(shù)。
下采樣比較
考慮到軟件之間引入的可變性，一個自然的問題是如何對這些差異的大小進(jìn)行基準(zhǔn)測試。為此，在生成過濾UMI矩陣之前，模擬reads和細(xì)胞的下采樣，并比較了沿下采樣分?jǐn)?shù)梯度引入的差異與全尺寸數(shù)據(jù)。我們使用各自軟件的默認(rèn)函數(shù)參數(shù)執(zhí)行分析的每個步驟，并且在每個步驟之前不對齊輸入數(shù)據(jù)，除了DE分析的那些步驟，這些步驟需要在每個步驟之前對齊輸入，以便比較相同cluster中的差異基因統(tǒng)計(差異基因選擇，logFC和調(diào)整的p值)。對于分析的每一步，除了生成如圖1所示的所有圖外，我們還選擇了一個單一的數(shù)字指標(biāo)，該指標(biāo)將捕獲組間差異的程度，如下所示:
? Cell filtering: Jaccard index of cell sets
? Gene filtering: Jaccard index of gene sets
? HVG selection: Jaccard index of HVG sets
? PCA: Mean of difference in corresponding PC loadings PCs 1-3
? KNN/SNN: Median magnitude of log of SNN degree ratio
? Clustering: ARI
? UMAP: Median Jaccard index of UMAP-derived KNN neighborhoods across all cells
? Marker gene selection: Jaccard index of significant marker gene sets
? Marker selection: Jaccard index of marker sets
? logFC: CCC
? Adjusted p-value: Fraction of adjusted p-values that flipped across the p=0.05 threshold between conditions
對reads和細(xì)胞進(jìn)行降采樣的最不穩(wěn)健的步驟是基因選擇;然而，考慮到hvg和差異基因與下采樣部分的全尺寸數(shù)據(jù)集的相似性，似乎基因集的差異主要在于不太重要的基因。
版本控制對DE分析具有重要意義
除了軟件選擇(例如，Seurat vs. Scanpy)之外，軟件版本也可以在結(jié)果的解釋中發(fā)揮作用。將Seurat v5與v4進(jìn)行比較，在重要差異基因、marker和logFC估計值集方面存在相當(dāng)大的差異。logFC計算的差異源于不同版本間偽計數(shù)應(yīng)用程序的變化。Marker選擇的差異完全來自于logFC計算和過濾參數(shù)的差異。將Scanpy v1.9與較早的v1.4進(jìn)行比較還揭示了重要marker基因和marker list的巨大差異，這是由于刪除了不同版本之間的marker過濾。這些版本之間的logFC計算和調(diào)整后的p值沒有差異。比較使用默認(rèn)設(shè)置的Cell Ranger軟件v7和Cell Ranger v6生成的計數(shù)矩陣也揭示了所有DE指標(biāo)之間的差異?？鏑ell Ranger版本的分析顯示，pipeline的所有步驟都存在相當(dāng)大的差異。
這些命令之間的主要區(qū)別在于v7中默認(rèn)包含基因計數(shù)矩陣中的內(nèi)含子計數(shù)，而v6中默認(rèn)排除內(nèi)含子計數(shù)。這種區(qū)別對UMI過濾和每個細(xì)胞的基因小提琴圖有影響，與cell Ranger v6有更多的限制。
Random seeds
涉及隨機(jī)化影響的工作流步驟有近似knn搜索、Louvain/Leiden圖形聚類和UMAP。為了對軟件或數(shù)據(jù)大小之間的差異程度進(jìn)行基準(zhǔn)測試，我們使用相同的輸入數(shù)據(jù)和軟件選擇運(yùn)行這些步驟，只改變應(yīng)用的隨機(jī)種子。在相同的PCA輸入條件下，相同算法間SNN鄰域的Jaccard指數(shù)中位數(shù)和對數(shù)度比的變化(Annoy為0.85和0.05,umap-learn/PyNNDescent為1和0)遠(yuǎn)低于Seurat和Scanpy的0.27和1.61對數(shù)度比，表明軟件之間的差異不能僅僅用隨機(jī)性來解釋。隨機(jī)種子間Louvain聚類后的ARI為0.96，顯著高于Seurat和Scanpy中Louvain實(shí)現(xiàn)間的ARI(0.85)。UMAP衍生的隨機(jī)UMAP種子間KNN的Jaccard指數(shù)中值，Seurat為0.41,Scanpy為0.47，顯著高于相同輸入為0.21的Seurat和Scanpy之間的UMAP圖。然而，對于Seurat和Scanpy，在隨機(jī)UMAP種子中對相同數(shù)據(jù)進(jìn)行Leiden聚類后的ARI為0.64，與Seurat和Scanpy計算的觀察到的ARI相似，給定相同的PCA和SNN輸入，UMAP為0.69。這表明，盡管在Seurat或Scanpy中隨機(jī)種子之間生成的UMAP圖與軟件之間生成的UMAP圖具有更高的相似性，但Leiden算法不能完全捕獲這種相似性。
總結(jié)
Seurat和Scanpy在使用默認(rèn)設(shè)置執(zhí)行分析的方式上存在相當(dāng)大的差異，這些差異只能通過調(diào)整函數(shù)參數(shù)來部分調(diào)和。這些差異相當(dāng)于當(dāng)降采樣讀數(shù)小于5%或降采樣細(xì)胞小于20%時引入的可變性。計數(shù)矩陣生成和分析中涉及的軟件的版本控制也會對下游分析產(chǎn)生影響，特別是在沒有仔細(xì)考慮跨版本行為變化的情況下。為了在scRNA-seq分析中實(shí)現(xiàn)準(zhǔn)確性和可重復(fù)性，必須進(jìn)行一致的封裝選擇、深思熟慮的參數(shù)選擇和有意的版本控制。
