hello，大家好，抓住這一周最后的小尾巴朴摊，給大家匯總一下我們單細(xì)胞數(shù)據(jù)尋找差異基因的22種方法，之前我分享的一篇文章10X單細(xì)胞（10X空間轉(zhuǎn)錄組）基因集打分方法匯總給大家匯總了各種基因集打分方法的匯總此虑，好了甚纲，彌補(bǔ)差異基因分析的這一塊拼圖。

參考的文章在Confronting false discoveries in single-cell differential expression朦前，是一篇方法論總結(jié)的文章介杆，發(fā)表在NC（Nature communications）鹃操，IF14分（上漲的還挺快）,好了，我們研究一下各個(gè)方法的優(yōu)缺點(diǎn)春哨。

圖片.png

先來總結(jié)一下結(jié)論荆隘。pseudobulk找差異基因的方法優(yōu)于單細(xì)胞直接找差異的方法，很多單細(xì)胞找差異基因的方法對高表達(dá)基因具有偏好性悲靴，主要是因?yàn)榘衙總€(gè)細(xì)胞當(dāng)成生物學(xué)復(fù)制的錯(cuò)誤判斷臭胜，從而造成了很多假陽性莫其，錯(cuò)誤的差異分析結(jié)果癞尚。值得我們注意。

單細(xì)胞轉(zhuǎn)錄組學(xué)中的差異表達(dá)分析能夠剖析細(xì)胞類型對疾病乱陡、創(chuàng)傷或?qū)嶒?yàn)操作等擾動的特異性反應(yīng)浇揩。 雖然許多統(tǒng)計(jì)方法可用于識別差異表達(dá)的基因，但區(qū)分這些方法及其性能的原則仍不清楚憨颠。在這里胳徽，研究表明這些方法的相對性能取決于它們解釋生物重復(fù)之間變異的能力。 忽略這種不可避免的變化的方法是有偏見的爽彤，并且容易出現(xiàn)錯(cuò)誤的發(fā)現(xiàn)养盗。 事實(shí)上，最廣泛使用的方法可以在沒有生物學(xué)差異的情況下發(fā)現(xiàn)數(shù)百個(gè)差異表達(dá)的基因（誤差很大）适篙。

介紹

RNA 種類的豐富程度反映了細(xì)胞和組織的過去往核、現(xiàn)在和未來狀態(tài)。 通過實(shí)現(xiàn) mRNA populations的完整量化嚷节，RNA 測序 (RNA-seq) 為了解生物樣品中活躍的分子過程提供了前所未有的途徑聂儒。 疾病、創(chuàng)傷和實(shí)驗(yàn)操作會擾亂這些過程硫痰，從而導(dǎo)致特定 mRNA 的表達(dá)發(fā)生變化衩婚。 從歷史上看，這些改變的 mRNA 是在非擾動與擾動組織中使用大量 RNA-seq 鑒定的效斑。 然而非春，生物組織由多種細(xì)胞類型組成，它們對擾動的反應(yīng)可能會有很大不同缓屠。 多細(xì)胞組織內(nèi) mRNA 豐度的變化被細(xì)胞類型的不同反應(yīng)和這些細(xì)胞類型的相對豐度的變化所混淆税娜。 因此，批量 RNA-seq 的分辨率不足以表征對生物擾動的多方面響應(yīng)藏研。

單細(xì)胞 RNA-seq (scRNA-seq) 能夠在單個(gè)細(xì)胞的分辨率下量化 RNA 豐度敬矩。單細(xì)胞技術(shù)的成熟現(xiàn)在能夠?qū)?fù)雜組織內(nèi)的細(xì)胞狀態(tài)進(jìn)行大規(guī)模比較，從而提供適當(dāng)?shù)姆直媛蕘砥饰黾?xì)胞類型對擾動的特異性反應(yīng)蠢挡。單細(xì)胞數(shù)據(jù)的稀疏性和異質(zhì)性最初鼓勵開發(fā)專門的統(tǒng)計(jì)方法來識別差異表達(dá)的 mRNA弧岳。用于差異表達(dá)分析的統(tǒng)計(jì)方法的激增促使研究人員query哪種方法產(chǎn)生的生物學(xué)結(jié)果最準(zhǔn)確凳忙。為了回答這個(gè)問題，研究人員轉(zhuǎn)向模擬禽炬，試圖創(chuàng)建一個(gè)可以對各種方法進(jìn)行基準(zhǔn)測試的基本事實(shí)涧卵。然而，模擬需要指定一個(gè)模型腹尖，從中生成差異表達(dá)的synthetic模式柳恐。

這些分歧強(qiáng)調(diào)了為單細(xì)胞數(shù)據(jù)中的差異表達(dá)建立健全的認(rèn)識論基礎(chǔ)的重要性。在這項(xiàng)工作中热幔，研究推斷開發(fā)這樣一個(gè)基礎(chǔ)需要在已知實(shí)驗(yàn)基礎(chǔ)事實(shí)的多個(gè)數(shù)據(jù)集中量化可用方法的性能乐设，并定義導(dǎo)致性能差異的原則。因此绎巨，首先建立了一個(gè)方法框架近尚，使我們能夠管理真實(shí)數(shù)據(jù)集資源。使用此資源场勤，我們對用于差異表達(dá)分析的各種可用方法進(jìn)行了明確的比較戈锻。發(fā)現(xiàn)這些方法的性能差異反映了某些方法無法解釋生物復(fù)制之間的內(nèi)在變異。對這一原理的理解使我們發(fā)現(xiàn)和媳，即使在沒有生物學(xué)差異的情況下格遭，最常用的方法也可以識別差異表達(dá)的基因。這些錯(cuò)誤的發(fā)現(xiàn)可能會誤導(dǎo)調(diào)查人員留瞳。然而拒迅，我們表明，使用解釋重復(fù)間變異的統(tǒng)計(jì)方法可以避免錯(cuò)誤發(fā)現(xiàn)撼港∑核總之，研究揭示了單細(xì)胞數(shù)據(jù)中有效差異表達(dá)分析的基礎(chǔ)原理帝牡，并集成了一個(gè)工具箱往毡，用于實(shí)現(xiàn)單細(xì)胞用戶的相關(guān)統(tǒng)計(jì)方法。

結(jié)果試驗(yàn)

A ground-truth resource to benchmark single-cell differential expression.

我們旨在比較可用的基于差異表達(dá) (DE) 分析的統(tǒng)計(jì)方法生成生物學(xué)準(zhǔn)確結(jié)果的能力靶溜。 我們推斷开瞭，在已知實(shí)驗(yàn)基本事實(shí)的真實(shí)數(shù)據(jù)集中執(zhí)行這種比較將忠實(shí)地反映這些方法的性能差異，同時(shí)避免模擬數(shù)據(jù)的缺點(diǎn)罩息。 我們假設(shè)嗤详，可以從匹配的批量和 scRNA-seq 中獲得與該基本事實(shí)最接近的近似值，該序列對相同的純化細(xì)胞群進(jìn)行瓷炮，暴露于相同的擾動葱色，并在相同的實(shí)驗(yàn)室中進(jìn)行測序。 對文獻(xiàn)的廣泛調(diào)查確定了總共 18 個(gè)符合這些標(biāo)準(zhǔn)的“黃金標(biāo)準(zhǔn)”數(shù)據(jù)集娘香。 該compendium使我們能夠在已知基本事實(shí)的實(shí)驗(yàn)環(huán)境中對 DE 方法進(jìn)行大規(guī)模比較苍狰。

圖片.png

• 注：A systematic benchmark of differential expression in single-cell transcriptomics. a Schematic overview of the eighteen ground-truth datasets analyzed in this study. b Statistical methods for DE analysis employed in 500 recent scRNA-seq papers. Grey bars represent DE analysis methods included in this study. “Other” includes methods used in two or fewer studies. Inset pie chart shows the total proportion of recent scRNA-seq papers that employed DE analysis methods included in this study. Source data are provided as a Source Data file. c Area under the concordance curve (AUCC) for fourteen DE methods in the eighteen ground-truth datasets shown in a. d Mean difference in the AUCC (Δ AUCC) between the fourteen DE methods shown in c. Asterisks indicate comparisons with a two-tailed t-test p-value less than 0.05. e AUCC of GO term enrichment, as evaluated using gene set enrichment analysis, in the eighteen ground-truth datasets shown in a. f Rank and statistical significance of the GO term GO:0043330 (“response to exogenous dsRNA”) in GSEA analyses of mouse bone marrow mononuclear cells stimulated with poly-I:C, a type of synthetic dsRNA, for four h, using the output of fourteen DE methods. Source data are provided as a Source Data file.

Pseudobulk methods outperform generic and specialized single-cell DE methods.

我們總共選擇了 14 種 DE 方法進(jìn)行比較办龄，這些方法代表了單細(xì)胞轉(zhuǎn)錄組學(xué)中使用最廣泛的統(tǒng)計(jì)方法（方法，“差異表達(dá)分析方法”）淋昭。 近 90% 的近期研究都使用了這些方法俐填。 我們根據(jù)批量與 scRNA-seq 數(shù)據(jù)集的 DE 結(jié)果之間的一致性評估了每種方法的相對性能。 為了量化這種一致性翔忽，我們計(jì)算了批量與 scRNA-seq 數(shù)據(jù)集結(jié)果之間的一致性曲線下面積 (AUCC)英融。

我們在 18 個(gè)黃金標(biāo)準(zhǔn)數(shù)據(jù)集的整個(gè)概要中比較了 14 種方法的性能。 該分析立即表明歇式，所有六種表現(xiàn)最佳的方法都具有共同的分析特性驶悟。 在應(yīng)用統(tǒng)計(jì)檢驗(yàn)之前，這些方法在生物復(fù)制品中聚集細(xì)胞贬丛，形成所謂的“假體”撩银。 相比之下给涕，比較單個(gè)細(xì)胞的方法表現(xiàn)不佳豺憔。 pseudobulk方法和單細(xì)胞方法之間的差異非常顯著，并且對用于量化一致性的方法具有魯棒性够庙。 此外恭应，與匹配的蛋白質(zhì)組學(xué)數(shù)據(jù)的比較表明，pseudobulk 方法還可以更準(zhǔn)確地預(yù)測蛋白質(zhì)豐度的變化

圖片.png

圖片.png

我們query DE 方法之間的差異是否也會影響轉(zhuǎn)錄組實(shí)驗(yàn)的功能解釋耘眨。 為此昼榛，我們比較了批量基因本體論 (GO) terms富集分析與 scRNA-seq DE。 我們發(fā)現(xiàn)pseudobulk方法再次更忠實(shí)地反映了真實(shí)情況剔难，如批量 RNA-seq 中捕獲的那樣胆屿。 例如，在比較用 poly(I:C)（一種合成雙鏈 RNA）刺激的小鼠吞噬細(xì)胞時(shí)偶宫，單細(xì)胞方法未能確定相關(guān)的 GO terms。

單細(xì)胞 DE 方法偏向于高度表達(dá)的基因。

pseudobulk方法的意想不到的優(yōu)越性迫使我們研究負(fù)責(zé)它們重述生物基本事實(shí)的能力的機(jī)制裁奇。 為了研究這些機(jī)制撼泛，我們制定并測試了幾個(gè)可能解釋這些性能差異的假設(shè)。

先前的研究表明吵冒，對于高度表達(dá)的基因纯命，關(guān)于 DE 的推斷通常更準(zhǔn)確。 單細(xì)胞中基因表達(dá)的測量本質(zhì)上是稀疏的痹栖。 通過在每個(gè)重復(fù)中聚集細(xì)胞亿汞，pseudobulk 方法顯著減少了數(shù)據(jù)中的零數(shù)量，特別是對于低表達(dá)的基因揪阿。 因此疗我，我們最初假設(shè)匙铡，pseudobulk 方法和單細(xì)胞方法之間的準(zhǔn)確性差異可以通過pseudobulk 方法在低表達(dá)基因中的優(yōu)越性能來解釋。

圖片.png

• Schematic illustration of the creation of ‘pseudobulks’ from single-cell data. Top, biological replicate from which each cell was obtained. Bottom, simulated gene expression matrix. Read counts for each gene are aggregated across cells of a given type within each biological replicate碍粥。

為了驗(yàn)證這一假設(shè)鳖眼，我們將基因分配到三個(gè)大小相同的bins中，包括低表達(dá)嚼摩、中度和高度表達(dá)的基因钦讳。 然后，我們重新計(jì)算了每個(gè) bin 內(nèi)批量和 scRNA-seq DE 之間的一致性枕面。 與我們的預(yù)測相反愿卒，我們觀察到低表達(dá)基因的pseudobulk方法和單細(xì)胞方法之間的差異很小。 相反潮秘，pseudobulk方法和單細(xì)胞方法之間最顯著的差異出現(xiàn)在高度表達(dá)的基因中琼开。

圖片.png

• 注： Mean AUCCs across eighteen ground-truth datasets after dividing the transcriptome into terciles of lowly, moderately, or highly expressed genes

這一意想不到的結(jié)果讓我們質(zhì)疑單細(xì)胞 DE 方法是否會對高表達(dá)基因產(chǎn)生系統(tǒng)性錯(cuò)誤。 為了探索這種可能性枕荞，我們仔細(xì)檢查了大量數(shù)據(jù)集柜候，以識別 scRNAseq 數(shù)據(jù)中的每種方法都錯(cuò)誤地稱為 DE 的基因。 我們發(fā)現(xiàn)由單細(xì)胞 DE 方法識別的假陽性比由pseudobulk方法識別的假陽性表達(dá)更高躏精。 相反渣刷，單細(xì)胞 DE 方法忽略的假陰性往往表達(dá)較低。 總之矗烛，這些發(fā)現(xiàn)暗示了單細(xì)胞方法將高度表達(dá)的基因識別為 DE 的系統(tǒng)趨勢辅柴，即使它們的表達(dá)保持不變。

圖片.png

為了通過實(shí)驗(yàn)驗(yàn)證這一結(jié)論瞭吃，我們分析了一個(gè)數(shù)據(jù)集碌嘀，其中將合成 mRNA populations摻入每個(gè)含有單個(gè)細(xì)胞的孔中。 因此歪架，這些單細(xì)胞文庫中的每一個(gè)都包含相同濃度的每種合成 mRNA股冗。 我們發(fā)現(xiàn)單細(xì)胞方法錯(cuò)誤地將許多豐富的spike-ins 識別為DE。 相比之下牡拇，pseudobulk方法避免了這種偏差魁瞪。

圖片.png

• d Spearman correlation between the mean expression of 80 ERCC spike-ins expressed in at least three cells and the –log10 p-value of differential expression assigned by each DE method. e Scatterplots of mean ERCC expression vs. –log10 p-value for exemplary single-cell and pseudobulk DE methods. Trend lines and shaded areas show local polynomial regression and the 95% confidence interval, respectively.

然后我們query這種偏差在單細(xì)胞轉(zhuǎn)錄組學(xué)中是否普遍存在。 我們匯編了 46 個(gè) scRNA-seq 數(shù)據(jù)集的compendium惠呼，其中包含不同的物種导俘、細(xì)胞類型、技術(shù)和生物擾動剔蹋。 我們發(fā)現(xiàn)單細(xì)胞 DE 方法在整個(gè)compendium中顯示出對高表達(dá)基因的系統(tǒng)偏好旅薄。

圖片.png

• 注：Rank and statistical significance of the GO term GO:0043330 (“response to exogenous dsRNA”) in GSEA analyses of mouse bone marrow mononuclear cells stimulated with poly-I:C, a type of synthetic dsRNA, for four h, using the output of fourteen DE methods. Source data are provided as a Source Data file.

Together, these experiments suggest that the inferior performance of single-cell methods can be attributed to their bias towards highly expressed genes.

圖片.png

圖片.png

DE analysis of single-cell data must account for biological replicates.

這些發(fā)現(xiàn)意味著pseudobulk方法擁有一個(gè)共同的分析特性，可以避免這種偏差。 我們進(jìn)行了一系列實(shí)驗(yàn)來確定這種機(jī)制少梁。

經(jīng)過多年的發(fā)展洛口，用于識別pseudobulk數(shù)據(jù)（即 edgeR、DESeq2 和 limma）中的 DE 基因的統(tǒng)計(jì)工具已經(jīng)得到改進(jìn)凯沪。 因此第焰，我們query這些方法是否結(jié)合了獨(dú)立于跨細(xì)胞聚合基因表達(dá)程序的固有優(yōu)勢。 為了測試這種可能性妨马，我們禁用了聚合過程并將pseudobulk方法應(yīng)用于單個(gè)細(xì)胞挺举。 引人注目的是，這個(gè)過程取消了pseudobulk方法的優(yōu)越性烘跺。 對高表達(dá)基因的偏見的出現(xiàn)與這種表現(xiàn)的下降并行湘纵。

圖片.png

• 注：a Schematic illustration of the experiment shown in b, in which the aggregation procedure was disabled and pseudobulk DE methods were applied to individual cells。b Left, AUCC of the original fourteen DE methods, plus six pseudobulk methods applied to individual cells, in the eighteen ground-truth datasets. Right, Spearman correlation between ERCC mean expression and –log₁₀ p-value assigned by six pseudobulk DE methods, before and after disabling the aggregation procedure.

這一結(jié)果提出了聚合過程本身直接導(dǎo)致pseudobulk方法優(yōu)越性的可能性滤淳。 為了評估這個(gè)概念梧喷，我們將聚合程序應(yīng)用于隨機(jī)的細(xì)胞組，這產(chǎn)生了一個(gè)由“pseudo-replicates”組成的pseudobulk矩陣脖咐。這個(gè)實(shí)驗(yàn)導(dǎo)致pseudobulk方法的性能出現(xiàn)類似的下降铺敌，并結(jié)合了 偏向于高表達(dá)基因.

圖片.png

圖片.png

圖片.png

我們試圖了解可以解釋這兩個(gè)實(shí)驗(yàn)中pseudobulk方法性能下降的共同因素。 我們認(rèn)識到這兩個(gè)實(shí)驗(yàn)都會丟失有關(guān)生物復(fù)制的信息文搂。 聚集隨機(jī)的細(xì)胞組形成pseudo-replicates适刀，或在單個(gè)細(xì)胞的比較中完全忽略復(fù)制秤朗，都引入了對高表達(dá)基因的偏見和相應(yīng)的性能損失煤蹭。

在相同的實(shí)驗(yàn)條件下，重復(fù)表現(xiàn)出基因表達(dá)的內(nèi)在差異取视，這反映了生物學(xué)和技術(shù)因素硝皂。 我們推斷，未能考慮到這些差異可能會導(dǎo)致方法將重復(fù)之間的固有可變性錯(cuò)誤地歸因于擾動的影響作谭。 為了研究這種潛在的機(jī)制稽物，我們比較了假體和假重復(fù)中每個(gè)基因表達(dá)的差異。 最初折欠，我們在用 poly-I:C12 刺激的骨髓單核細(xì)胞數(shù)據(jù)集中進(jìn)行了這種比較贝或。 我們發(fā)現(xiàn)重排復(fù)制品會導(dǎo)致基因表達(dá)方差的系統(tǒng)性降低，影響 98.2% 的基因锐秦。 我們接下來測試了這種方差的減少是否在我們的 46 個(gè)數(shù)據(jù)集概要中系統(tǒng)地發(fā)生咪奖。 每次比較都顯示出基因表達(dá)方差的相同減少

圖片.png

• 注：e Variance of gene expression in pseudobulks formed from biological replicates and pseudo-replicates in mouse bone marrow mononuclear cells stimulated with poly-I:C. Shuffling the replicate associated with each cell produced a systematic decrease in the variance of gene expression. Right, pie chart shows the proportion of genes with increased or decreased variance in pseudo-replicates, as compared to biological replicates. f Decreases in the variance of gene expression in pseudo-replicates as compared to biological replicates across 46 scRNA-seq datasets. Points show the mean variance in biological replicates; arrowheads show the mean variance in pseudo-replicates。

基因表達(dá)方差的減少導(dǎo)致統(tǒng)計(jì)測試將基因表達(dá)的微小變化歸因于擾動的影響酱床。 例如羊赵，在 poly-I:C 數(shù)據(jù)集中，未能考慮重復(fù)中 Txnrd3 的可變表達(dá)導(dǎo)致將該基因誤認(rèn)為差異表達(dá)扇谣。 此外昧捷，我們發(fā)現(xiàn)高表達(dá)基因在假重復(fù)中表現(xiàn)出最大的方差下降闲昭，從而解釋了單細(xì)胞方法對高表達(dá)基因的偏見

圖片.png

• g Left, expression of the gene c in biological replicates (points) and pseudo-replicates (arrowheads) from unstimulated cells and cells stimulated with poly-I:C, with the range of possible pseudo-replicate expression values shown as a density. Right, mean (horizontal line) and variance (shaded area) of Txnrd3 expression in biological replicates (left) and pseudo-replicates (right). P-values were calculated by edgeR-LRT.

這一系列實(shí)驗(yàn)共同揭示了pseudobulk方法出乎意料的優(yōu)越性的原理。 差異表達(dá)的統(tǒng)計(jì)方法必須考慮生物復(fù)制的內(nèi)在變異性靡挥，以在單細(xì)胞數(shù)據(jù)中產(chǎn)生生物學(xué)上準(zhǔn)確的結(jié)果序矩。 考慮到這種可變性，pseudobulk 方法可以正確識別由生物擾動引起的基因表達(dá)變化跋破。 相比之下贮泞，未能解釋生物復(fù)制會導(dǎo)致單細(xì)胞方法系統(tǒng)地低估基因表達(dá)的差異。 這種對方差的低估使單細(xì)胞方法偏向于高度表達(dá)的基因幔烛，從而損害了它們產(chǎn)生生物學(xué)準(zhǔn)確結(jié)果的能力啃擦。

False discoveries in single-cell DE.

我們意識到，如果不考慮生物學(xué)重復(fù)之間的差異可能會在存在真實(shí)生物學(xué)擾動的情況下產(chǎn)生錯(cuò)誤發(fā)現(xiàn)饿悬，那么在沒有任何生物學(xué)差異的情況下也可能出現(xiàn)錯(cuò)誤發(fā)現(xiàn)令蛉。為了測試這種可能性，我們在沒有任何組間差異的情況下模擬了重復(fù)之間具有不同程度異質(zhì)性的單細(xì)胞數(shù)據(jù)狡恬。我們將每個(gè)重復(fù)隨機(jī)分配到人工“對照”或“處理”組珠叔，并在兩種條件之間測試 DE。引人注目的是弟劲，單細(xì)胞方法在沒有任何擾動的情況下鑒定了數(shù)百個(gè) DE 基因祷安。此外，根據(jù)我們對單細(xì)胞 DE 方法失敗背后機(jī)制的理解兔乞，被錯(cuò)誤稱為 DE 的基因是那些在重復(fù)之間表達(dá)變化最大的基因汇鞭。 Pseudobulk 方法消除了 DE 基因的錯(cuò)誤檢測。然而庸追，創(chuàng)建pseudo-replicates導(dǎo)致虛假 DE 基因的重新出現(xiàn)霍骄，進(jìn)一步證實(shí)了準(zhǔn)確 DE 分析的要求。當(dāng)額外的重復(fù)被引入數(shù)據(jù)集時(shí)淡溯，錯(cuò)誤發(fā)現(xiàn)的數(shù)量減少了读整。相比之下，在模擬數(shù)據(jù)中引入額外的細(xì)胞只會加劇潛在的問題

這些發(fā)現(xiàn)迫使我們調(diào)查在真實(shí)的單細(xì)胞數(shù)據(jù)中是否會出現(xiàn)類似的錯(cuò)誤發(fā)現(xiàn)咱娶。 為了探索這種可能性米间，我們最初分析了暴露于干擾素 5 的人外周血單核細(xì)胞 (PBMC) 的數(shù)據(jù)集。 我們提取了未接觸過干擾素的對照樣品膘侮，并將它們隨機(jī)分成兩組屈糊。 然后我們進(jìn)行了 DE 分析。 未能解釋生物復(fù)制的內(nèi)在可變性在隨機(jī)分配的復(fù)制之間產(chǎn)生了數(shù)百個(gè) DE 基因

Unsettled by this appearance of false discoveries喻喳，我們query這一觀察是否反映了一個(gè)普遍的陷阱另玖。 為了全面解決這個(gè)問題，我們確定了總共 14 個(gè)數(shù)據(jù)集，其中在控制條件下至少包含 6 個(gè)重復(fù)谦去。 與之前的實(shí)驗(yàn)一樣慷丽，我們將這些未受干擾的樣本隨機(jī)分為合成對照組和處理組，然后在這兩組之間進(jìn)行 DE 分析鳄哭。 該系統(tǒng)分析證實(shí)要糊，與pseudobulk方法相比，單細(xì)胞方法產(chǎn)生了系統(tǒng)性過多的假陽性妆丘。 盡管完全沒有生物擾動锄俄，但產(chǎn)生的 DE 基因富集了數(shù)百個(gè)基因本體論 (GO) terms。 此外勺拣，我們再次證實(shí)奶赠，被錯(cuò)誤識別為 DE 的基因?qū)?yīng)于重復(fù)之間具有最高變異性的基因 。

總之药有，這些實(shí)驗(yàn)揭示了單細(xì)胞轉(zhuǎn)錄組學(xué)中 DE 分析的一個(gè)基本缺陷毅戈。 然而，我們的直覺是愤惰，這個(gè)陷阱可能會影響任何從每個(gè)生物復(fù)制中獲得許多觀察結(jié)果的技術(shù)苇经。 例如，我們預(yù)計(jì)空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)中也會出現(xiàn)錯(cuò)誤發(fā)現(xiàn)宦言。 為了測試這一預(yù)測扇单，我們分析了一個(gè)空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)集，該數(shù)據(jù)集對來自肌萎縮側(cè)索硬化 (ALS) 模型的脊髓進(jìn)行了分析奠旺。 我們將來自對照小鼠的數(shù)據(jù)隨機(jī)分為兩組蜘澜，并在脊髓的每個(gè)區(qū)域內(nèi)進(jìn)行 DE。 未能解釋生物復(fù)制之間變異性的統(tǒng)計(jì)方法在每個(gè)區(qū)域內(nèi)鑒定了數(shù)千個(gè) DE 基因凉倚。 相比之下兼都，pseudobulk方法消除了這些錯(cuò)誤發(fā)現(xiàn)。

這些實(shí)驗(yàn)表明稽寒，生物復(fù)制之間的變異性會混淆受生物擾動影響的基因的識別。在動物模型中趟章，可以最大限度地減少在重復(fù)之間產(chǎn)生這種變異性的許多因素杏糙，包括遺傳背景、環(huán)境蚓土、生物擾動的強(qiáng)度和時(shí)間以及樣品處理宏侍。相比之下，在涉及人類受試者的實(shí)驗(yàn)中蜀漆，這些變異來源本質(zhì)上更難以控制谅河。這種區(qū)別提出了一種可能性，即人體組織的單細(xì)胞研究在生物復(fù)制之間表現(xiàn)出更大的變異性，因此更容易受到錯(cuò)誤發(fā)現(xiàn)的影響绷耍。為了系統(tǒng)地評估這種可能性吐限，我們計(jì)算了 41 個(gè)人和小鼠 scRNA-seq 數(shù)據(jù)集中重復(fù)之間的變異性。與我們的假設(shè)一致褂始，我們檢測到人類數(shù)據(jù)集中重復(fù)之間的變異性明顯更大诸典。雖然我們表明考慮生物重復(fù)對于任何 DE 分析都是至關(guān)重要的，但這一結(jié)果強(qiáng)調(diào)了在人體組織的單細(xì)胞研究中解決這個(gè)問題的最重要的意義崎苗。

圖片.png

圖片.png

圖片.png

True and false discoveries in the injured mouse spinal cord

我們最終試圖證明 DE 分析可以在以前未探索的生物組織中產(chǎn)生真假發(fā)現(xiàn)的程度狐粱。 為此，我們描述了脊髓損傷 (SCI) 對損傷下方細(xì)胞中基因表達(dá)的影響胆数。 我們特別關(guān)注腰椎脊髓肌蜻，因?yàn)樵搮^(qū)域經(jīng)歷了多方面的變化，導(dǎo)致神經(jīng)元功能不可逆轉(zhuǎn)的退化必尼。

我們在遭受中胸脊髓嚴(yán)重挫傷的小鼠身上進(jìn)行了實(shí)驗(yàn)宋欺。 全身運(yùn)動學(xué)的多因素量化揭示了小鼠產(chǎn)生運(yùn)動能力的嚴(yán)重?fù)p害。 我們發(fā)現(xiàn)損傷引發(fā)了整個(gè)腰段新突觸的異常生長胰伍，并伴有異常節(jié)段反射的出現(xiàn)齿诞。 SCI 下方回路的這種混亂重組與痙攣和神經(jīng)元功能障礙有關(guān)

然后，我們收獲了慢性 SCI 和未受傷對照小鼠的腰脊髓骂租，并對這些組織進(jìn)行了單核 RNA-seq (snRNA-seq)祷杈。我們總共對 19,237 個(gè)細(xì)胞進(jìn)行了測序，這些細(xì)胞涵蓋了腰脊髓的所有主要細(xì)胞類型 .

我們最初的目標(biāo)是確定轉(zhuǎn)錄最受損傷干擾的細(xì)胞類型渗饮。 為了回答這個(gè)問題但汞，我們使用 Augur 進(jìn)行了細(xì)胞類型優(yōu)先排序。 這種無偏見的分析表明互站，內(nèi)皮細(xì)胞在損傷下方的脊髓中經(jīng)歷了最深刻的轉(zhuǎn)錄變化

這一意外發(fā)現(xiàn)促使我們研究這種優(yōu)先排序背后的特定轉(zhuǎn)錄變化私蕾，以及不同統(tǒng)計(jì)方法揭示這些變化的能力。 為此胡桃，我們使用代表性的單細(xì)胞和假體方法對受傷和未受傷的內(nèi)皮細(xì)胞進(jìn)行了 DE 分析踩叭。 我們選擇 Wilcoxon 秩和檢驗(yàn)作為單細(xì)胞方法，因?yàn)樵摍z驗(yàn)已成為單細(xì)胞轉(zhuǎn)錄組學(xué)領(lǐng)域中使用最廣泛的方法翠胰，而 edgeR-LRT 因其高性能水平而作為pseudobulk方法容贝。 這些方法確定了大部分不同的 DE 基因組，兩種方法之間只有四個(gè)基因重疊之景。 相反斤富，Wilcoxon 秩和檢驗(yàn)和 edgeR-LRT 分別指定另外 44 個(gè)和 12 個(gè)基因作為 DE

迄今為止，我們的結(jié)果表明锻狗，未能解釋重復(fù)之間的差異會導(dǎo)致單細(xì)胞 DE 方法產(chǎn)生錯(cuò)誤的發(fā)現(xiàn)满力。因此焕参，我們懷疑該數(shù)據(jù)集中通過 Wilcoxon 秩和檢驗(yàn)確定的一些其他基因可能代表假陽性。為了闡明這些基因在受傷脊髓中的真實(shí)表達(dá)油额，我們進(jìn)行了系統(tǒng)的體內(nèi)篩選叠纷。我們獲得了僅通過兩種方法之一鑒定的 19 個(gè)假定 DE 基因的 RNAscope 探針，并量化了這些基因在受傷和未受傷小鼠的內(nèi)皮細(xì)胞中的表達(dá)悔耘。 RNAscope 驗(yàn)證了被 edgeRLRT 稱為 DE 的六個(gè)基因中的五個(gè)讲岁。與此形成鮮明對比的是，Wilcoxon 秩和檢驗(yàn)稱為 DE 的十三個(gè)基因中只有三個(gè)可以得到證實(shí)（p < 0.05衬以，χ2 檢驗(yàn)缓艳；圖 5f-h）。一些經(jīng)過驗(yàn)證的 edgeR-LRT 基因看峻，包括 Slc7a11 和 Igfbp6阶淘，參與內(nèi)皮細(xì)胞對缺氧的反應(yīng)，支持在腰脊髓中建立慢性缺氧狀態(tài)互妓。與慢性缺氧的預(yù)期后果一致溪窒，我們檢測到損傷水平以下存在大量萎縮性血管

圖片.png

Together, these observations illustrate the potential for singlecell DE methods to produce false discoveries. Conversely, valid single-cell DE analysis that accounted for variation between biological replicates yielded reproducible conclusions that could be validated in vivo.

DE analysis with mixed models.

我們的實(shí)驗(yàn)表明，生物復(fù)制之間的差異決定了單細(xì)胞 DE 方法的性能冯勉。因此澈蚌，我們對線性混合模型令人不滿意的性能感到困惑。通過對生物復(fù)制內(nèi)部和之間的變化進(jìn)行明確建模灼狰，與pseudobulk方法相比宛瞄，混合模型應(yīng)該受益于更高的統(tǒng)計(jì)能力。為了澄清這種差異交胚，我們評估了八個(gè)額外的泊松或負(fù)二項(xiàng)式廣義線性混合模型 (GLMMs)份汗。在 25-50 個(gè)細(xì)胞的數(shù)據(jù)集中，GLMMs 可以在非常特定的參數(shù)組合下產(chǎn)生準(zhǔn)確的結(jié)果蝴簇。然而杯活，在包含 500 個(gè)或更多細(xì)胞的數(shù)據(jù)集中，它們的性能收斂于pseudobulk DE方法熬词。此外旁钧，擬合最佳GLMMs所需的計(jì)算資源是巨大的。即使在下采樣數(shù)據(jù)集中荡澎，單個(gè)細(xì)胞類型的DE分析平均需要13.5小時(shí)均践。相比之下，pseudobulk方法只需要在我們的 46 個(gè)數(shù)據(jù)集的概要中摩幔，每種細(xì)胞類型的分鐘數(shù)。這些觀察表明鞭铆，在實(shí)踐中或衡，pseudobulk 方法為單細(xì)胞 DE 分析提供了速度和準(zhǔn)確性之間的極好權(quán)衡焦影。

圖片.png

圖片.png

圖片.png

Discussion

需要在單細(xì)胞轉(zhuǎn)錄組學(xué)中進(jìn)行準(zhǔn)確的 DE 分析，以剖析對疾病封断、創(chuàng)傷和實(shí)驗(yàn)操作的多方面反應(yīng)背后的轉(zhuǎn)錄程序斯辰。 盡管統(tǒng)計(jì)方法對 DE 分析很重要，但決定其性能的原則仍然難以捉摸坡疼。 在這里彬呻，我們證明了有效 DE 分析的核心原則是統(tǒng)計(jì)方法能夠解釋生物復(fù)制的內(nèi)在變異性。 考慮到這種可變性柄瑰，決定了統(tǒng)計(jì)方法的生物學(xué)準(zhǔn)確性闸氮。 相反，在沒有任何生物學(xué)差異的情況下教沾，未能解釋生物復(fù)制變異性的方法可能會產(chǎn)生數(shù)百個(gè)錯(cuò)誤發(fā)現(xiàn)蒲跨。

研究人員研究單細(xì)胞以了解對生物擾動反應(yīng)的更一般原理。澄清這些原則需要統(tǒng)計(jì)推斷授翻，這些推斷超出構(gòu)成任何特定數(shù)據(jù)集的單個(gè)細(xì)胞或悲。我們的結(jié)果表明，通過在單個(gè)細(xì)胞水平上進(jìn)行統(tǒng)計(jì)推斷堪唐，單細(xì)胞 DE 方法將生物復(fù)制之間的變異性與生物擾動的影響混為一談巡语。重復(fù)之間存在差異是不可避免的，可以歸因于技術(shù)因素和內(nèi)在生物學(xué)差異淮菠。之前已經(jīng)認(rèn)識到男公，將重復(fù)之間的變異性與感興趣的生物學(xué)效應(yīng)混為一談的可能性可能導(dǎo)致虛假的發(fā)現(xiàn)18,34。然而兜材，這些研究幾乎完全依賴于合成數(shù)據(jù)理澎，并輔以一些說明性案例研究。因此曙寡，已發(fā)表的單細(xì)胞數(shù)據(jù)分析中錯(cuò)誤發(fā)現(xiàn)的普遍性以及這些錯(cuò)誤發(fā)現(xiàn)影響研究生物學(xué)結(jié)論的傾向仍不清楚糠爬。

在這里，我們表明錯(cuò)誤發(fā)現(xiàn)的出現(xiàn)是一種普遍現(xiàn)象举庶。利用具有實(shí)驗(yàn)基礎(chǔ)事實(shí)的 18 個(gè)單細(xì)胞數(shù)據(jù)集的集合执隧，我們證明使用不適當(dāng)?shù)慕y(tǒng)計(jì)方法會產(chǎn)生錯(cuò)誤的發(fā)現(xiàn)，從而損害單細(xì)胞實(shí)驗(yàn)的生物學(xué)解釋户侥。這些錯(cuò)誤的發(fā)現(xiàn)有可能浪費(fèi)時(shí)間镀琉、精力和財(cái)力來追求誤導(dǎo)性假設(shè)。例如蕊唐，我們通過對受傷小鼠脊髓的系統(tǒng)體內(nèi)篩選表明屋摔，通過最常用的統(tǒng)計(jì)方法鑒定的大多數(shù) DE 基因都是錯(cuò)誤的發(fā)現(xiàn)。此外替梨，我們闡明了未能解釋生物學(xué)和技術(shù)變異性使某些基因不成比例地被錯(cuò)誤地識別為 DE 的機(jī)制的進(jìn)展钓试。我們在額外的 46 個(gè)單細(xì)胞 RNA-seq 研究的多方面數(shù)據(jù)集中證明了這些機(jī)制的普遍性装黑。了解這些機(jī)制使我們發(fā)現(xiàn)，同樣的基本問題會影響其他高維分析弓熏，包括空間轉(zhuǎn)錄組學(xué)恋谭，并且最有可能出現(xiàn)在人體組織的研究中，這表明生物復(fù)制水平的推斷對于理解細(xì)胞和人類疾病的分子基礎(chǔ)挽鞠。

我們的結(jié)果表明疚颊，單細(xì)胞 DE 方法可能會產(chǎn)生錯(cuò)誤的發(fā)現(xiàn)。這種理解揭示了該領(lǐng)域的巨大風(fēng)險(xiǎn)信认。我們的研究結(jié)果表明材义，許多已發(fā)表的研究結(jié)果可能是錯(cuò)誤的。此外狮杨，如果不解決母截，可能會分配大量的研究資金來跟進(jìn)這些錯(cuò)誤的發(fā)現(xiàn)，從而損害科學(xué)橄教。然而清寇，使用解釋重復(fù)之間變異性的 DE 方法可以直接糾正這種相關(guān)的可能性。其中护蝶，我們發(fā)現(xiàn)pseudobulk方法在轉(zhuǎn)錄組华烟、蛋白質(zhì)組和功能解釋的水平上實(shí)現(xiàn)了對實(shí)驗(yàn)基本事實(shí)的最高保真度。因此持灰，我們認(rèn)為迫切需要對用于單細(xì)胞數(shù)據(jù) DE 分析的統(tǒng)計(jì)方法進(jìn)行范式轉(zhuǎn)變盔夜。我們的觀察強(qiáng)調(diào)了這種轉(zhuǎn)變的必要性，即過去兩年發(fā)表的大多數(shù)研究都使用了不適當(dāng)?shù)慕y(tǒng)計(jì)方法進(jìn)行 DE 分析堤魁。此外喂链，該領(lǐng)域使用最廣泛的分析包目前默認(rèn)采用易于錯(cuò)誤發(fā)現(xiàn)的 DE 方法。多條件數(shù)據(jù)集的日益流行強(qiáng)調(diào)了采用適當(dāng)?shù)慕y(tǒng)計(jì)方法來防止錯(cuò)誤發(fā)現(xiàn)擴(kuò)散的重要性妥泉。為了促進(jìn)這種轉(zhuǎn)變椭微，我們在 R 包中實(shí)現(xiàn)了這里測試的所有方法

關(guān)注一下方法

圖片.png

看看示例代碼

Libra is an R package to perform differential expression on single-cell data. Libra implements a total of 22 unique differential expression methods that can all be accessed from one function. These methods encompass traditional single-cell methods as well as methods accounting for biological replicate including pseudobulk and mixed model methods. The code for this package has been largely inspired by the Seurat and Muscat packages. Please see the documentation of these packages for further information.

System requirements

Libra relies on functions from the following R packages:

    dplyr (>= 0.8.0),
    purrr (>= 0.3.2),
    tibble (>= 2.1.3),
    magrittr (>= 1.5),
    tester (>= 0.1.7),
    Matrix (>= 1.2-14),
    pbmcapply (>= 1.5.0),
    lmtest (>= 0.9-37),
    tidyselect (>= 0.2.5),
    DESeq2 (>= 0.4.0),
    Seurat (>= 3.1.5),
    blme (>= 1.0-4),
    edgeR (>= 3.28.1),
    glmmTMB (>= 1.0.2.1),
    limma (>= 3.1-3),
    lme4 (>= 1.1-25),
    lmerTest (>= 3.1-3),
    matrixStats (>= 0.57.0),
    methods,
    stats,
    Rdpack (>= 0.7)

In addition, the Seurat, monocle3, or SingleCellExperiment packages must be installed for Libra to take Seurat, monocle, SingleCellExperiment objects as input, respectively. Methods that require additional packages may also require additional installs (e.g., MAST).

Libra has been tested with R version 3.6.0 and higher.

Installation

To install Libra, first install the devtools package, if it is not already installed:

> install.packages("devtools")

Libra additionally requires the following installations to perform differential expression testing:

> if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

> BiocManager::install(c("edgeR", "DESeq2", "limma"))

If Seurat is not installed this will be needed for single-cell methods.

> install.packages("Seurat")

Finally, install Libra from GitHub:

> devtools::install_github("neurorestore/Libra")

This should take no more than a few minutes.

Usage

The main function of Libra, run_de, takes as input a preprocessed features-by-cells (e.g., genes-by-cells for scRNA-seq) matrix, and a data frame containing metadata associated with each cell, minimally including the cell type annotations, replicates, and sample labels to be predicted. This means that in order to use Libra, you should have pre-processed your data (e.g., by read alignment and cell type assignment for scRNA-seq) across all experimental conditions.

Libra provides a universal interface to perform differential expression using 22 discrete methods. These methods are summarized as follows:

Single cell methods

• Wilcoxon Rank-Sum test
• Likelihood ratio test
• Student's t-test
• Negative binomial linear model
• Logistic regression
• MAST

Pseudobulk methods

• edgeR-LRT
• edgeR-QLF
• DESeq2-LRT
• DESeq2-Wald
• limma-trend
• limma-voom

Mixed model methods

• Linear mixed model
• Linear mixed model-LRT
• Negative binomial generalized linear mixed model
• Negative binomial generalized linear mixed model-LRT
• Negative binomial generalized linear mixed model with offset
• Negative binomial generalized linear mixed model with offset-LRT
• Poisson generalized linear mixed model
• Poisson generalized linear mixed model-LRT
• Poisson generalized linear mixed model with offset
• Poisson generalized linear mixed model with offset-LRT

By default Libra will use a pseudobulk approach, implementing the edgeR package with a likelihood ratio test (LRT) null hypothesis testing framework. Each of the 22 tests can be accessed through three key variables of the run_de function: de_family, de_method, and de_type. Their precise access arguments are summarized in the below table.

If batch effects are present in the data, these should be accounted for, e.g., using Seurat or Harmony, to avoid biasing differential expression by technical differences or batch effects.

To run Libra with default parameters on a genes-by-cells scRNA-seq matrix expr, and an accompanying data frame meta, with cell_type, replicate, and label columns containing cell types, replicates, and experimental conditions, respectively, use the run_de function:

> library(Libra)
> DE = run_de(expr, meta = meta)

If your columns have different names, you can specify these using the cell_type_col, replicate_col, and label_col arguments:

> DE = run_de(expr, meta = meta, cell_type_col = "cell.type", label_col = "condition")

If you would like to store the pseudobulk matrices in a variable, before running differential expression, you can do the following:

> matrices = to_pseudobulk(expr, meta = meta)

Libra can also run directly on a Seurat object. For a Seurat object sc, with the sc@meta.data data frame containing cell_type and label columns, simply do:

> DE = run_de(sc)

The same code can be used if sc is a monocle3 or SingleCellExperiment object instead.

Demonstration

To see Libra in action, load the toy single-cell RNA-seq dataset that is bundled with the Libra package:

> data("hagai_toy")

This dataset consists of 600 cells, distributed evenly between six replicates and two conditions.

The hagai_toy object is a Seurat object with columns named cell_type and label in the meta.data slot, meaning we can provide it directly as input to Libra:

> head(hagai_toy@meta.data)

                      orig.ident nCount_RNA nFeature_RNA                                       sample replicate label                                  cell_type
2-CGAACATGTATAATGG SeuratProject         18           11 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-ATTCTACAGTGGTAGC SeuratProject         23           12 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-GCATACACAAACTGTC SeuratProject          3            3 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-GAAGCAGAGATGCCAG SeuratProject         14            8 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-ATCACGAGTCTAGCGC SeuratProject          3            3 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-CTGCCTATCTGTCCGT SeuratProject         28           13 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes

We run run_de, and inspect the differential expression:

> DE = run_de(hagai_toy)
> head(DE)

We can also use a different statistical framework, for example DESeq2:

> DE = run_de(hagai_toy, de_family = 'pseudobulk', de_method = 'DESeq2', de_type = 'LRT')
> head(DE)

Alternatively, we can use a mixed-model approach, which by default will use a negative binomial model structure:

> DE = run_de(hagai_toy, de_family = 'mixedmodel')
> head(DE)

However, this can be adjusted using the de_method argument:

> DE = run_de(hagai_toy, de_family = 'mixedmodel', de_method = 'linear', de_type = 'LRT')
> head(DE)

Running this example on a MacBook should be instantaneous. However, analyzing >20 real single-cell RNA-seq datasets, we found Libra takes a median of ~5 minutes. In general, runtime scales close to linearly with the number of cell types and cells. If using mixed models, by default, Libra runs on four cores, with each gene analyzed on a different core. To change the number of cores, use the n_threads argument. For example, running Libra on eight threads:

> DE = run_de(hagai_toy, de_family = 'mixedmodel', de_method = 'linear', de_type = 'LRT', n_threads = 8)

... will run about twice as fast.

Calculating delta variance

We recently showed that statistical methods for differential expression must account for the intrinsic variability of biological replicates to generate biologically accurate results in single-cell data (Squair et al., 2021, Biorxiv; https://www.biorxiv.org/content/10.1101/2021.03.12.435024v1). Within the same experimental condition, replicates exhibit inherent differences in gene expression, which reflect both biological and technical factors. We reasoned that failing to account for these differences could lead methods to misattribute the inherent variability between replicates to the effect of the perturbation. To study this possibility, we compare the variance in the expression of each gene in pseudobulks and pseudo-replicates. We call this measure 'delta variance'. Users can use the calculation of delta variance to inform their differential expression results. For example, genes identified as differentially expressed by methods that do not account for biological replicate (i.e., 'singlecell' methods) that have a high delta variance should be treated with caution as they are likely to be false positives. Delta variance can be calculated as follows:

> DV = calculate_delta_variance(hagai_toy)

This function will return a list of vectors, one for each cell type, each of which contains the delta variance for the genes present in the input expression matrix.

生活很好，有你更好