去除RNA-seq数据批次效应

整合不同的表达谱数据时，往往需要查看和去除批次效用。可以根据样本的信息：批次，建库方法、测序方法等，作图查看这些因素是否有影响。去除 batch effect 后再研究不同处理下的差异表达基因，可以减少假阳性。

1. 载入有批次效应的数据

rm(list=ls())
setwd("/Users/xxx/scripts/R_scripts")## 1. 载入有批次效应的数据
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
#
# BiocManager::install("pasilla",force= TRUE)library(pasilla)dataFilesDir = system.file("extdata", package = "pasilla", mustWork=TRUE)
pasillaSampleAnno = read.csv(file.path(dataFilesDir, "pasilla_sample_annotation.csv"))
# pasillaSampleAnnocount_df <- read.table(file.path(dataFilesDir,"pasilla_gene_counts.tsv"),sep="\t",header = TRUE, row.names=1)
count_matrix <- as.matrix(count_df)#样品信息数据框，行名为表达矩阵的列名
colData = data.frame(condition = as.factor(c(rep("untreated",4),rep("treated",3))),type = as.factor(c("SR","SR","PE","PE","SR","PE","PE")))
rownames(colData) <- colnames(count_matrix)#dir(system.file("extdata", package = "pasilla", mustWork=TRUE), pattern = ".txt$")
#gffFile = file.path(dataFilesDir, "Dmel.BDGP5.25.62.DEXSeq.chr.gff")# library("DEXSeq")
# dxd = DEXSeqDataSetFromHTSeq(
#   countfiles = file.path(dataFilesDir, paste(pasillaSampleAnno$file, "txt", sep=".")),
#   sampleData = pasillaSampleAnno,
#   design= ~ sample + exon + condition:exon,
#   flattenedfile = gffFile)
# dxd

2. 可视化查看批次效应

## 2. 可视化查看批次效应
require(DESeq2)
library(ggplot2)
dds = DESeqDataSetFromMatrix(countData = count_matrix,colData = colData,design= ~ condition+type)pcaData <- plotPCA(DESeqTransform(dds), intgroup=c("condition", "type"), returnData=TRUE)# plotPCA {DESeq2}
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape = type)) +geom_point(size=3) +# xlim(-4e+05, 5e+05) +# ylim(-4e+05, 2e+05) +xlab(paste0("PC1: ",percentVar[1],"% variance")) +ylab(paste0("PC2: ",percentVar[2],"% variance")) +geom_text(aes(label=name),vjust=2)#ggsave("myPCAWithBatchEffect.png")

3.使用sva包的ComBat 函数去除批次效应

## 3.使用sva包的ComBat 函数去除批次效应
library(sva)
count_mat <- ComBat(dat=count_matrix, batch=colData$type)
# ComBat:Adjust for batch effects using an empirical Bayes framework
count_mat[1:3,1:5]

4. DESeq2包消除批次效应方法

## 4. DESeq2包消除批次效应
# design 设计矩阵中加入引起批次效应的因素(SR or PE)
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_matrix,colData = colData,design = ~ condition+ type)
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, name="condition_untreated_vs_treated")summary(res)
head(res)
resOdered <- res[order(res$padj),]
deg <- as.data.frame(resOdered)
#deg <- na.omit(deg)
dim(deg)
#write.csv(deg,file= "diff_deseq2.csv")# 看看批次效应导致的差异表达基因,
# 这里是测序方法SR(single-end)和PE (paired-end)
batch_res <-  results(dds, name="type_SR_vs_PE")
summary(batch_res)
head(batch_res)
batchResOdered <- res[order(batch_res$padj),]
batchDeg <- as.data.frame(batchResOdered)
batchDeg <- na.omit(batchDeg)
dim(batchDeg)
head(batchDeg)

5. 用limma的removeBatchEffect去除批次效应

## 5. 用limma的removeBatchEffect去除批次效应### 不推荐
require(limma)
require(edgeR)  # DGEList来自edgeR包# design <- model.matrix(~0+colData$condition)
# expr_mat<- removeBatchEffect(count_matrix, batch=colData$type,
#                   design= design)expr_mat<- removeBatchEffect(count_matrix, batch=colData$type)
# table(expr_mat<0) #会出现负数！
expr_mat[which(expr_mat<0)]=0 # 不知是否合理！expr_mat <- round(expr_mat) # read 为整数# removeBatchEffect:This function is not intended to be used
# prior to linear modelling. For linear modelling,
# it is better to include the batch factors in the linear model.# 查看去掉批次效应之后的数据
count_matrix[1:3,1:5]
expr_mat[1:3,1:5]# 接着，查看数据PCA的分布是否消除了批次效应，
# 再用去除批次效应的expr_mat做下游分析
# dds = DESeqDataSetFromMatrix(countData = expr_mat,
#                              colData = colData,
#                              design= ~ condition +type)
#
# pcaData <- plotPCA(DESeqTransform(dds), intgroup=c("condition", "type"),
#                    returnData=TRUE)
#
# # plotPCA {DESeq2}
# percentVar <- round(100 * attr(pcaData, "percentVar"))
# ggplot(pcaData, aes(PC1, PC2, color=condition, shape = type)) +
#   geom_point(size=3) +
#   # xlim(-4e+05, 5e+05) +
#   # ylim(-4e+05, 2e+05) +
#   xlab(paste0("PC1: ",percentVar[1],"% variance")) +
#   ylab(paste0("PC2: ",percentVar[2],"% variance")) +
#   geom_text(aes(label=name),vjust=2)

去除RNA-seq数据批次效应相关推荐

batchelor包去除单细胞RNA-seq数据批次效应
1. 数据集的加载 # ### batchelor ## # if (!require("BiocManager", quietly = TRUE)) # install.pack ...
BatchQC包可视化分析去除组学数据批次效应
BatchQC优势:网页方式对数据操作并进行可视化展示,从而确定是否有批次效应以及去除后的数分布等.初步分析包括数据可视化.层次聚类.主成分分析(PCA)和显著性检验等. 1. 加载包 # 安装 if ...
基因表达数据批次效应去除方法的研究进展
基因表达数据批次效应去除方法的研究进展李飒 , 赵毅强摘要:在组学和大数据时代,整合分析材料相同但时间.平台.方法.技术和实验室等不同批次的表达数据集将成为常态.但是,不同批次数据集由于非生物因素 ...
Brief Bioinform | 农科院深圳基因组所王怡雯组提出一种去除微生物组数据中批次效应的多元算法框架...
PLSDA-batch:去除微生物组数据中批次效应的多元算法框架 PLSDA-batch: a multivariate framework to correct for batch effects ...
MultiBaC包消除不同组学数据之间的批次效应
MultiBaC包: 能够消除在不同批次中生成的不同组学之间的批次效应.因此,MultiBaC是一种利用现有数据集跨组学技术进行元分析的有效工具. 1. 包和演示数据集的加载 # if (!requi ...
高通量数据中批次效应的鉴定和处理（六）- 直接校正表达矩阵
生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 ...
高通量数据中批次效应的鉴定和处理（五）- 预测并校正可能存在的混杂因素...
生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 ...
高通量数据中批次效应的鉴定和处理（三）- 如何设计尽量避免批次影响
生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 ...
高通量数据中批次效应的鉴定和处理（一）
生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 ...

去除RNA-seq数据批次效应

1. 载入有批次效应的数据

2. 可视化查看批次效应

3.使用sva包的ComBat 函数去除批次效应

4. DESeq2包消除批次效应方法

5. 用limma的removeBatchEffect去除批次效应

去除RNA-seq数据批次效应相关推荐

最新文章

热门文章