生物信息学习的正确姿势

NGS系列文章包括NGS基础、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集）等内容。

高通量数据中批次效应的鉴定和处理（一）

高通量数据中批次效应的鉴定和处理（二）

高通量数据中批次效应的鉴定和处理（三）- 如何设计尽量避免批次影响

高通量数据中批次效应的鉴定和处理（四）- 在差异基因鉴定过程中移除批次效应

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

# 获取标准化后的表达矩阵并移除低表达基因dat  <- counts(dds, normalized = TRUE)idx  <- rowMeans(dat) > 1dat  <- 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:  2Iteration (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:  2Iteration (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        Sampleuntrt_N61311       untrt     N61311  1.0211325  0.2678603 -0.5299095  untrt_N61311untrt_N052611      untrt    N052611  1.1803986  0.1588371  0.1732030 untrt_N052611untrt_N080611      untrt    N080611  1.1796083 -0.5942965 -0.0832029 untrt_N080611untrt_N061011      untrt    N061011  0.9232642  0.1930937  0.3640127 untrt_N061011trt_N61311           trt     N61311  0.8939275  0.2529656 -0.5620218    trt_N61311trt_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_.trt398                                 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_.trt497                                      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_.trt450                                    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_.trtENSG00000163803    Before_batch_untrt._lowerThan_.trtENSG00000168811    After_sva_batch_untrt._higherThan_.trtENSG00000149218    After_known_batch_untrt._lowerThan_.trtENSG00000168811    After_known_batch_untrt._higherThan_.trtENSG00000118689    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;TypeENSG00000143850;SVA_batch_specificENSG00000065809;SVA_batch_uncorrect_commonENSG00000109689;Uncorrect_specificENSG00000124762;Known_batch_uncorrect_commonENSG00000172061;Known_batch_specificENSG00000163394;Known_batch_SVA_batch_commonENSG00000178695;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;TypeENSG00000144649;SVA_batch_specificENSG00000187134;SVA_batch_uncorrect_commonENSG00000137124;Uncorrect_specificENSG00000151690;Known_batch_uncorrect_commonENSG00000180914;Known_batch_specificENSG00000221866;Known_batch_SVA_batch_commonENSG00000152583;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个去探索的问题，欢迎留言或投稿讨论：

1. 如果不设置只返回两个混杂因素，实际SVA会判断出存在3个混杂因素，全部混杂因素都考虑进去结果会有什么变化呢？

2. 上面是取了单个基因查看其表达模式，还可以进一步比较不同子集的基因表达水平、差异倍数、FDR、差异倍数方差的整体分布，分析受影响的主要是哪些类型的基因？