生物信息學(xué)習(xí)的正確姿勢

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測序分析?(重磅綜述：三萬字長文讀懂單細(xì)胞RNA測序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

高通量數(shù)據(jù)中批次效應(yīng)的鑒定和處理（一）

高通量數(shù)據(jù)中批次效應(yīng)的鑒定和處理（二）

高通量數(shù)據(jù)中批次效應(yīng)的鑒定和處理（三）- 如何設(shè)計(jì)盡量避免批次影響

高通量數(shù)據(jù)中批次效應(yīng)的鑒定和處理（四）- 在差異基因鑒定過程中移除批次效應(yīng)

高通量數(shù)據(jù)中批次效應(yīng)的鑒定和處理（五）- 預(yù)測并校正可能存在的混雜因素

直接校正表達(dá)矩陣

處理批次因素最好的方式還是如前面所述將其整合到差異基因鑒定模型中，降低批次因素帶來的模型殘差的自由度。但一些下游分析，比如數(shù)據(jù)可視化，也需要直接移除效應(yīng)影響的數(shù)據(jù)來展示，這時(shí)可以使用ComBat或removeBatchEffect函數(shù)來處理。

輸入數(shù)據(jù)，標(biāo)準(zhǔn)化且log轉(zhuǎn)換后的數(shù)據(jù) ehbio.simpler.DESeq2.normalized.rlog.xls

id untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011 trt_N61311 trt_N052611 trt_N080611 trt_N061011 ENSG00000115414 18.02 18.62 17.83 18.45 17.95 18.54 18.15 18.50 ENSG00000011465 17.79 18.36 17.96 18.52 17.79 18.23 17.84 18.47 ENSG00000091986 17.15 17.74 16.41 17.59 17.27 17.79 16.88 17.76 ENSG00000103888 15.56 16.90 15.88 16.42 15.94 17.43 17.38 17.05

包含已知批次信息和預(yù)測的批次信息的樣本屬性文件 sampleFile2

Samp conditions individual sizeFactor SV1 SV2 SV3 untrt_N61311 untrt N61311 1.0211325 -0.10060313 -0.4943517 -0.31643389 untrt_N052611 untrt N052611 1.1803986 0.01827734 -0.1701068 0.58841464 untrt_N080611 untrt N080611 1.1796083 -0.42949195 0.3756338 -0.08929556 untrt_N061011 untrt N061011 0.9232642 0.53452392 0.2413738 -0.17649091 trt_N61311 trt N61311 0.8939275 -0.12535603 -0.4956603 -0.36550102 trt_N052611 trt N052611 0.6709229 0.03588273 -0.151201 0.5914179 trt_N080611 trt N080611 1.3967665 -0.46668403 0.4413431 -0.07016903 trt_N061011 trt N061011 0.9462307 0.53345114 0.2529692 -0.16194213

加載需要的包

# 若缺少YSX包，則安裝 # BiocManager::install("Tong-Chen/YSX", update=F) suppressMessages(library(DESeq2)) suppressMessages(library("RColorBrewer")) suppressMessages(library("gplots")) suppressMessages(library("amap")) suppressMessages(library("ggplot2")) suppressMessages(library("BiocParallel")) suppressMessages(library("YSX")) suppressMessages(library(sva)) suppressMessages(library(ggfortify)) suppressMessages(library(patchwork)) suppressMessages(library(ggbeeswarm)) suppressMessages(library(limma))

讀入標(biāo)準(zhǔn)化后的表達(dá)矩陣和樣品信息表

expr_File <- "ehbio.simpler.DESeq2.normalized.rlog.xls" expr_mat <- sp_readTable(expr_File, row.names=1)head(expr_mat)sampleFile <- "sampleFile2.txt" metadata <- sp_readTable(sampleFile, row.names=1) head(metadata)

使用ComBat校正

ComBat校正時(shí)考慮生物分組信息

biological_group = "conditions" batch = "individual"metadata[[biological_group]] <- factor(metadata[[biological_group]]) metadata[[batch]] <- factor(metadata[[batch]])# 模型中引入關(guān)注的生物變量和其它非批次變量，保留生物差異和非批次差異 modcombat = model.matrix(as.formula(paste('~', biological_group, sep=" ")), data= metadata)# ComBat需要的是matrix expr_mat_batch_correct <- ComBat(dat=as.matrix(expr_mat), batch=metadata[[batch]], mod=modcombat) expr_mat_batch_correct <- as.data.frame(expr_mat_batch_correct)expr_mat_batch_correct[1:3,1:4]untrt_N61311 untrt_N052611 untrt_N080611 untrt_N061011 ENSG00000109906 5.107739 5.201200 5.293888 5.092441 ENSG00000152583 7.002840 7.199577 7.281516 7.103559 ENSG00000224114 6.428427 6.433242 6.236464 6.235308

校正后的PCA結(jié)果顯示在PC1軸代表的差異變大了，PC2軸代表的差異變小了，不同來源的樣本在PC2軸的分布沒有規(guī)律了 (或者說成鏡像分布了)。

sp_pca(expr_mat_batch_correct[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

ComBat校正時(shí)不考慮分組信息，也可以獲得一個(gè)合理的結(jié)果，但是一部分組間差異被抹去了。

# ComBat需要的是matrix expr_mat_batch_correct2 <- ComBat(dat=as.matrix(expr_mat), batch=metadata[[batch]], mod=NULL) expr_mat_batch_correct2 <- as.data.frame(expr_mat_batch_correct2)sp_pca(expr_mat_batch_correct2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

關(guān)于運(yùn)行ComBat時(shí)是否應(yīng)該添加關(guān)注的生物分組信息，即mod變量，存在一些爭議。反對(duì)添加mod的人的擔(dān)心是這么處理后，是否會(huì)強(qiáng)化生物分組之間的差異。支持添加mod的人是擔(dān)心如果不添加mod那么去除批次時(shí)可能也會(huì)去除樣本組之間的差異，尤其是實(shí)驗(yàn)設(shè)計(jì)不合理時(shí)。

如果是非平衡實(shí)驗(yàn)，類似我們在實(shí)驗(yàn)設(shè)計(jì)部分時(shí)提到的方案2，則沒有辦法添加mod變量，程序會(huì)報(bào)出Design matrix is not full rank類似的錯(cuò)誤，這時(shí)是不能區(qū)分差異是來源于批次還是來源于樣本，強(qiáng)行移除批次時(shí)，也會(huì)移除一部分或者全部樣本分組帶來的差異。這個(gè)在第一篇帖子處有兩位朋友的留言討論可以參考。

ComBat只能處理批次信息為l離散型分組變量的數(shù)據(jù)，不能處理sva預(yù)測出的連續(xù)性混雜因素。

使用limma校正

如果批次信息有多個(gè)或者不是分組變量而是類似SVA預(yù)測出的數(shù)值混雜因素，則需使用limma的removeBatchEffect (這里使用的是SVA預(yù)測出的全部3個(gè)混雜因素進(jìn)行的校正。)。

樣品在PC1和PC2組成的空間的分布與ComBat結(jié)果類似，只是PC1能解釋的差異略小一些。

SV = metadata[,c("SV1","SV2","SV3")] expr_mat_batch_correct_limma1 <- removeBatchEffect(expr_mat, covariates = SV, design=modcombat) sp_pca(expr_mat_batch_correct_limma1[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

removeBatchEffect運(yùn)行時(shí)也可以不提供生物分組信息。

expr_mat_batch_correct_limma1 <- removeBatchEffect(expr_mat, covariates = SV) sp_pca(expr_mat_batch_correct_limma1[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

removeBatchEffect也可以跟ComBat一樣，對(duì)給定的已知一個(gè)或多個(gè)定性批次信息進(jìn)行校正。

expr_mat_batch_correct_limma2 <- removeBatchEffect(expr_mat, batch=metadata[[batch]], design=modcombat) sp_pca(expr_mat_batch_correct_limma2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

不指定目標(biāo)分組變量，結(jié)果也不受影響。

expr_mat_batch_correct_limma2 <- removeBatchEffect(expr_mat, batch=metadata[[batch]]) sp_pca(expr_mat_batch_correct_limma2[1:5000,], metadata,color_variable="conditions", shape_variable = "individual") +aes(size=1) + guides(size = F)

同時(shí)考慮批次、混雜因素和生物分組信息進(jìn)行校正，校正后差異就全部集中在生物分組信息水平 (PC1)上了 (PC1 variance=100)，應(yīng)該是過擬合了，每組樣本的基因表達(dá)都一致了。

expr_mat_batch_correct_limma3 <- removeBatchEffect(expr_mat, batch=metadata[[batch]], covariates = SV, design=modcombat) sp_pca(expr_mat_batch_correct_limma3[1:5000,], metadata,color_variable="conditions", shape_variable = "individual", coord_fixed_ratio=0) +aes(size=1) + guides(size = F)

封面來源于：https://www.pexels.com/zh-cn/photo/264537/

往期精品(點(diǎn)擊圖片直達(dá)文字對(duì)應(yīng)教程)

后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的高通量数据中批次效应的鉴定和处理（六）- 直接校正表达矩阵的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。