高通量数据中批次效应的鉴定和处理（五）- 预测并校正可能存在的混杂因素

生信宝典

发布于 2020-08-06 14:54:40

7570

发布于 2020-08-06 14:54:40

预测并校正可能存在的混杂因素

# 获取标准化后的表达矩阵并移除低表达基因
dat  <- counts(dds, normalized = TRUE)
idx  <- rowMeans(dat) > 1
dat  <- dat[idx, ]
# 根据关键生物表型构建设计矩阵
mod  <- model.matrix(as.formula(paste0("~ ", design)), colData(dds))
# 构建对照设计矩阵
mod0 <- model.matrix(~ 1, colData(dds))
# 指定混杂因素的数目为 2，也可以让 sva 自己预测
svseq <- svaseq(dat, mod, mod0, n.sv = 2)

# 每一行为一个样品，每一列为对应样品不同的混杂因素及其影响程度
svseq$sv

Number of significant surrogate variables is:  2
Iteration (out of 5 ):1  2  3  4  5
      [,1]      [,2]
[1,]  0.2678603 -0.52990953
[2,]  0.1588371  0.17320301
[3,] -0.5942965 -0.08320290
[4,]  0.1930937  0.36401274
[5,]  0.2529656 -0.56202177
[6,]  0.1750282  0.27665277
[7,] -0.6236673 -0.03396788
[8,]  0.1701789  0.39523355

下面的方式也可以 (svaseq 是在 sva 的基础上对数据做了一个 log 转换；如果处理的是芯片数据，通常已经做过 log 换，直接使用 sva 即可)。

# 获取标准化后的表达矩阵
dat <- normexpr$rlog
# 根据关键生物表型构建设计矩阵
mod  <- model.matrix(as.formula(paste0("~ ", design)), colData(dds))
# 构建对照设计矩阵
mod0 <- model.matrix(~ 1, colData(dds))
# 指定混杂因素的数目为 2，也可以让 sva 自己预测
svseq2 <- sva(dat, mod, mod0, n.sv = 2)
svseq2$sv

Number of significant surrogate variables is:  2
Iteration (out of 5 ):1  2  3  4  5
      [,1]      [,2]
[1,]  0.2500285 -0.52880173
[2,]  0.1689412  0.21277897
[3,] -0.5718486 -0.06104912
[4,]  0.1763780  0.34073904
[5,]  0.2397509 -0.58308079
[6,]  0.1921676  0.26715457
[7,] -0.6476571 -0.02618922
[8,]  0.1922394  0.37844828

添加预测出的Surrogate variable属性到dds对象

dds$SV1 <- svseq$sv[,1]
dds$SV2 <- svseq$sv[,2]

design(dds) <- as.formula(paste("~ SV1 + SV2 +", design))
# 基于预测出的混杂因素再次进行分析
dds <- DESeq(dds)

可视化展示预测出的Surrogate variable属性与已知的批次信息的关系

plot_data <- as.data.frame(colData(dds))
plot_data$Sample <- rownames(plot_data)
head(plot_data)

              conditions individual sizeFactor        SV1        SV2        Sample
untrt_N61311       untrt     N61311  1.0211325  0.2678603 -0.5299095  untrt_N61311
untrt_N052611      untrt    N052611  1.1803986  0.1588371  0.1732030 untrt_N052611
untrt_N080611      untrt    N080611  1.1796083 -0.5942965 -0.0832029 untrt_N080611
untrt_N061011      untrt    N061011  0.9232642  0.1930937  0.3640127 untrt_N061011
trt_N61311           trt     N61311  0.8939275  0.2529656 -0.5620218    trt_N61311
trt_N052611          trt    N052611  0.6709229  0.1750282  0.2766528   trt_N052611

从下图可以看出，预测出的混杂因素SV1, SV2与样品来源的个体信息 (individual)还是比较一致的 (N052611与N061011的区分不明显)。差异最大的是N080611，这与之前分析的PCA结果也是吻合的。

ggplot(plot_data, aes(x=SV1, y=SV2, color=conditions, shape=individual)) +
    geom_point() + geom_text_repel(aes(label=Sample))

基于预测出的混杂因素再次进行差异分析，获得差异基因文件ehbio.simpler.sva_batch.DESeq2.all.DE和其它可视化图表（暂时忽略）。

multipleGroupDEgenes(dds, design=design, output_prefix=output_prefix, padj=padj, log2FC=log2FC)

比较批次校正前、已知批次校正后和预测的批次校正后差异基因变化

根据已知批次信息校正后差异基因数目变多了，上调多了99个，下调多了61个。根据预测的混杂因素校正后，上调多了52个，下调少了1个。

de_before_batch = sp_readTable("ehbio.simpler.DESeq2.all.DE", header=F)
de_before_batch$V2 = paste("Before_batch",de_before_batch$V2,sep="_")
table(de_before_batch$V2)

Before_batch_untrt._higherThan_.trt  Before_batch_untrt._lowerThan_.trt
                                398                                 466

de_after_known_batch = sp_readTable("ehbio.simpler.batch.DESeq2.all.DE", header=F)
de_after_known_batch$V2 = paste("After_known_batch",de_after_known_batch$V2,sep="_")
table(de_after_known_batch$V2)

After_known_batch_untrt._higherThan_.trt  After_known_batch_untrt._lowerThan_.trt
                                     497                                      527

de_after_sva_batch = sp_readTable("ehbio.simpler.sva_batch.DESeq2.all.DE", header=F)
de_after_sva_batch$V2 = paste("After_sva_batch",de_after_sva_batch$V2,sep="_")
table(de_after_sva_batch$V2)

After_sva_batch_untrt._higherThan_.trt  After_sva_batch_untrt._lowerThan_.trt
                                   450                                    465

画个Venn图，看下哪些基因是新增的差异基因，哪些基因批次校正后没差异了。

all_de = rbind(de_before_batch, de_after_batch)
# 随机查看6行，信息代表更全面
all_de[sample(1:nrow(all_de),6),]
# 结果存储到文件中
sp_writeTable(all_de, file="Compare_de_gene_beofore_and_after_batch.txt", keep_rownames = F, col.names = F)

ENSG00000260455    After_known_batch_untrt._higherThan_.trt
ENSG00000163803    Before_batch_untrt._lowerThan_.trt
ENSG00000168811    After_sva_batch_untrt._higherThan_.trt
ENSG00000149218    After_known_batch_untrt._lowerThan_.trt
ENSG00000168811    After_known_batch_untrt._higherThan_.trt
ENSG00000118689    After_known_batch_untrt._lowerThan_.trt

一个方式是采用代码，直接出图

suppressMessages(library(VennDiagram))
suppressMessages(library(grid))
sp_vennDiagram(all_de, label1="Before_batch_untrt._higherThan_.trt",
               label2="After_known_batch_untrt._higherThan_.trt",
               label3="After_sva_batch_untrt._higherThan_.trt")

这里还是采用在线工具http://www.ehbio.com/test/venn/#/ 来做，能直接获得每个子集的基因，准备在线工具所需的文件，一个两列格式的文件：第一列为基因名，第二列为基因的上下调状态。

拷贝文件数据到网站数据输入处：

从untrt上调基因Venn图可以看出，校正已知批次信息后既有新增的untrt上调差异基因，又丢失了之前的一部分untrt上调差异基因；校正预测的混杂因素后，大部分新增差异基因都与已知批次信息校正后的结果一致，但新增untrt上调差异基因少。

从untrt下调基因Venn图可以看出，校正预测的混杂因素后，新增39个差异基因；批次校正前鉴定为存在差异的40个基因在校正后被认为是非差异显著基因。

下面还是从这些基因的表达模式上看是否可以找到一些线索？

下图比对绘出了7种不同类型untrt上调的差异基因中随机选取1个绘制的表达模式比较图。

untrt_up_genes <- "Name;Type
ENSG00000143850;SVA_batch_specific
ENSG00000065809;SVA_batch_uncorrect_common
ENSG00000109689;Uncorrect_specific
ENSG00000124762;Known_batch_uncorrect_common
ENSG00000172061;Known_batch_specific
ENSG00000163394;Known_batch_SVA_batch_common
ENSG00000178695;All_common"

untrt_up_genes <- read.table(text=untrt_up_genes, sep=";", header=T, row.names=NULL)

untrt_up_genes_expr <- merge(untrt_up_genes, normexpr$rlog, by.x="Name", by.y=0, all.x=T)

untrt_up_genes_expr_long <- reshape2::melt(untrt_up_genes_expr, id_vars=c("Name","Type"),
                                 variable.name="Sample", value.name = "Expr")

head(untrt_up_genes_expr_long)

metadata$Sample = rownames(metadata)

sp_boxplot(untrt_up_genes_expr_long, melted=T, metadata=metadata,
           xvariable = "conditions", yvariable = "Expr", jitter_bp = T,
           group_variable_for_line = "individual",
           facet = "Type", scales="free_y", legend.position = c(0.5,0.1),
           x_label="",manual_color_vector = "Set2") +
  theme(legend.direction = "horizontal")

All_common代表了差异倍数特别大的基因，不论是否校正都可以检测出差异；不同类型批次信息校正后被检测视为差异的基因都有表达的本底差异；Uncorrect_specific的基因本底表达无固定模式。

下图比对绘出了7种不同类型untrt下调的差异基因表达分布，基本结论与上图类似。

untrt_down_genes <- "Name;Type
ENSG00000144649;SVA_batch_specific
ENSG00000187134;SVA_batch_uncorrect_common
ENSG00000137124;Uncorrect_specific
ENSG00000151690;Known_batch_uncorrect_common
ENSG00000180914;Known_batch_specific
ENSG00000221866;Known_batch_SVA_batch_common
ENSG00000152583;All_common"

untrt_down_genes <- read.table(text=untrt_down_genes, sep=";", header=T, row.names=NULL)

untrt_down_genes_expr <- merge(untrt_down_genes, normexpr$rlog, by.x="Name", by.y=0, all.x=T)

untrt_down_genes_expr_long <- reshape2::melt(untrt_down_genes_expr, id_vars=c("Name","Type"),
                                 variable.name="Sample", value.name = "Expr")

head(untrt_down_genes_expr_long)

metadata$Sample = rownames(metadata)

sp_boxplot(untrt_down_genes_expr_long, melted=T, metadata=metadata,
           xvariable = "conditions", yvariable = "Expr", jitter_bp = T,
           group_variable_for_line = "individual",
           facet = "Type", scales="free_y", legend.position = c(0.5,0.1),
           x_label="",manual_color_vector = "Set2") +
  theme(legend.direction = "horizontal")

额外的一个信息是SVA_batch_speific中红色和绿色个体本地表达区分不明显。这可能是基于SVA预测的混杂因素与已知的批次因素校正后结果有差异的一个原因 (这两个个体的SV值很接近)。

另外一个导致SVA预测的批次与已知的批次效应校正后结果不同的原因也可能是我们只让SVA预测了2个混杂因素。留下2个去探索的问题，欢迎留言或投稿讨论：

如果不设置只返回两个混杂因素，实际SVA会判断出存在3个混杂因素，全部混杂因素都考虑进去结果会有什么变化呢？
上面是取了单个基因查看其表达模式，还可以进一步比较不同子集的基因表达水平、差异倍数、FDR、差异倍数方差的整体分布，分析受影响的主要是哪些类型的基因？

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-08-03，如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自生信宝典微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度