按基因在染色体上的顺序画差异甲基化热图
你好啊，我是Christine，这里是周一数据挖掘专栏，欢迎和我们一起学习！
本周我们学习画一个差异甲基化热图，数据分析流程和基因表达芯片没有什么差别，特点在于：
甲基化数据太大，对笔记本电脑很不友好，促使我学习了一下怎样用命令行按列筛选数据（PS：甲基化数据分析还是尽量别用个人电脑吧 ~）
学习了给基因注释染色体臂信息，并按基因在染色体上的位置画图
文章介绍
文章标题为Germline BRCA2 mutations drive prostate cancers with distinct evolutionary trajectories，2017年1月发表于Nature Communication
在潜在可治愈的局限性前列腺癌（localized prostate cancer）中，高达30%的患者在接受放疗或外科手术后会发生扩散和转移。2016年1月，加拿大研究人员在Nature发表了题为Genomic hallmarks of localized，non-indolent prostate cancer的文章，解释了这种现象的基因指纹图谱，数据包括：200个全基因组测序、277个外显子测序数据，分析了SNV, 甲基化、基因组重排。同年，该团队发表了这篇文章，阐述为什么BRCA2缺陷的前列腺癌侵袭性更高。
BRCA2-mutant PCa：具有BRCA2遗传变异的前列腺癌
sporadic Pca：散发性前列腺癌
文章主要思路：
分析14名BRCA2-mutant PCa患者的突变、拷贝数变异、甲基化、基因组重排
比较BRCA2-mutant PCa和sporadic Pca的基因组特征，发现BRCA2-mutant PCa具有多种侵袭性疾病的标志
根据CNA和SNV构建克隆进化树，发现BRCA2-mutant PCa的进化轨迹与sporadic Pca不同
本次我们重复的图是Supplementary Figure 9，BRCA2-mutant(n = 10) 和sporadic PCa(n = 100) 的差异甲基化探针。
Step1.下载数据
甲基化芯片数据存放在GSE80685，平台是Illumina HumanMethylation450 BeadChip
样本信息：
BRCA2突变样本（n=12）= 10个样本 + 2个重复
sporadic样本（n=160）= 104 unique samples + 34个样本重复2次 + 8个样本重复3次 + 2个样本重复4次
有104个sporadic样本，图上的意思是作者选择了年龄>50的100名患者，但是我没找到样本的年龄信息，就用了全部104个样本
rm(list = ls())
options(stringsAsFactors = F)
## 下载甲基化芯片数据
library(GEOquery)
gset <- getGEO("GSE80685", destdir="./raw_data/",
AnnotGPL = F,
getGPL = F)
### 查看一下数据
class(gset) #查看数据类型
length(gset) #此处只有1个芯片平台
class(gset[[1]])
gset #芯片平台是GPL13534，172 samples，0 features
### 获得表型信息
a=gset[[1]]
pd <- pData(a)
colnames(pd)
#brca2和sporadic样本(10 BRCA2-mutant, 100 sporadic PCa)
table(pd$`disease state:ch1`) #12个BRCA2，160个Sporadic，因为样本有重复
# 直接去掉重复样本
s <- unique(pd$`^sample:ch1`)
pd <- pd[pd$title %in% s,]
table(pd$`disease state:ch1`) #10个BRCA2，104个Sporadic
save(pd,file = "phenotype.Rdata"）
不知道是不是因为数据太大，getGEO没有下载到甲基化矩阵信息，从网页下载GSE80685_Matrix_processed.txt.gz，压缩文件有594M，解压后会有1.5G，文件很大，先简单探索一下：
# 查看行数，435378行
zcat GSE80685_Matrix_processed.txt.gz | wc -l
# 查看前6行
zcat GSE80685_Matrix_processed.txt.gz |head |less -S
# 查看列数，345列
zcat GSE80685_Matrix_processed.txt.gz | head -1 |awk '{print NF}'
用命令行按列筛选样本
一共172个样本，却有345列，因为每个样本都含有“Detection_Pval”，删掉：
# 为了方便读取，这里先把文件解压
gunzip GSE80685_Matrix_processed.txt.gz
# 先把ID列保存
awk '{print $1}' GSE80685_Matrix_processed.txt > meth_id.txt
# 再取所有的偶数列（其实这一步可以不要，但这个知识点可能以后会用到，这里记录一下）
awk '{for(i=1;i<=NF;i+=2)$i=""}{print $0}' GSE80685_Matrix_processed.txt > meth_data.txt
因为pvalue大部分都是0，筛选之后还有1.3G，还是太大。
数据中有重复样本，名字中含有"."，这里直接把它们删掉：
## 重复样本的名字中含有".",根据样本名只含有字母和数字将其过滤
awk 'NR==1{for(i=1;i<=NF;i++)if($i~/^[A-Z0-9]*$/)a[++c]=i}{for(i=1;i<=c;i++)printf $a[i]"\t";printf "\n"}' meth_data.txt > uniq_matrix.txt
# 把第一列ID信息加回去
paste meth_id.txt uniq_matrix.txt > meth_matrix.txt
筛选完之后还剩855M，读到R里面应该没问题。
rm(list = ls())
options(stringsAsFactors = F)
# 读入甲基化信息
dat <- read.delim("./raw_data/meth_matrix.txt",header = T, stringsAsFactors = F)
dim(dat)
rownames(dat) <- dat[,1]
dat <- dat[,2:115]
dat[1:4,1:4]
## 表型信息
load("phenotype.Rdata")
Step2. 差异分析
芯片数据一般选择limma包
### limma做差异分析
library(limma)
group_list <- pd$`disease state:ch1`
exprSet <- dat
#分组矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
#比较矩阵
cont.matrix <- makeContrasts(contrasts="BRCA2-Sporadic",levels = design)
#差异表达分析
fit <- lmFit(exprSet, design)
fit2 <-contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
#提取结果
tempOutput <- topTable(fit2, coef="BRCA2-Sporadic", n=Inf)
tempOutput <- na.omit(tempOutput)
DEG_limma <- tempOutput[,c(1,5)]
colnames(DEG_limma) <- c('log2FoldChange','p.adj')
logFC_cutoff <- 0.1
DEG_limma$change <- ifelse(DEG_limma$p.adj < 0.05 & abs(DEG_limma$log2FoldChange) > logFC_cutoff,ifelse(DEG_limma$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
table(DEG_limma$change) # UP 2109, down 5612
DEG_limma <- DEG_limma[DEG_limma$change != "NOT",]
diff_matrix <- dat[rownames(DEG_limma),] #保存差异探针的甲基化矩阵
save(diff_matrix,file = "diff_matirx.Rdata")
差异基因数为7721，与文章中写的7445接近，接下来可以画图了。
Step3. 画热图
和我们平时画的热图行列排布相反，文中的热图是以探针为行，以样本为列，且探针按照在染色体上的位置排列，样本按照表型信息和甲基化β值排列。
rm(list=ls())
## 注释探针对应的染色体位置
### 作者用的IlluminaHumanMethylation450kanno.ilmn12.hg19包，用芯片对应的GPL13534注释也是可以的
BiocManager::install("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
head(ann450k)
load("diff_matirx.Rdata")
load("phenotype.Rdata")
diff_matrix[1:4,1:4]
# 按染色体位置顺序给探针排序
df <- as.data.frame(ann450k[rownames(diff_matrix), c("chr","pos")])
## 染色体臂信息：注释信息里没有，Jimmy老师给了个文件
cytoBand <- read.delim("C:/Users/lx/Desktop/week4/raw_data/cytoBand.txt.gz",
header=FALSE, stringsAsFactors=FALSE)
colnames(cytoBand) <- c("chr","start","stop","arms","gene_pos")
cytoBand$arms <- substr(cytoBand$arms,1,1)
library(dplyr)
df_arms <- group_by(cytoBand,chr,arms) %>%
summarise(start=min(start),stop=max(stop))
df$probe <- rownames(df)
df <- merge(df,df_arms,all=T)
df <- df[df$pos > df$start & df$pos < df$stop,]
df$chr <- factor(df$chr,levels = paste0("chr",c(1:22,"X","Y")))
df <- df[order(df$chr,df$pos),] #先按染色体排序，再按染色体上的位置排序
col_g <- df[,c("chr","arms")]
rownames(col_g) <- df$probe
# 按表型和beta值给sample排序
pd$beta <- apply(diff_matrix,2,mean) #计算每个样本中探针β值的平均数
pd <- pd[order(pd$`disease state:ch1`,pd$beta),] #按样本类型和β值给样本排序
col_s <- data.frame("BRCA2_state"=pd$`disease state:ch1`)
rownames(col_s) <- pd$title
## 画热图
library(pheatmap)
# 原文中是以列为基因，以行为样本，表达矩阵需要转置
choose_matrix <- t(diff_matrix[rownames(col_g),rownames(col_s)])
choose_matrix[1:4,1:4]
chr_color <- rainbow(24)
names(chr_color) <- paste0("chr",c(1:22,"X","Y"))
ann_colors = list(arms = c(p="black",q="lightgrey"),
BRCA2_state = c(BRCA2 = "black", Sporadic = "lightgrey"),
chr = chr_color
)
png(width = 800, height = 550,filename = "heatmap3.png")
pheatmap(choose_matrix,
cluster_rows = F,
cluster_cols = F,
show_colnames = F,
show_rownames = F,
annotation_row = col_s,
annotation_col = col_g,
annotation_colors = ann_colors,
color = colorRampPalette(c( "firebrick", "white","navy"))(100)
)
dev.off()
和原文的图比较：
因为在表型信息中没找到Gleason score，这里没画。
数据太多，好像也看不出什么，感觉我的图上甲基化程度要高一点（蓝色多一点）
一点心得：
一个看起来很简单的热图，重复起来比我想象得要麻烦，但过程中也学到许多新东西。学海无涯，欢迎你与我们同行 ~
如果你对上面的图表完全无法理解，那么你可能需要下面的课程：