对转录组测序的counts矩阵去除批次效应

RNA-seq我们在生信技能树应该是至少推出了400篇教程，而且是我们全国巡讲的标准品知识点，其中还有一个阅读量过两万的综述翻译及其细节知识点的补充：

相信大家听完了我B站的RNA-seq分析流程后，对这个数据的应用方向都不陌生。最近连续收到好几个求助，都是关于转录组测序的counts矩阵去除批次效应，值得写推文解答一下咯！

首先安装R包pasilla

options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos
options()$BioC_mirror
BiocManager::install("pasilla",ask = F,update = F)
 BiocManager::install("DESeq",ask = F,update = F)

简单了解

其实我们只需要使用表达矩阵和表型信息即可

library(pasilla) 
data("pasillaGenes") 

> pasillaGenes
CountDataSet (storageMode: environment)
assayData: 14470 features, 7 samples 
  element names: counts 
protocolData: none
phenoData
  sampleNames: treated1fb treated2fb ... untreated4fb (7 total)
  varLabels: sizeFactor condition type
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
  pubMedIds: 20921232 

Processed data from NCBI Gene Expression Omnibus under accession numbers GSM461176 to GSM461181.

cts=counts(pasillaGenes)
coldata=Biobase::pData(pasillaGenes)
coldata <- coldata[,c("condition","type")]
rownames(coldata) <- sub("fb", "", rownames(coldata))
coldata$type <- sub("-", "_", coldata$type)
coldata$type <- as.factor(coldata$type) 
colnames(cts)=rownames(coldata)
coldata
cts[1:4,1:4]

表达矩阵和表型信息如下：

> coldata
           condition        type
treated1     treated single_read
treated2     treated  paired_end
treated3     treated  paired_end
untreated1 untreated single_read
untreated2 untreated single_read
untreated3 untreated  paired_end
untreated4 untreated  paired_end
> cts[1:4,1:4]
            treated1 treated2 treated3 untreated1
FBgn0000003        0        0        1          0
FBgn0000008       78       46       43         47
FBgn0000014        2        0        0          0
FBgn0000015        1        0        1          0
> 

很明显，我们关心的是treated和untreated两个组的差异而单端测序和双端测序是我们想要抹去的批次效应，因为这篇文章是2010的数据，十年过去了，那个年代NGS是很稀罕的事情，测序仪还在不停的进化，单端测序和双端测序混在一起是容易理解的事情。

这个时候差异分析需要考虑批次效应

很多人以为去除批次效应是要改变你的表达矩阵，新的表达矩阵然后去走差异分析流程，其实大部分的差异分析流程包里面，人家内置好了考虑你的批次效应这样的混杂因素的函数用法设计，比如：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ type+condition)
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, name="condition_untreated_vs_treated")
summary(res)

做出来的差异分析结果是：

out of 11836 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 544, 4.6%
LFC < 0 (down)     : 458, 3.9%
outliers [1]       : 0, 0%
low counts [2]     : 3629, 31%
(mean count < 5)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

可以拿它去绘制火山图，差异基因的热图等等下游分析。

因为我们把批次这个因素，写在了DESeq2里面，所以它在帮我们做差异分析的时候，其实就考虑了，得到的差异分析结果，就是去除了批次效应的。

如果你确实一定要亲眼看看批次效应到底是如何影响这个表达矩阵的，就需要看PCA啦

使用limma的removeBatchEffect去除批次效应

知道注意的是limma的removeBatchEffect这个时候肯定会改变你的counts值矩阵，改变后就没办法走DESeq2差异分析流程啦，仅仅是为了拿到去除批次效应前后对比的表达矩阵而已。

library(ggplot2)
vsd <- vst(dds, blind=FALSE)
pcaData <- plotPCA(vsd, intgroup=c("condition", "type"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape = type)) +
  geom_point(size=3) +
  xlim(-12, 12) +
  ylim(-10, 10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  geom_text(aes(label=name),vjust=2)
ggsave("myPCAWithBatchEffect.png")

assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$type)

pcaData <- plotPCA(vsd, intgroup=c("condition", "type"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape = type)) +
  geom_point(size=3) +
  xlim(-12, 12) +
  ylim(-10, 10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  geom_text(aes(label=name),vjust=2)
ggsave("myPCABatchEffectRemoved.png")

效果如下：

还可以绘制样本相关性热图啦：

library("pheatmap")
library("RColorBrewer")
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pdf("heatmap.pdf", height = 4, width = 5)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)
dev.off()

以上示例代码来源于：https://biohpc.cornell.edu/doc/RNA-Seq-2019-exercise3.html

但是我觉得这些并没有我自己的代码好：免费的数据分析付费的成品代码

总结一下（以下是付费内容）