多芯片分析

发表于: 2021-09-17

因为目前合并多个测序、芯片数据集这一块并没有完全统一的标准，方法大概有五六种。公说公有理婆说婆有理，对于我这样的新手来说，最简单的是跟随顶级文章的文章思路或者分析流程和步骤。

于是我选取了一篇欧洲泌尿外科的顶级文章，从这篇文章的补充材料可以看出来：

移除批次效应前

移除批次效应后

作者将所有的芯片测序数据分为affy和illumina两类;
将affy公司的芯片通过标准的affy公司流程分析及combat;
illumina同样是通过标准的illumina公司流程分析及combat;
最后将两类芯片combat到一起，上一步是消除batch effect是不同的cohort，而这一步是affy和illumina两个平台的区别。

加载包

###step1 加载包################################
# =======================================================

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("bladderbatch")

rm(list=ls()) 

library("sva")

options(stringsAsFactors=FALSE)

setwd('D:\\train\\data')

library(bladderbatch)

加载用于分析的数据

# =======================================================
###step2 加载用于分析的数据################################
#包括一个表达量矩阵和一个分组文件
# =======================================================
data(bladderdata)

edata = exprs(bladderEset)


edata[1:5,1:5]

pheno = pData(bladderEset)

pheno[1:5,1:4]

dt <- data.frame(cel=rownames(pheno), pheno)

dt[1:5,1:5]

绘制图片显示以前聚类效果

# =======================================================
###step3 绘制图片显示合并以前的聚类结果##################
# =======================================================

PCA_raw <- prcomp(t(edata), scale = FALSE)

dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
                     Disease = dt$cancer,
                     batch = dt$batch)

dataGG$batch <- as.factor(dataGG$batch)


library(ggplot2)

(qplot(PC1, PC2, data = dataGG, color = Disease,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))


library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = batch,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))


boxplot(edata, target = "core",
        main = "Boxplots of log2-intensities for the raw data")


dist_mat <- dist(t(edata))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow=c(2,1))
plot(clustering, labels = pheno$batch)
plot(clustering, labels = pheno$cancer)

从这些图片看出有明显批次效应，批次效应出现在以下情形:

一个实验不同部分在不同时间完成；
一个实验不同部分不同人完成；
试剂用量、芯片、实验仪器不同；
自己的数据与网站数据混用

Combat参数（prarametric）或者非参数（non-prarametric）的经验贝叶斯框架进行批次效应校正。

移除批次效应

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# =======================================================
####step4 combat移除批次效应 ############################
# =======================================================

library(bladderbatch)

data(bladderdata)

pheno = pData(bladderEset)
# add fake age variable for numeric
pheno$age = c(1:7, rep(1:10, 5))

# write.table(data.frame(cel=rownames(pheno), pheno),
#             row.names=F, quote=F, sep="\t",
#             file="bladder-pheno.txt")

edata = exprs(bladderEset)

# write.table(edata, row.names=T, 
#             quote=F, sep="\t", file="bladder-expr.txt")
# use dataframe instead of matrix

mod = model.matrix(~as.factor(cancer) + age, data=pheno)

t = Sys.time()

# parametric adjustment

cdata = ComBat(dat=edata, batch=as.factor(pheno$batch), 
               mod=mod,par.prior=TRUE, prior.plots=TRUE)

print(Sys.time() - t)

print(cdata[1:5, 1:5])

# write.table(cdata, "r-batch.txt", sep="\t", quote=F)

绘制聚类后的图片

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

# =======================================================
####step5 绘制合并后聚类图片############################
# =======================================================

dt <- data.frame(cel=rownames(pheno), pheno)


PCA_raw <- prcomp(t(cdata), scale = FALSE)


dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
                     Disease = dt$cancer,
                     batch = dt$batch)

dataGG$batch <- as.factor(dataGG$batch)

library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = Disease,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = batch,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

boxplot(cdata, target = "core",
        main = "Boxplots of log2-intensities for the raw data")

dist_mat <- dist(t(cdata))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow=c(2,1))
plot(clustering, labels = pheno$batch)
plot(clustering, labels = pheno$cancer)

前后对比，可以发现不同batch基本不再是泾渭分明，而我们呢关心的cancer和normal泾渭分明了，便于下一步分析。

实例验证 数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

# =======================================================
####step6 实例验证 ################################
# =======================================================
rm(list=ls()) 
library(dplyr)
library(tibble)
library(tidyr)
library("sva")
options(stringsAsFactors=FALSE)
setwd('D:\\train\\data\\multi_micro')
data1 <- read.table('GSE68799.txt', header = T, sep = '\t')
data2 <- read.table('GSE102349.txt', header = T, sep = '\t')
data3 <- read.table(file = 'GSE118719.tsv', sep = '\t', header = TRUE)
data1[1:5,1:5]
data2[1:5,1:5]
data3[1:5,1:5]
data1 <- data1 %>%
  tidyr::separate(GeneID, into=c('geneSymbol', 'chr'), sep='\\_')%>%
  dplyr::select(-chr)
data1[1:5,1:5]
data2$geneSymbol <- data2$Gene.Symbol
data2$Gene.Symbol <- NULL
data3$Geneid <- NULL
a <- intersect(data2$geneSymbol, data1$geneSymbol)
a <- intersect(a, data3$geneSymbol)
data1 <- data1[which(data1$geneSymbol %in% a), ]
data2 <- data2[which(data2$geneSymbol %in% a), ]
data3 <- data3[which(data3$geneSymbol %in% a), ]

data1 <- data1 %>%
  distinct(geneSymbol, .keep_all = T)
data2 <- data2 %>%
  distinct(geneSymbol, .keep_all = T)
data3 <- data3 %>%
  distinct(geneSymbol, .keep_all = T)%>%
  column_to_rownames(var = 'geneSymbol')
data3 <- log2(data3 +1)
data3[1:5,1:5]
data3 <- data3 %>%
  rownames_to_column(var = 'geneSymbol')
data <- merge(data1, data2, by='geneSymbol')

data <- merge(data, data3, by='geneSymbol')%>%
  column_to_rownames(var = 'geneSymbol')
data[1:5,1:5]
metadat <- as.data.frame(colnames(data))
# write.csv(metadat, file = 'metadat.csv')

write.csv(data, file = 'data.csv')

# =======================================================
####step7 移除批次效应 ################################
# =======================================================
setwd('D:\\train\\data\\multi_micro')
rm(list=ls()) 
library("sva")
options(stringsAsFactors=FALSE)
library(bladderbatch)
data(bladderdata)
edata  <- read.csv('data.csv', header = T, row.names = 1)
pheno <- read.csv('metadat.csv', header = T)
dt <- pheno
PCA_raw <- prcomp(t(edata), scale = FALSE)
dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
                     Disease = dt$cancer,
                     batch = dt$batch)

dataGG$batch <- as.factor(dataGG$batch)

library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = Disease,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = batch,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

dist_mat <- dist(t(edata))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow=c(2,1))
plot(clustering, labels = pheno$batch)
plot(clustering, labels = pheno$cancer)


#再做一个分组列，用于批次效应中排除项。
pheno$hasCancer <- as.numeric(pheno$cancer == "T")  
#分组模型
model <- model.matrix(~hasCancer, data = pheno)
t = Sys.time()
cdata <- ComBat(dat = as.matrix(edata),
                       batch = pheno$batch, mod = model)

print(Sys.time() - t)

print(cdata[1:5, 1:5])

PCA_raw <- prcomp(t(cdata), scale = FALSE)

dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
                     Disease = dt$cancer,
                     batch = dt$batch)

dataGG$batch <- as.factor(dataGG$batch)
library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = Disease,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

library(ggplot2)
(qplot(PC1, PC2, data = dataGG, color = batch,
       main = "PCA plot of the raw data (log-transformed)", size = I(2),
       asp = 1.0)
  + scale_colour_brewer(palette = "Set2"))

dist_mat <- dist(t(cdata))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow=c(2,1))
plot(clustering, labels = pheno$batch)
plot(clustering, labels = pheno$cancer)

question：
同属于illumina平台的batch1和batch2的移除批次效应结果很满意，但是affy效果一般；
1.可能是affy数据表达矩阵为原始数据，而illumina为标准化后的数据，所以差距大；
2.这篇顶刊提示我么为什么对affy标准化后才和illumina合并

转载于: https://www.jianshu.com/p/4bc1f842964c