（未测试）TCGA数据库的癌症甲基化芯片数据重分析

最近生信技能树更新了一系列甲基化数据分析教程，本周我们一起来完成其中的一个作业：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肿瘤样本的。

总结：

1. TCGA的甲基化数据一般都很大，可以的话，在读入R之前先进行一下样本筛选
2. 差异分析之前一定要进行质控检验，本例中就发现了异常样本
3. 更完整的甲基化芯片信号值矩阵差异分析代码见Jimmy老师的github：https://github.com/jmzeng1314/methy_array/