其他
TCGA数据库的癌症甲基化芯片数据重分析
你好啊，我是Christine，这里是生信菜鸟团周一数据挖掘专栏，欢迎和我们一起学习！
最近生信技能树更新了一系列甲基化数据分析教程，本周我们一起来完成其中的一个作业：TCGA数据库的各个癌症甲基化芯片数据重新分析
作业来自文献：Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma，下载TCGA数据库中口腔癌的甲基化芯片信号值矩阵，然后挑选有N-T配对的32个病人的数据进行差异分析，画火山图和热图。
准备数据
Step1. 从UCSC xena(https://xenabrowser.net/datapages/)中下载头颈癌HNSC的甲基化数据HumanMethylation450.gz（数据很大，如果下载失败，那应该是网络的问题），表型数据HNSC_clinicalMatrix.gz，生存数据HNSC_survival.txt.gz。
Step2. 挑选有Normal-Tumor配对的32个OSCC病人
# UCSC的xena浏览器下载指定HNSC的phenotype数据
pd.all <- read.delim("./raw_data/HNSC_clinicalMatrix", header = T, stringsAsFactors = F)
colnames(pd.all)
pd <- pd.all[,c("sampleID","anatomic_neoplasm_subdivision","sample_type")]
table(pd$anatomic_neoplasm_subdivision)
# 挑出OSCC样本（参考文章methods）
oscc <- c("Oral Cavity","Oral Tongue","Buccal Mucosa","Lip","Alveolar Ridge","Hard Palate","Floor of mouth")
pd <- pd[pd$anatomic_neoplasm_subdivision %in% oscc,]
# 挑出normal-tumor配对样本
pd$patient <- substr(pd$sampleID,1,12)
pd <- pd[pd$patient %in% pd$patient[pd$sample_type=="Solid Tissue Normal"],]
rownames(pd) <- pd$sampleID
pd$sample_type <- ifelse(pd$sample_type=="Primary Tumor","Tumor","Normal")
# 保存配对样本（此时有48对）
write.table(pd$sampleID, file = "./raw_data/samples.txt",quote = F, row.names = F, col.names = F)
Step3.缩减甲基化信号矩阵：因为整个HNSC甲基化矩阵有580个样本，读进R会很慢，而且没有必要，所以先用linux命令根据sample.txt筛选出N-T配对的OSCC样本。
# 解压甲基化信号矩阵文件
$ gunzip HumanMethylation450.gz
# 文件内容是这样的
$ head -5 samples.txt
TCGA-CV-5436-01
TCGA-CV-5436-11
TCGA-CV-5442-01
TCGA-CV-5442-11
TCGA-CV-5966-01
$ cat HumanMethylation450 | cut -f 1-6 |head -5
sample TCGA-CN-6022-01 TCGA-CQ-5327-01 TCGA-CV-5435-11 TCGA-HD-7753-01 TCGA-CQ-5325-01
cg13332474 0.0220 0.0264 0.0307 0.0263 0.0218
cg00651829 0.0209 0.4514 0.0214 0.0231 0.0248
cg17027195 0.0443 0.2414 0.0330 0.0361 0.0416
cg09868354 0.0490 0.0546 0.0469 0.0754 0.0516
## 筛选N-T配对的OSCC样本
# 1.获得OSCC样本所在的列
$ cols=($(sed '1!d;s/\t/\n/g' HumanMethylation450 | grep -nf samples.txt | sed 's/:.*$//'))
# 2.筛选出OSCC样本的甲基化信号矩阵
$ cut -f 1$(printf ",%s" "${cols[@]}") HumanMethylation450 > OSCC_NT_methy450k.txt
载入数据为champ对象
## 安装ChAMP包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ChAMP") #有点漫长，可能会报错，把报错的包手动安装一下应该没什么问题
library("ChAMP")
## 读取具有N-T配对的OSCC甲基化信号矩阵
library(data.table)
a <- fread("./raw_data/OSCC_NT_methy450k.txt", data.table = F )
a[1:4,1:4]
rownames(a)=a[,1]
a=a[,-1]
beta=as.matrix(a)
# beta信号值矩阵里面不能有NA值
beta=impute.knn(beta)
betaData=beta$data
betaData=betaData+0.00001
sum(is.na(betaData))
#表型数据
pd <- pd[colnames(betaData),]
table(pd$patient) #此时有些normal样本没有配对的tumor
NT.s <- pd$patient[duplicated(pd$patient)] #得到32对样本，和文章一致了
pd <- pd[pd$patient %in% NT.s,]
betaData <- betaData[,pd$sampleID]
dim(betaData)
# 载入为ChAMP对象
library(ChAMP)
myLoad=champ.filter(beta = betaData ,pd = pd) #这一步已经自动完成了过滤
dim(myLoad$beta)
save(myLoad,file = './Rdata/step1-output.Rdata')
甲基化矩阵质控
rm(list = ls())
load('./Rdata/step1-output.Rdata')
# 数据归一化
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
dim(myNorm)
QC.GUI(beta=myNorm,arraytype="450K") # 显示有NA值
num.na <- apply(myNorm,2,function(x)(sum(is.na(x))))
hist(num.na) # 看来是有样本有问题，个别样本校正后NA值过高
table(num.na) # 61个样本没有NA，3个样本含有超过250k个NA
myNorm <- myNorm[,which(num.na < 250000)] # 删掉异常样本
save(myNorm,file = './Rdata/step2-champ_myNorm.Rdata')
# 主成分分析
library("FactoMineR")
library("factoextra")
dat <- t(myNorm)
pd <- myLoad$pd[colnames(myNorm),] #去掉异常样本
group_list=pd$sample_type
table(group_list)
dat.pca <- PCA(dat, graph = FALSE)
fviz_pca_ind(dat.pca,
geom.ind = "point",
col.ind = group_list,
addEllipses = TRUE,
legend.title = "Groups")
# 热图
cg=names(tail(sort(apply(myNorm,1,sd)),1000))
library(pheatmap)
ac=data.frame(group=group_list)
rownames(ac)=colnames(myNorm)
pheatmap(myNorm[cg,],show_colnames =F,show_rownames = F,
annotation_col=ac,filename = 'heatmap_top1000_sd.png')
dev.off()
# 相关关系矩阵热图
# 组内的样本的相似性应该是要高于组间的！
pheatmap::pheatmap(cor(myNorm[cg,]),
annotation_col = ac,
show_rownames = F,
show_colnames = F)
# 去除异常样本
# 图中看出TCGA-CV-5971-01，TCGA-CV-6953-11，TCGA-CV-6955-11这3个样本有点异常
myNorm <- myNorm[,!(colnames(myNorm) %in% c("TCGA-CV-5971-01","TCGA-CV-6953-11","TCGA-CV-6955-11"))]
pd <- pd[colnames(myNorm),]
save(pd,myNorm,file = "./Rdata/filtered.Rdata")
从主成分分析图可以看出，其实甲基化信号并没有将tumor和normal分得很开，top1000sd的信号信号矩阵和相关关系矩阵热图显示有3个样本看着很奇怪，我把它们删除了。
差异分析
rm(list=ls())
load("./Rdata/filtered.Rdata")
# 差异分析
group_list <- pd$sample_type
myDMP <- champ.DMP(beta = myNorm,pheno=group_list)
head(myDMP[[1]])
save(myDMP,file = './Rdata/step3-output-myDMP.Rdata')
画图
根据差异分析结果画火山图，因为作者的差异分析用的是paired t test，和我的结果差异有点大，我自己选了差异甲基化的cutoff。
df_DMP <- myDMP$Tumor_to_Normal[,1:5]
logFC_cutoff <- 0.45
df_DMP$change <- ifelse(df_DMP$adj.P.Val < 10^-15 & abs(df_DMP$logFC) > logFC_cutoff,
ifelse(df_DMP$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
table(df_DMP$change)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(df_DMP[df_DMP$change =='UP',]) ,
'\nThe number of down gene is ',nrow(df_DMP[df_DMP$change =='DOWN',]))
library(ggplot2)
g <- ggplot(data=df_DMP,
aes(x=logFC, y=-log10(adj.P.Val),
color=change)) +
geom_point(alpha=0.4, size=1) +
theme_set(theme_set(theme_bw(base_size=10)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=10,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red'))
print(g)
作者用了Cox回归从差异分析结果中获得了影响生存的CpGs，画热图
### 挑选影响生存的CpG位点，画热图
# 提取差异甲基化CpGs
choose_gene <- rownames(df_DMP[df_DMP$change != "NOT",])
choose_matrix <- myNorm[choose_gene,]
## 批量生存分析
# 载入生存信息
phe <- read.delim("./raw_data/HNSC_survival.txt.gz",header = T, stringsAsFactors = F)
phe <- phe[phe$sample %in% colnames(choose_matrix),]
phe <- phe[substr(phe$sample,14,15)=="01",]
choose_matrix <- choose_matrix[,phe$sample]
library(survival)
logrankP <- apply(choose_matrix, 1, function(x){
#x <- choose_matrix[1,]
phe$group <- ifelse(x>mean(x),"High","Low")
res <- coxph(Surv(OS.time, OS)~group, data=phe)
beta <- coef(res)
se <- sqrt(diag(vcov(res)))
p.val <- 1 - pchisq((beta/se)^2, 1)
})
names(logrankP) <- rownames(choose_matrix)
table(logrankP<0.05) #151个CpG位点
surv_gene <- names(sort(logrankP))[1:20] #挑了20个生存最显著的DMP，用于画图
# 热图
plot_matrix <- myNorm[surv_gene,]
annotation_col <- data.frame(Sample=myLoad$pd$sample_type ) #热图数据的bar
rownames(annotation_col) <- colnames(plot_matrix)
ann_colors = list(Sample = c(Normal="green", Tumor="red"))
library(pheatmap)
pheatmap(plot_matrix,show_colnames = F,annotation_col = annotation_col,
border_color=NA,
color = colorRampPalette(colors = c("white","navy"))(50),
annotation_colors = ann_colors)
本次作业主要是学习用ChAMP包处理TCGA甲基化数据的流程，用的方法和原文不同，结果差异有点大。原文筛选配对样本的目的是为了进行“配对样本t检验”，我们这里其实可以保留全部OSCC肿瘤样本的。
总结：
TCGA的甲基化数据一般都很大，可以的话，在读入R之前先进行一下样本筛选
差异分析之前一定要进行质控检验，本例中就发现了异常样本
更完整的甲基化芯片信号值矩阵差异分析代码见Jimmy老师的github：https://github.com/jmzeng1314/methy_array/
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
全国巡讲全球听（买一得五），第二期 ，你的生物信息学入门课
生信技能树的2019年终总结 ，你的生物信息学成长宝藏
如果你也想加入我们的甲基化芯片数据处理讨论群，欢迎查看：一个甲基化芯片信号值矩阵差异分析的标准代码，找到小助手申请入群！