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

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測(cè)序分析?(重磅綜述：三萬字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)、DNA甲基化分析、重測(cè)序分析、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)

預(yù)測(cè)并校正可能存在的混雜因素

# 獲取標(biāo)準(zhǔn)化后的表達(dá)矩陣并移除低表達(dá)基因 dat <- counts(dds, normalized = TRUE) idx <- rowMeans(dat) > 1 dat <- dat[idx, ] # 根據(jù)關(guān)鍵生物表型構(gòu)建設(shè)計(jì)矩陣 mod <- model.matrix(as.formula(paste0("~ ", design)), colData(dds)) # 構(gòu)建對(duì)照設(shè)計(jì)矩陣 mod0 <- model.matrix(~ 1, colData(dds)) # 指定混雜因素的數(shù)目為 2，也可以讓 sva 自己預(yù)測(cè) svseq <- svaseq(dat, mod, mod0, n.sv = 2)# 每一行為一個(gè)樣品，每一列為對(duì)應(yīng)樣品不同的混雜因素及其影響程度 svseq$svNumber 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 的基礎(chǔ)上對(duì)數(shù)據(jù)做了一個(gè) log 轉(zhuǎn)換；如果處理的是芯片數(shù)據(jù)，通常已經(jīng)做過 log 換，直接使用 sva 即可)。

# 獲取標(biāo)準(zhǔn)化后的表達(dá)矩陣 dat <- normexpr$rlog # 根據(jù)關(guān)鍵生物表型構(gòu)建設(shè)計(jì)矩陣 mod <- model.matrix(as.formula(paste0("~ ", design)), colData(dds)) # 構(gòu)建對(duì)照設(shè)計(jì)矩陣 mod0 <- model.matrix(~ 1, colData(dds)) # 指定混雜因素的數(shù)目為 2，也可以讓 sva 自己預(yù)測(cè) svseq2 <- sva(dat, mod, mod0, n.sv = 2) svseq2$svNumber 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

添加預(yù)測(cè)出的Surrogate variable屬性到dds對(duì)象

dds$SV1 <- svseq$sv[,1] dds$SV2 <- svseq$sv[,2]design(dds) <- as.formula(paste("~ SV1 + SV2 +", design)) # 基于預(yù)測(cè)出的混雜因素再次進(jìn)行分析 dds <- DESeq(dds)

可視化展示預(yù)測(cè)出的Surrogate variable屬性與已知的批次信息的關(guān)系

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

從下圖可以看出，預(yù)測(cè)出的混雜因素SV1, SV2與樣品來源的個(gè)體信息 (individual)還是比較一致的 (N052611與N061011的區(qū)分不明顯)。差異最大的是N080611，這與之前分析的PCA結(jié)果也是吻合的。

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

基于預(yù)測(cè)出的混雜因素再次進(jìn)行差異分析，獲得差異基因文件ehbio.simpler.sva_batch.DESeq2.all.DE和其它可視化圖表（暫時(shí)忽略）。

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

比較批次校正前、已知批次校正后和預(yù)測(cè)的批次校正后差異基因變化

根據(jù)已知批次信息校正后差異基因數(shù)目變多了，上調(diào)多了99個(gè)，下調(diào)多了61個(gè)。根據(jù)預(yù)測(cè)的混雜因素校正后，上調(diào)多了52個(gè)，下調(diào)少了1個(gè)。

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 466de_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 527de_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

畫個(gè)Venn圖，看下哪些基因是新增的差異基因，哪些基因批次校正后沒差異了。

all_de = rbind(de_before_batch, de_after_batch) # 隨機(jī)查看6行，信息代表更全面 all_de[sample(1:nrow(all_de),6),] # 結(jié)果存儲(chǔ)到文件中 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

一個(gè)方式是采用代碼，直接出圖

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")

這里還是采用在線工具h(yuǎn)ttp://www.ehbio.com/test/venn/#/ 來做，能直接獲得每個(gè)子集的基因，準(zhǔn)備在線工具所需的文件，一個(gè)兩列格式的文件：第一列為基因名，第二列為基因的上下調(diào)狀態(tài)。

拷貝文件數(shù)據(jù)到網(wǎng)站數(shù)據(jù)輸入處 （操作就不演示了看上一篇文章）：

從untrt上調(diào)基因Venn圖可以看出，校正已知批次信息后既有新增的untrt上調(diào)差異基因，又丟失了之前的一部分untrt上調(diào)差異基因；校正預(yù)測(cè)的混雜因素后，大部分新增差異基因都與已知批次信息校正后的結(jié)果一致，但新增untrt上調(diào)差異基因少。

從untrt下調(diào)基因Venn圖可以看出，校正預(yù)測(cè)的混雜因素后，新增39個(gè)差異基因；批次校正前鑒定為存在差異的40個(gè)基因在校正后被認(rèn)為是非差異顯著基因。

下面還是從這些基因的表達(dá)模式上看是否可以找到一些線索？

下圖比對(duì)繪出了7種不同類型untrt上調(diào)的差異基因中隨機(jī)選取1個(gè)繪制的表達(dá)模式比較圖。

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代表了差異倍數(shù)特別大的基因，不論是否校正都可以檢測(cè)出差異；不同類型批次信息校正后被檢測(cè)視為差異的基因都有表達(dá)的本底差異；Uncorrect_specific的基因本底表達(dá)無固定模式。

下圖比對(duì)繪出了7種不同類型untrt下調(diào)的差異基因表達(dá)分布，基本結(jié)論與上圖類似。

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")

額外的一個(gè)信息是SVA_batch_speific中紅色和綠色個(gè)體本地表達(dá)區(qū)分不明顯。這可能是基于SVA預(yù)測(cè)的混雜因素與已知的批次因素校正后結(jié)果有差異的一個(gè)原因 (這兩個(gè)個(gè)體的SV值很接近)。

另外一個(gè)導(dǎo)致SVA預(yù)測(cè)的批次與已知的批次效應(yīng)校正后結(jié)果不同的原因也可能是我們只讓SVA預(yù)測(cè)了2個(gè)混雜因素。留下2個(gè)去探索的問題，歡迎留言或投稿討論：

如果不設(shè)置只返回兩個(gè)混雜因素，實(shí)際SVA會(huì)判斷出存在3個(gè)混雜因素，全部混雜因素都考慮進(jìn)去結(jié)果會(huì)有什么變化呢？

上面是取了單個(gè)基因查看其表達(dá)模式，還可以進(jìn)一步比較不同子集的基因表達(dá)水平、差異倍數(shù)、FDR、差異倍數(shù)方差的整體分布，分析受影響的主要是哪些類型的基因？

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

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

總結(jié)

以上是生活随笔為你收集整理的高通量数据中批次效应的鉴定和处理（五）- 预测并校正可能存在的混杂因素...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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