Part 2 - 不用矯正的方法來處理批次效應
2020/10/24
在上一篇文章中，作者曾寫過過度處理批次效應的後果。這裡可能會有些爭議，他將對如何在不使用不知道具體原理的「批次效應矯正」工具的情況下處理批次效應給與一些建議。為什麼要這麼做？批次效應矯正工具真的那麼糟糕嗎？
並不是作者認為批次效應矯正方法不好。而是它們並不能百分百清楚的告訴用戶到底怎麼處理了數據，處理到什麼程度。因此，先自己處理數據，如果批次效應實在是很影響後續分析，再來使用所謂的批次矯正，因為並不是所有的批次效應都可以或者應該被矯正。
我們的公眾號生信技能樹之前也寫了很多關於批次效應的推文，例如
大家可以當做是背景知識讀一下！
去除已知的批次效應哈哈哈作者這句話很有意思，我的數據被詛咒了所以不好分析..其實是可能有一些批次效應作祟。
批次效應主要是想要去除實驗不同步驟之間的系統差異。需要保持實驗中的變量除外的環境都一致，這樣才會突出變量的差異。有些是生物學差異，但最主要的差異是由技術效應導致的。最常見的情況是，批次效應由許多大大小小的差異導致，每個差異都造成了總體不必要的差異的一部分。
也許是我們對不同的樣品使用了不同批次的試劑、用過不同化學成分的試劑盒、完全不同的單細胞測序平臺、一個樣本必須在冰上保存一夜、改變了組織分離方法、不同的樣品在不同的測序平臺。這些改變中的每一個都可能導致不同的結果。
雖然其中一些差異將導致很難甚至不可能預測和逆轉的轉錄變化，還有一些則可以輕易去除。因此，對有批次效應的數據，可以做的事就是先去除已知的批次效應。
這裡將重點介紹 10X 單細胞表達數據。
去除 swapped barcodes如圖所示，在新型號的 Illumina 測序機上，流動池中的連接和DNA擴增同時進行，各步驟之間不需要對流動池進行任何清洗。樣本標記仍保留在溶液中，並且可以使用具有不同標記的庫中的DNA分子作為模板擴展。錯誤標記的分子在流動池的納米孔之間的轉移導致錯誤標記的DNA分子的聚集和測序。
首先去除 swapped barcodes. 用這個R包 DropletUtils 來處理很簡單~ (作者給與了友情提示說這個包不好下載，結果嗖的一下就下完了..eee).
library(DropletUtils)
cleaned = swappedDrops('data/outs/molecule_info.h5')大多數情況下，這一步並不能起到多大作用。但曾經有過這樣的情況，在新型的以Nano-Well 為特點的Patterned Flow Cell Technology（PFCT）的測序儀上運行的數據得到了極大的改善。在較新的 10X kits 中引入雙重標籤應該可以避免這個問題。
去除環境RNAAmbient gene expression（環境基因表達）是指不是來自 barcode 細胞，而是來自其他溶解細胞的count，這些細胞的 mRNA 在文庫構建之前汙染了細胞懸液。這些增加的環境計數會扭曲下遊分析，如標記基因鑑定或其他差異表達檢測，尤其是當樣本之間的水平變化時。在基於液滴的 scRNA-seq 數據集中校正這些影響是可能的，由於大量的空液滴，可用於模擬環境RNA表達譜。SoupX使用這種方法直接糾正計數數據。
第二件事是去除環境RNA汙染。這一點非常重要，因為環境RNA看起來與實驗中所有細胞的平均值非常相似。下圖顯示了這一點，它將背景汙染中每個基因的平均表達與實驗中所有細胞的平均表達進行了比較。
這樣做的結果是，如果有兩個(或更多)具有不同混合細胞的數據，最終看起來會有所不同，因為背景環境不同。去除環境汙染可以使用作者寫的這個 SoupX ，但也有其他方法。環境RNA汙染會產生一種幾乎總是存在的批次效應，去除它可以明顯提高整合度。此圖顯示了將 SoupX 應用於某些數據之前和之後的批處理熵。熵越高，表明去除背景時混合效果越好。
在這篇論文 SoupX paper 裡對環境RNA汙染及其影響有更完整的解釋。
處理細胞周期細胞處於細胞周期的不同階段並不是嚴格意義上的批次效應。但是它的關聯性足夠強，並且可以對數據聚類的方式產生很大影響，所以這裡還是要討論一下。周期細胞可以使真正不相關的細胞看起來是相關的，因為它們都處於 S、G2 或 M 期。大多數人對作者的擔憂的回應都是這樣的：「沒關係，可以回歸的去除細胞周期的影響。」
下面的圖表可以清楚地看到，處於 S 期( PCNA 陽性)和 G2M 期( MKI67 期)的細胞與它們來自哪個簇都是不同的(這些細胞是 this paper 中的乳頭狀腎細胞癌細胞)。在這裡，使用 PCNA 作為 S 期標記物，MKI67 作為 G2M 期標記物，作者認為觀察標記物的表達比 Seurat 的 CellCycleScoring 函數更可靠。（這就是生物學背景知識了，在大家處理數據的過程中要儘可能多記憶一些重要的細節知識點！！！）
去除細胞周期基因
但是如果我們丟棄所有的周期基因呢？這應該會使它們和 G1 期細胞一樣吧？在Seurat我們可以很容易地做到這一點
phaseS = cnts['PCNA',]>0
phaseM = cnts['MKI67',]>0
altCnts = cnts[!rownames(cnts) %in% unlist(cc.genes.updated.2019),]
#Usual stuff...
altCnts@meta.data$Phase = ifelse(phaseS,'S',ifelse(phaseM,'G2M','G1'))
DimPlot(altCnts,group.by='Phase')顯然，處於周期的哪個階段的定義有點粗略，但概念相當清晰。去除主要的細胞周期基因並沒有起到多大作用。如果你能找到一份更完整的細胞周期相關基因列表，應該會更好的去除這些基因。但是，考慮到基因在不同環境中具有不同的功能，這種方法發揮作用的機率很小。（生物學知識確定有點雜亂無章）
回歸的去除細胞周期的影響
使用 Seurat vignette 來看看 regress cell cycle out 會怎樣：
reg = CreateSeuratObject(cnts)
reg = NormalizeData(reg)
reg = CellCycleScoring(reg,s.features=cc.genes.updated.2019$s.genes,g2m.features=cc.genes.updated.2019$g2m.genes)
reg = ScaleData(reg,vars.to.regress=c('S.Score','G2M.Score'))
reg = FindVariableFeatures(reg)
reg = RunPCA(reg)
reg = RunUMAP(reg,dims=seq(50))
DimPlot(reg,group.by='Phase')這看起來好多了，但顯然還是有區別的。也許這種差異是潛在的細胞狀態之間的某種真正的差異，這種狀態不是細胞周期，而是細胞本身固有的。但是，除非你真的缺乏細胞，或者有一些細胞類型只存在於S/G2M期，否則根本不需要這麼麻煩，只需去除周期細胞，這就不再是問題了。
去除無關緊要的基因
在某些情況下，在分析之前去除不需要的基因可以有效防止批次效應的產生。
例如，有一些已知的基因會因組織處理而發生變化(例如，熱休克基因)。在這種情況下，你很可能會懷疑和這些基因相關的一些簇是否是和基因相關，或者是實驗處理過程中的變化。那麼為什麼不乾脆把他們排除在外呢？要排除的基因通常是線粒體基因和熱休克基因。但這還取決於實驗設計。
顯然，這些基因中的信息有時仍然有用，這個想法只是為了防止它們聚類。可能最穩健的方法是將它們從輸入計數矩陣中排除，但將信息作為元數據保留，以便不時之需。例如,
mtGenes=grep('^MT-',rownames(cnts),value=TRUE)
redCnts = cnts[!rownames(cnts) %in% mtGenes,]
mtFrac = colSums(cnts[mtGenes,])/colSums(cnts)
red = CreateSeuratObject(redCnts)
red@meta.data$mtFrac = mtFrac這將消除MT基因對分析的任何影響，但是如果您希望使用作為元數據提供的信息。另一種選擇是將它們從用於計算PC和下遊分析的高度可變基因中排除。例如。
reg@assays$RNA@var.features = setdiff(reg@assays$RNA@var.features,mtGenes)其實處理到最後，會感覺我們重要的單細胞轉錄組其實並不是全局的基因表達量啦，有點類似於特定的gene panel，也就是說目前主流的10X單細胞儀器其實可以被BD平臺的單細胞取代？
管它幹啥~批次效應在 tSNE 或 UMAP 上看起來很難看，可能至少會有一個編輯抱怨它。
主要取決於想要在分析中實現什麼目標，以及通過改變數據來消除批次效應是否真的會進一步推動這一目標。
例如，單細胞分析的一個共同目標是在實驗中定義特定細胞類型特有的基因。從表面上看，將你的細胞類型分成多個簇似乎會阻礙這一努力。但是，所有這些實際上都會使注釋數據的任務變得更加困難。作者從來沒有遇到過不想合併兩個按算法選擇的簇的注釋的單細胞分析。
雖然沒有任何確鑿的證據，但作者認為如果提供的數據沒有通過批次效應處理，那麼標記查找方法將會做得和處理後的數據一樣好，甚至更好。當然，如果您將標記查找定義為一項差異表達任務 (which I don’t think you should)，那麼使用 edgeR 或 DESeq2 會得到更好的結果，這些方法需要輸入 raw counts。
結論有時，你的數據有很多技術差異，即使在完成上述所有操作之後，這些差異也不會去除。
或者，您的主要興趣可能是尋找細胞間表達的連續變化(例如，偽時間)。在這些情況下，別無選擇，只能使用自動批次效應矯正工具。
但希望您可以嘗試一下上面的方法去除，而不是默認運行批次效應矯正直接去除。
翻譯自：https://constantamateur.github.io/2020-10-24-scBatch2/更多教程見其博客：https://constantamateur.github.io/