最近在進(jìn)行單細(xì)胞相關(guān)的數(shù)據(jù)的處理，由于樣本是不同的病理時期的，料想它們在測序的時候也是不同步進(jìn)行的茴厉，因此認(rèn)為批次效應(yīng)也肯定是存在的旅掂，但是之前一直沒有詳細(xì)分享過這方面的問題赏胚，在網(wǎng)上看了很多技術(shù)帖，有一點(diǎn)基礎(chǔ)的了解商虐，再進(jìn)一步的看看文獻(xiàn)原著觉阅，希望深入的了解單細(xì)胞的批次效應(yīng)，以及在不同的樣本數(shù)據(jù)情況下秘车，用何種批次處理效果好典勇，以及批次處理效果的標(biāo)準(zhǔn)如何評定。

在這里插入圖片描述

A benchmark of batch-effect correction methods for single-cell RNA sequencing data

用于普通的RNAseq數(shù)據(jù)整合時的批次矯正分析叮趴，如ComBat 和limma割笙，其也可以用于單細(xì)胞批次矯正處理。但是由于單細(xì)胞數(shù)據(jù)的drop out以及RNA 擴(kuò)增和捕獲過程中的失敗眯亦，因此需要新的批次處理的技術(shù)伤溉。

然而近期由Haghverdi提出的一個流行般码，且受歡迎的方法MNN(mutual nearest neighbors)可用于2個數(shù)據(jù)集之間的分析，在2個數(shù)據(jù)集中乱顾，配對的細(xì)胞列表用于計算翻譯向量板祝，進(jìn)而在共享空間比對。這種方法的優(yōu)點(diǎn)是數(shù)據(jù)進(jìn)行了歸一化處理走净，因而可用于下游的分析扔字。然而這種方法需要在高維基因表達(dá)空間計算臨近細(xì)胞列表，因而費(fèi)時費(fèi)力温技。需要進(jìn)一步的更進(jìn)革为。
基于上述的MNN 方法，研究者進(jìn)一步的開發(fā)了新的算法即FastMNN舵鳞，這種方法通過使用PCA分析在亞空間里使用MNN的方法震檩，進(jìn)而在運(yùn)行時間和精確性上進(jìn)一步得到提升。
其他的2種方法蜓堕，Scanorama and BBKNN也通過尋找亞空間的方法來運(yùn)用MNN抛虏。

對于我們最常見的單細(xì)胞分析的軟件Seurat來說，它在 進(jìn)行細(xì)胞的批次處理和數(shù)據(jù)整合的時候套才，進(jìn)行了更新?lián)Q代迂猴，最開始的版本2017年，Satjia lab開發(fā)了Multi canonical correlation analysis (CCA)的算法來對數(shù)據(jù)進(jìn)行降維背伴，然后再最相關(guān)的數(shù)據(jù)特征沸毁，然后進(jìn)行批次數(shù)據(jù)的比對。在Seurat3.0以上的版本中即Seurat Intergation傻寂，其對批次矯正的算法進(jìn)行了更新息尺，主要是2步，首先通過CCA 的算法對數(shù)據(jù)進(jìn)行降維得到亞空間疾掰，然后進(jìn)一步的通過 MNN的算法對CCA亞空間的細(xì)胞進(jìn)行計算搂誉，MNN的算法類似于一個anchors 的作用，將2個數(shù)據(jù)集的細(xì)胞進(jìn)行錨定静檬，如圖炭懊。

在這里插入圖片描述

除了上述提到的方法，近期還提到了其他的方法拂檩，比如說Harmony侮腹，其首先利用PCA進(jìn)行降維，在低維PCA的空間里广恢，通過迭代消除當(dāng)前的批次效應(yīng)凯旋，在每次迭代的過程中呀潭，他將不同批次的相似細(xì)胞聚類钉迷，最大化每個聚類批次的多樣性至非，然后計算應(yīng)用于每個細(xì)胞的批次矯正因子，這種方法能夠很快和準(zhǔn)確的進(jìn)行檢測糠聪。
LIGER是一種新的算法荒椭，主要是用于預(yù)知其他方法的缺點(diǎn)，它假設(shè)每個數(shù)據(jù)集之間的差異完全是由技術(shù)誤差導(dǎo)致的舰蟆，而非生物學(xué)效應(yīng)趣惠，那么這種假設(shè)將讓我們考慮完全除去這些技術(shù)誤差。LIGER通過使用非負(fù)矩陣分解首先獲取輸入數(shù)據(jù)的低維表示身害，該表示法由2部分組成味悄，一組特定批次因子和一組共有因子。此后執(zhí)行聚類塌鸯，然后使用共享因子鄰域圖搜索具有相似鄰域的細(xì)胞侍瑟，在已經(jīng)鑒定的細(xì)胞簇，然后將因子加載分位數(shù)歸一化去匹配選定的參考數(shù)據(jù)集丙猬，從而完成矯正涨颜。

機(jī)器學(xué)習(xí)領(lǐng)域深度神經(jīng)網(wǎng)絡(luò)近期也被用于單細(xì)胞的批次矯正過程。
在這篇文章中茧球，作者對以下14個算法MNN庭瑰，F(xiàn)ast MNN以及MultiCCASeurat 和Seurat3以及MMD-ResNet，Harmony和Scanorama 以及BBKNN和scGen Combat 和LIGER抢埋，limma 和scMerge 以及ZINB-WaVE等算法進(jìn)行了測試弹灭，在五種不同的場景下測試了十個數(shù)據(jù)集，五種場景分別是不同測序技術(shù)批次的已知細(xì)胞類型揪垄；含有未知細(xì)胞的批次鲤屡；多批次；超過百萬細(xì)胞的大數(shù)據(jù)福侈；用于基因差異表達(dá)的數(shù)據(jù)集模擬酒来。

結(jié)果

通過5項(xiàng)指標(biāo)來評估14中批次矯正的方法對十個數(shù)據(jù)集處理的效果

在這里插入圖片描述

本文中用于批次效應(yīng)矯正效果評估的五大指標(biāo)，

1.KBET ( k-Nearest neighbor batch-effect test)主要用于評估批次矯正后肪凛，不同的批次之間在局部的混合效果堰汉，本方法是基于奇異值降解的基礎(chǔ)上進(jìn)行降維。當(dāng)拒絕率低的時候伟墙，表明不同批次混合的結(jié)果好(其實(shí)他的拒絕率是一個零和假設(shè)翘鸭，也即所有的批次都能夠很好的混合，局部分布和整體分布的相似足夠大戳葵，而不被拒絕說批次混合不好)就乓，整個的拒絕率是0-1之間，當(dāng)拒絕率接近0的時候，就意味著不同批次的細(xì)胞批次效應(yīng)處理的好生蚁，細(xì)胞能夠很好的混合噩翠。在這種kbet 的評估中呢，主要是通過批次矯正后的產(chǎn)生PC成分或者是批次降維后的表達(dá)矩陣邦投，選取Top 的PC成分伤锚，由于K臨近值的數(shù)目對于對KBET 有很大的影響，因此作者運(yùn)行KBET通過K值預(yù)測的方法志衣，根據(jù)KBET文章中的方法屯援，作者選取K輸入值分別為樣本數(shù)量的5%，10%和15%和20%以及25%念脯，運(yùn)行KBET獲得所有的KBET 的拒絕值狞洋，然后取中位數(shù)，作為最終的KBET值绿店，用于對每種方法進(jìn)行評估徘铝。 最后通過Wilcoxon statistical test和BH 的矯正方法來確認(rèn)整合輸出的結(jié)果相比于其他的方法具有顯著的統(tǒng)計學(xué)意義。
2.LISCI(Local inverse Simpson’s index) 局部辛普森指數(shù) 惯吕，也是用來評估局部的水平惕它，來評估細(xì)胞的批次和細(xì)胞類型混合的結(jié)果。其包含iLISI(即整合局部辛普森指數(shù))和相比于KBET固定的K值废登，LISI主要通過具有固定的復(fù)雜度的局部距離分布選擇臨近的鄰居淹魄。臨近值的選擇然后用于計算辛普森指數(shù)多樣性，這是這個鄰域呈現(xiàn)出來的有效類型數(shù)堡距，指數(shù)通過用來計算批次標(biāo)簽甲锡，接近于期待批次數(shù)的分?jǐn)?shù)，表明批次混合效果好羽戒。其中iLISI分?jǐn)?shù)只計算所有批次中都出現(xiàn)的細(xì)胞缤沦，對于細(xì)胞類型的LISI，這些指數(shù)用于計算所有細(xì)胞類型的標(biāo)簽易稠，分?jǐn)?shù)接近1表明這種類型的細(xì)胞純度非常高缸废。在本文中作者為每個細(xì)胞計算了它們的cLISI值和iLISI值，然后進(jìn)一步的計算它們的中位數(shù)驶社，接著通過掃描最大和最小的分?jǐn)?shù)對中位數(shù)進(jìn)行scale企量，為了進(jìn)一步的聯(lián)合分析細(xì)胞的純度和批次的混合效應(yīng)，作者通過計算cLISI 和iLISI的調(diào)和平均數(shù)計算得到F值亡电，高的F值表明批次效應(yīng)矯正的越好届巩，Wilcoxon統(tǒng)計學(xué)檢驗(yàn)和Benjamini Hochberg 矯正用于對iLISI 和cLISI的結(jié)果進(jìn)行統(tǒng)計學(xué)顯著性檢測。

在這里插入圖片描述

3.ASW (Average silhouette width) 平均輪廓寬度
也是通過聯(lián)合細(xì)胞類型和批次標(biāo)簽來評估批次矯正的效果份乒，數(shù)據(jù)點(diǎn)的輪廓分?jǐn)?shù)恕汇，通過它距相鄰的Cluster的平均距離減去其與相同cluster 中的其他成員的平均距離腕唧，然后再除以2個值中的較大者，這種結(jié)果的值在-1到1之間瘾英，值越大表明其在當(dāng)前的Cluster 中矯正的越好枣接，值越低表明矯正的效果越差。同一個label 內(nèi)的所有數(shù)據(jù)點(diǎn)的平均值用來衡量細(xì)胞的類型和批次混合的效果方咆。在作者的這項(xiàng)研究中，作者首先將數(shù)據(jù)集隨機(jī)抽樣到原始數(shù)據(jù)集的80%蟀架，然后對通過向下取樣的方法獲得的數(shù)據(jù)集計算其前20個PC作為輸入數(shù)據(jù)來計算數(shù)據(jù)之間的距離以獲得ASW得分瓣赂，為了確保這個ASW得分的穩(wěn)定性，作者重復(fù)此過程20次片拍，為細(xì)胞類型混合和批次效應(yīng)混合獲取20個ASW煌集，然后將批次和細(xì)胞類型的中位數(shù)用于后續(xù)的進(jìn)一步的計算，作者逆轉(zhuǎn)了批次ASW值(越高越好)捌省，然后進(jìn)一步的將2個ASE值進(jìn)行歸一化到0和1之間苫纤，為了進(jìn)一步的聯(lián)合分析細(xì)胞類型純度和批次混合效應(yīng)，作者也進(jìn)一步的計算了批次和細(xì)胞類型SAW值的調(diào)和平均數(shù)纲缓，得到F1值卷拘，F(xiàn)1值越大(批次的ASW值越小，細(xì)胞類型的ASW值越高)表明批次矯正的效果越好祝高，Wilcoxon統(tǒng)計學(xué)檢驗(yàn)和Benjamini Hochberg 矯正用于對 ASW結(jié)果的檢驗(yàn)栗弟，來分析此種方法的ASW是否在統(tǒng)計學(xué)上優(yōu)于其他的方法。

在這里插入圖片描述

4.ARI (Adjusted rand index)能夠通過細(xì)胞類型和批次混合結(jié)果評估批次矯正的結(jié)果工闺，ARI用于評估兩個標(biāo)記列表里的匹配細(xì)胞的百分比乍赫，在作者的工作中，作者對數(shù)據(jù)集進(jìn)行二次取樣陆蟆，取原始數(shù)據(jù)的80%雷厂，然后對二次取樣后矯正過的數(shù)據(jù)獲取前20 的PC成分，通過使用K值函數(shù)進(jìn)行K值聚類叠殷，K是每個特有的細(xì)胞類型數(shù)目改鲫，為了使用ARI評估每個細(xì)胞類型的純度細(xì)胞類型標(biāo)記的結(jié)果與adjustedRandIndex 函數(shù)聚類的結(jié)果進(jìn)行比較，高的ARI的結(jié)果對應(yīng)的是高純度細(xì)胞類型林束，為了批次混合評估钩杰，只有出現(xiàn)在所有細(xì)胞簇的細(xì)胞用于后續(xù)的評估，批次標(biāo)記和K值聚類的結(jié)果標(biāo)記進(jìn)行比較诊县，低的ARI結(jié)果表明好的混合效果讲弄。和輪廓系數(shù)相似，為了產(chǎn)生穩(wěn)定的ARI結(jié)果依痊，作者重復(fù)了20次隨機(jī)二次取樣避除，為每個i細(xì)胞類型和批次獲得20個ARI分值怎披，隨后將ARI的中位數(shù)值進(jìn)行歸一化到0和1，通過調(diào)和平均數(shù)來獲取F1值瓶摆，綜合評估批次矯正的效果凉逛。高的F1值來源于一個低的批次混合的ARI和一個高的細(xì)胞類型的ARI，最后通過Wilcoxon統(tǒng)計學(xué)檢驗(yàn)和Benjamini Hochberg 矯正來檢驗(yàn)這種方法是否是統(tǒng)計學(xué)最顯著的批次矯正方法群井。

在這里插入圖片描述

5.DEG (Differential gene expression analysis)
為了實(shí)現(xiàn)差異基因表達(dá)分析状飞，作者在Seurat2中采用FindMarkers 函數(shù)通過似然比來檢測的那細(xì)胞基因表達(dá)情況，采用Bonferroni correction(BF) 的方法來獲取adjust 的p值书斜，且閾值設(shè)置為0.05 诬辈。對原始的數(shù)據(jù)進(jìn)行DEG 分析，它們包含批次效應(yīng)荐吉，也對通過其他批次矯正后的矩陣進(jìn)行DEG 分析焙糟。

DEG準(zhǔn)確性分析
主要是TP,TN,FP,FN 和precision 的計算

接下來作者通過對不同的場景進(jìn)行批次矯正處理
1.已知細(xì)胞類型，不同的技術(shù)
在這種情境下样屠，在包含不同批次處理的有相同細(xì)胞類型的數(shù)據(jù)集中進(jìn)行批次矯正檢驗(yàn)穿撮。
那么首先 作者對數(shù)據(jù)集2小鼠細(xì)胞圖譜進(jìn)行批次分析和細(xì)胞類型分析，并通過UMAP進(jìn)行可視化展示痪欲。結(jié)果如下

在這里插入圖片描述

可以明顯的看到悦穿，通過批次處理之后，數(shù)據(jù)集2中的數(shù)據(jù)能夠充分混合重疊的是經(jīng)過 Seurat 2, Seurat 3, Harmony, fastMNN, MNN Correct, scGen, Scanorama, scMerge,and LIGER successfully mixed the common cells這幾個軟件业踢，同樣的初看上去 這幾個軟件處理后的細(xì)胞類型聚類 也是這幾個軟件較好咧党，能夠相對較好的把每個細(xì)胞類型分開。ComBat 陨亡，limma傍衡，MMD-ResNet以及ZINB-WaVE以及BBKNN在相近的批次中能夠產(chǎn)生相似的細(xì)胞類型。

比較局部辛普森指數(shù)(iLISI)發(fā)現(xiàn)scMerge 是最好的方法用于批次混合负蠕，LIGER是第二好的方法蛙埂，這兩種方法的cLISI分?jǐn)?shù)也非常好(1-cLISI)的值大于0.96，對于KBET遮糖，Harmony批次混合效果也非常好绣的，接著是LIGER和scGen，p值都小于0.001欲账。使用ASW來評估的話屡江，Seurat3和Harmony 在批次混合和細(xì)胞聚類2方面都表現(xiàn)很好，其他的方法在批次混合效應(yīng)上也很好赛不，1-ASW的值>0.9惩嘉，在ARI用于批次混合分析的時候，所有方法都獲得了較好的ARI值踢故，Harmony的p值為0.001文黎，ARI的分?jǐn)?shù)是0.67惹苗，在多種方法里，Harmony 排名第一耸峭，接著是MNN和Seurat3桩蓉。

在這里插入圖片描述

對于第5個細(xì)胞集存在2對相似的細(xì)胞群，分別是CD4和CD8以及單核CD14和FCGR3A細(xì)胞劳闹，幾乎所有的方法都無法分出這兩個細(xì)胞簇院究，但是在Seurat2和3以及Harmony 和FastMNN 和MNN矯正過的數(shù)據(jù)里面，CD4和CD8亞群存在最小的細(xì)胞混合本涕，在這些方法處理過后业汰，這些亞群的分離是明顯可見的。而scGen和MMD-ResNet以及LIGER也能混合批次偏友，但是CD4和CD8細(xì)胞的混合也非常多蔬胯。Scanorama ,ZINB-WaVE以及scMerge 不僅混合了CD4和CD8細(xì)胞对供，而且也批次也混合了位他。

cLISI特點(diǎn)，發(fā)現(xiàn)細(xì)胞類型純度高達(dá)0.98产场，但是cLISI這一特點(diǎn)只用于分析局部的細(xì)胞純度鹅髓，特異細(xì)胞類型邊緣的混合效果不好。LIGER在批次混合上位于第一京景，而Seurat2和3位于其后窿冯，這與KBET計算的結(jié)果一致。除了ASW特點(diǎn)外确徙，所有的批次混合值都大于0.95醒串，而Harmony 和Seurat3都排名于top,且p值為0.183，接著是MNN鄙皇。和ARI一樣芜赌，就細(xì)胞類型的純度而言，Harmony 排名第一伴逸，接著是fastMNN和Seurat3和MNN缠沈，折4種方式的ARI的批次分?jǐn)?shù)都大于0.97，使用rank sum .Harmony和Seurat3和LIGER排名前三错蝴。

在這里插入圖片描述

對于所有的數(shù)據(jù)集來說洲愤，Harmony 都是排名第一，其次作者推薦的幾個方法是Seurat3和LIGER和FastMNN顷锰。

場景二
不確定的細(xì)胞類型

在這里插入圖片描述

多數(shù)批次矯正算法悼吱，都要求至少有一個細(xì)胞群在2個數(shù)據(jù)集間共享慎框，從而去指導(dǎo)數(shù)據(jù)的比對，其中MNN和FastMNN和Seurat3和Scanorama都通過尋找MNN來尋找2個數(shù)據(jù)共有的細(xì)胞類群后添，而當(dāng)2個數(shù)據(jù)間沒有共有的細(xì)胞亞群的時候MNN的錯誤匹配以及不正確的比對將會發(fā)生笨枯，尤其是2個數(shù)據(jù)集中的細(xì)胞亞群高度相似的時候，這也是為什么批次矯正的方法都將CD1C和CD141細(xì)胞整合到了一起遇西。

就KBET分?jǐn)?shù)來說馅精，LIGER 和Seurat2對數(shù)據(jù)集的批次整合效果都很好，LIGER和seurat2 都能夠很好的進(jìn)行批次整合粱檀，就iLISI 檢測結(jié)果來說洲敢，LIGER和Seurat 2都能夠得到很高的分值。

數(shù)據(jù)6只包含2種細(xì)胞類型茄蚯，而2/3的細(xì)胞都只含有一種細(xì)胞類型压彭，而t-SNE和UMAP表明scGen，scMerge 以及BBKNN都能夠產(chǎn)生2個大型的特異的細(xì)胞類群293T和Jurkat,Harmoy 也能夠很好的進(jìn)行批次混合渗常，能夠?qū)深惣?xì)胞分為2組壮不，LIGER也能夠產(chǎn)生批次混合的細(xì)胞簇，但是在一些細(xì)胞類型上存在混合皱碘，但是在Seurat2的輸出結(jié)果里询一，2種細(xì)胞類型的點(diǎn)靠得很近，難以分開癌椿。從LISI特點(diǎn)表明Harmony,scMerge 以及scGene都是最好的方法健蕊，從批次混合何細(xì)胞類型純度來說

在這里插入圖片描述

293T和Jurkat細(xì)胞的混合和聚類情況

在這里插入圖片描述

場景3多批次

在這里插入圖片描述

場景4大數(shù)據(jù)
BBKNN,ComBat, Harmony, LIGER, limma, MMD-ResNet, Scanorama,
scGen, Seurat 3, and ZINB-WaVE這幾個算法都能夠跑全庫，而其他的幾種算法不能踢俄，其中FastMNN 和scMerge以及Seurat2等算法缩功，需要的內(nèi)存高達(dá)2。27TB褪贵，而MNN 算法進(jìn)行批次矯正的時長約要48h以上掂之。
以上的那些算法中,LIGER能夠很好的維持細(xì)胞類型的分離，從而實(shí)現(xiàn)批次混合脆丁，Seurat 3,Harmony, ZINB-WaVE, scGen, and MMD-ResNet效果相對較差世舰，而Limma 和Combat 以及Scanorama，BBKNN這幾種算法實(shí)現(xiàn)批次混合的效果很差槽卫，幾乎沒有批次的混合跟压。總之歼培，scGen是最好的方法震蒋，批次混合的p值小于0.001茸塞，并于LIGER一起細(xì)胞類型純度達(dá)到0.34。

場景5數(shù)據(jù)的模擬
批次矯正分析查剖，主要是為了后續(xù)下游的偽時序以及基因差異分析钾虐。在此方法張作者通過綜合分析8種方法批次矯正后的基因差異分析。作者選用具有批次
效應(yīng)樣本的全部基因或者高變異基因作為批次矯正的輸入數(shù)據(jù)笋庄，批次校正后效扫，進(jìn)行DEG分析，最后通過TP,TN以及FP,FN來計算DEG 的準(zhǔn)確率直砂。那么作者通過計算發(fā)現(xiàn)菌仁，在使用HVG作為批次矯正的關(guān)鍵基因時，進(jìn)行批次矯正后能夠得到更真實(shí)的DEG静暂，具體的見圖20C济丘。

在這里插入圖片描述

那么如何解釋HVG對于批次矯正后，能夠更好的進(jìn)行后續(xù)差異分析洽蛀，主要是因?yàn)镠VG能夠更好的保留細(xì)胞類型自身的特點(diǎn)摹迷，因此使用選擇的HVG進(jìn)行批次矯正能夠提供更好的精確性，但是使用HVG 限制了后續(xù)的分析辱士，因?yàn)榛驍?shù)目有限泪掀。因此有必要通過全部的基因進(jìn)行批次矯正處理听绳，然后在得出結(jié)論的時候仔細(xì)小心以防出現(xiàn)假陽性或者假陰性FP或者FN颂碘。
通過場景5作者推薦用于后續(xù)的DEG分析的方法主要包括ComBat,MNN以及ZINB-WaVE以及scMerge。

討論
1.批次矯正的方法
2.批次整合評估方法包括KBET iLISI椅挣，ARI和ASW头岔。對于前2者來說，兩者的一致性非常高主要是二者都是基于鄰域的鼠证，而對于ARI的準(zhǔn)確性來說峡竣，主要依賴于聚類算法和簇的數(shù)目
3.運(yùn)行時間和內(nèi)存需求 推薦最佳的方法還是Harmony,LIGER以及Seurat3，對于大數(shù)據(jù)集的處理量九，Harmony 處理時間低于1h

結(jié)論
通過14中批次矯正的方法和5種情景适掰，作者發(fā)現(xiàn)綜合考慮，在進(jìn)行單細(xì)胞的批次處理的時候荠列，LIGER类浪，Harmony 以及Seurat3都能夠很好的進(jìn)行批次處理。Harmony 無論是對常見的細(xì)胞類型以及不同的技術(shù)和運(yùn)行時間上都是非常不錯的方法肌似，同樣LIGER也是费就，尤其在未知細(xì)胞類型的時候，對于LIGER來說主要的缺陷是運(yùn)行時間長川队，Seurat3也能處理大數(shù)據(jù)集力细，但是其運(yùn)行時間筆LIGER還要長20%-50%睬澡，然而對于若是進(jìn)行下游的DEG分析，作者推薦scMerge進(jìn)行批次矯正眠蚂。

全文用到的批次矯正方法

在這里插入圖片描述

生活很好煞聪，等你超越