TCGA数据辅助甲基化区域的功能研究
你好啊，我是Christine，这里是周一数据挖掘专栏，欢迎和我们一起学习 ~
如果你不想错过我们的精彩教程，请置顶我们：没看到通知？是不是五行缺星？
如果你不想漏掉我们往期教程，请学会搜索：历史宝藏这样找
本周我们学习的文章是Methylome sequencing in triple-negative breast cancer reveals distinct methylation clusters with prognostic value，2015年发表于Nature Communication。
文章内容概述：
用MBDCap-Seq得到了19个TNBC样本及6个配对正常样本的差异甲基化区域（DMRs）
MBDCap-Seq：特异性捕获与MBD2（Methyl-CpG Binding Domain Protein 2）蛋白结合的DNA区域并测序。
对这些区域相关的基因进行功能注释，用TCGA数据解释这些基因的突变情况、甲基化对表达量的影响
进一步筛选出TCGA-TNBC中特异的甲基化区域
根据甲基化图谱将TCGA-TNBC分为生存差异显著的组，进一步鉴定出17个与生存相关的区域
这篇文章用自己的数据获得一些感兴趣的甲基化区域，后面很大篇幅都是在用TCGA的数据进行验证。（所以非常值得大家学习）
本周我们要重复的图是：Fig3a,b：MBDCap-Seq得到的高甲基化区域将TCGA-TNBC分为3组，且生存差异显著。
Step1. 将MBDCap-Seq高甲基化区域与HM450K探针对应
原文用的是hg18参考基因组，而UCSC xena上用到的甲基化探针是基于hg19基因组，所以先要将补充数据中给的hg18位置信息转为hg19版本。
不同版本的参考基因组差别真的很大，我一开始直接用MBDCap-seq的区域与hg19版本的HM450k做mapping，只得到100多个探针，折腾了很久找不到原因，问了Jimmy老师后才发现是参考坐标不同 (ー_ー)!!
（所以，有一个好的老师很重要！）
rm(list = ls())
options(stringsAsFactors = F)
## MBDCap-Seq的高甲基化区域
# 文件来自文章的补充表1
library(readxl)
hyper_id <- read_excel("./raw_data/ncomms6899-s2.xls",
sheet = "DMRs hyper", skip = 2)[,2:4]
hyper_id$pos <- paste0(hyper_id$chr,":",hyper_id$start,"-",hyper_id$end)
write.table(hyper_id$pos,quote = F,row.names = F,col.names = F,file = "./hyper_pos.txt")
将hyper_pos.txt文件上传到https://genome.ucsc.edu/cgi-bin/hgLiftOver网站，将hg18转为hg19，得到一个bed文件，染色体信息格式为chrN:start-end，为了方便后面处理，把":"改成"-"。
cat hyper_hg19.bed | sed 's/:/-/' > hyper_hg19.txt
找到MBDCap-Seq高甲基化区域对应的HM450K探针
rm(list = ls())
hyper_id <- read.delim(file = "./raw_data/hyper_hg19.txt",header = F,sep = "-",stringsAsFactors = F)
colnames(hyper_id) <- c("chr","start","end")
m450_id <- read.delim("./raw_data/illuminaMethyl450_hg19_GPL16304_TCGAlegacy",header = T,
stringsAsFactors = F)
m450_id <- m450_id[,c(1,3,4,5)]
colnames(m450_id) <- c("probe","chr","start","end")
m450_id <- m450_id[m450_id$chr != "",]
# 计算交集
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("GenomicRanges")
library(GenomicRanges)
#822个hypermethylation区域对应多少个HM450k探针
hyper <- GRanges(hyper_id)
m450 <- GRanges(m450_id)
h2m <- findOverlaps(m450,hyper)
probes <- m450_id[h2m@from,1] #4728个探针，文中为4987
write.table(probes,row.names = F,col.names = "sample", quote = F,file = "./probes.txt")
根据这些探针筛选TCGA-BRCA甲基化矩阵：
从UCSC xena下载TCGA-BRCA的450K甲基化矩阵HumanMethylation450.gz（共746M），解压后提取目标探针数据。
gunzip HumanMethylation450.gz
grep -f probes.txt HumanMethylation450 > hyper_methy450.txt
gzip hyper_methy450.txt
友情提示：如果用windows的小伙伴这里grep没有结果，可以试一试用notepad++打开probes.txt，选择”编辑“ → ”文档格式转换“ →“转换为UNIX格式”
jimmy友情提示，这里的grep -f 用法并不可取，有更高效的方法
Step2. 将TCGA-TNBC根据甲基化聚类
从UCSC xena下载TCGA-BRCA表型数据，提取TNBC样本
rm(list = ls())
# 读取临床表型数据
phe <- read.table("./raw_data/BRCA_clinicalMatrix.gz", sep = "\t", header = T, stringsAsFactors = F)
colnames(phe) #挑选需要的表型信息，这里选择了ER,PR,HER2状态，生存状况及生存时间
phe <- phe[,c(1,8,10,22,32,33)]
rownames(phe) <- gsub("-",".",phe$sampleID)
# 筛选TNBC样本
TNBC_phe <- phe[phe$ER_Status_nature2012 == "Negative" &
phe$HER2_Final_Status_nature2012 == "Negative" &
phe$PR_Status_nature2012 == "Negative",] #123个样本
# 读入甲基化数据
methy <- read.delim("./raw_data/hyper_methy450.txt.gz",header = T,stringsAsFactors = F)
methy[1:4,1:4]
s <- intersect(rownames(TNBC_phe),colnames(methy)) #73个TNBC样本
TNBC_phe <- TNBC_phe[s,]
rownames(methy) <- methy$sample
methy <- methy[,s]
hist(apply(methy,1, function(x){sum(is.na(x))}))
methy[is.na(methy)] <- 0
methy <- as.matrix(methy)
按文章的方法，用ConsensusClusterPlus进行无监督聚类，选择最适合的聚类个数：
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("ConsensusClusterPlus")
library(ConsensusClusterPlus)
res <- ConsensusClusterPlus(d=methy,maxK = 4, reps = 1000,
pItem = 0.8, pFeature = 0.8,
clusterAlg = "km", distance = "euclidean") #按文章中的参数设置
res
calcICL(res,title="untitled_consensus_cluster",plot=NULL,writeTable=FALSE) #计算聚类的一致性
运行结束后会给出一些图，反映每种聚类的效果，分为3类较合适。
Step3.热图及生存分析
#提取聚类结果
clusters <- res[[3]][["consensusClass"]]
table(clusters)
TNBC_phe$group <- clusters[rownames(TNBC_phe)]
#甲基化热图
library(pheatmap)
mat <- t(methy)
TNBC_phe <- TNBC_phe[order(TNBC_phe$group),]
anno <- data.frame(KEY=as.character(TNBC_phe$group))
rownames(anno) <- rownames(TNBC_phe)
mat <- mat[rownames(anno),]
anno_color <- list(KEY=c("1" ="red","2"="orange","3"="blue"))
pheatmap(mat,cluster_rows = F,
clustering_method = "average",
show_colnames = F, show_rownames = F,
annotation_row = anno,
annotation_colors = anno_color,
annotation_legend = F,
color = colorRampPalette(c("RoyalBlue", "white", "Firebrick1"))(100))
每组的样本数和原文不大一致，但从热图中能够粗略地看出3组的甲基化程度确实存在差异（红色-高甲基化，蓝色-低甲基化，橙色-中等）。
# 生存分析
library(survival)
library(survminer)
TNBC_phe$OS.time <- TNBC_phe$OS.time/30
mySurv_OS=with(TNBC_phe,Surv(OS.time, OS))
sfit=survfit(mySurv_OS~group,data=TNBC_phe)
ggsurvplot(sfit,
pval =TRUE,
xlab ="Time in months",
risk.table = T,
palette = c("orange","blue",'red'),
ggtheme =theme_light())
看上去3组的生存是有些差异的，但是不如原文明显，但高、低甲基化组的生存确实能明显区分。
总结一下，虽然没有完全重复出作者的结果，但大致趋势是一致的。过程中也学到了一些新知识：
不同版本间参考基因组坐标转换网站：https://genome.ucsc.edu/cgi-bin/hgLiftOver
用
GenomicRanges包的findOverlaps()函数判断2个染色体区域的重叠情况，这个包主要用于基因组间信息，它的其它功能介绍见R- GenomicRanges使用用
ConsensusClusterPlus依据重抽样方法验证聚类的合理性，用于评估无监督聚类到底分几类合适，结果会给出一些直观的图，方便我们判断。
最后，感谢你看完全文，欢迎加入我们一起学习 ~
千万不要错过了我们下个星期在杭州的巡讲哦！（错过等一年哈）
全国巡讲约你
第1-11站北上广深杭，西安，郑州， 吉林，武汉，成都，港珠澳（全部结束）
一年一度的生信技能树单细胞线下培训班（已结束）