肿瘤异质性+免疫浸润细胞数据挖掘(可能是最简单的3分文章了)
你好啊，我是Christine，这里是生信菜鸟团周一数据挖掘专栏，欢迎和我们一起学习！
随着肿瘤研究的不断深入，人们逐渐发现肿瘤内部的异质性是个不容忽视的问题，如今大热的单细胞测序技术是研究这个问题的好办法，但是之前那么多bulk测序的数据总不能浪费啊，于是就有了各种推断肿瘤异质性的工具。今天我们要学习的文章Tumor Heterogeneity Correlates with Less Immune Response and Worse Survival in Breast Cancer Patients 就是分析了肿瘤异质性和免疫浸润细胞的关系，今年4月发表在Annals of Surgical Oncology。
全文的分析非常少，而且很容易移植到TCGA的其它癌症，感觉没啥意义，不过数据分析技巧值得学习和训练一下。
内容提要：
根据体细胞突变数据推断肿瘤的异质性
怎样做出显著的生存分析图
背景介绍
基于拷贝数变异、甲基化、突变等数据有各种推断肿瘤异质性的方法。本文使用的方法是MATH (Mutant allele tumor heterogeneity)，是一种根据突变等位基因分布评估肿瘤内部异质性的方法，原理如下图所示，简单来说就是，如果一群细胞内的每种突变的频率都一样，那说明它很可能是由同一种细胞组成的；反之，如果有的突变频率很高，有的很低，分布非常不一致，那么这些细胞很可能彼此间差异非常大，甚至由不同的细胞类群构成。
详细的原理见Intra-tumor heterogeneity in head and neck cancer and its clinical implications
本周我们重现文章中的Fig1：
MATH值与肿瘤突变负荷，KI67基因的表达及肿瘤内异质性的的相关性（图A）
MATH值在TCGA-BRCA中的生存分析（图B,C）
计算MATH值
Step1. 从TCGAmutations包获得TCGA-BRCA 的体细胞突变maf文件
rm(list=ls())
options(stringsAsFactors = F)
# TCGAmutations包整合了TCGA中全部样本的maf文件
if(F){
devtools::install_github(repo = "PoisonAlien/TCGAmutations")
}
library(TCGAmutations)
tcga_available() #查看可用的数据
# 载入TCGA-BRCA maf文件，并保存
tcga_load(study = "BRCA") # 默认加载经过校正后的MC3 maf文件
save(tcga_brca_mc3, file = "./Rdata/TCGA_BRCA_MC3_maf.Rdata")
Step2. 利用maftools包的inferHeterogeneity函数计算MATH值
因为拷贝数变异（CNV）会带来突变频率的异常变化，该函数可以设置参数忽略CNV区域的变异，当然这样也会使变异数目大大减少，这篇文章没有提到这个问题，我默认它没有忽略CNV区域。
# 计算等位基因突变频率，inferHeterogeneity默认识别"t_vaf"列
library(maftools)
laml <- tcga_brca_mc3
laml@data$t_vaf <- laml@data$t_alt_count/laml@data$t_depth
#删除一些突变数小于10的样本
samples <- laml@variants.per.sample[laml@variants.per.sample$Variants > 10,]
samples <- as.character(samples$Tumor_Sample_Barcode)
# 推断肿瘤异质性：计算MATH值和克隆数
# 如果需要忽略CNV区域，加上segFile参数
# seg <- "./raw_data/CNV_segment" #segment文件目录
BRCA_het <- lapply(samples,function(tsb){
# tsb <- samples[1]
# segFile=file.path(seg, paste0(tsb,".seg.txt"))
het=inferHeterogeneity(maf=laml, tsb = tsb,vafCol="t_vaf")
mathd=het$clusterMeans
mathd=mathd[mathd$cluster !='outlier',] # 去除异常值
#mathd=mathd[mathd$cluster !="CN_altered",] # 去除由拷贝数变异产生的cluster
mathd=as.data.frame(table(mathd$Tumor_Sample_Barcode)) # 计算cluster数目
mathd$MATH <- het$clusterData$MATH[1]
return(mathd)
})
BRCA_het <- do.call(rbind,BRCA_het)
colnames(BRCA_het) <- c("Sample_ID","clusters","MATH")
save(BRCA_het,file="./Rdata/BRCA_MATH_contained_CNV.Rdata")
Boxplot
MATH值高低与突变负荷、KI67基因表达、肿瘤内异质性的boxplot
我没有看到文章中的mutation burden和intra-tumor heterogeneity具体是怎么定义的，所以按我自己的理解来做的
# 文中mutation burdern为样本的突变数
df <- tcga_brca_mc3@variants.per.sample
colnames(df) <- c("Sample_ID","Mutation_burden")
df$Mutation_burden <- log2(df$Mutation_burden)
df <- merge(BRCA_het,df)
# 按MATH值中位数分组
df$group <- ifelse(df$MATH > median(df$MATH),"High","Low")
## heterogeneity信息: 从pan-cancer文章的补充数据下载(来自ABSOLUTE软件的评估结果)
df_h <- read.delim("./raw_data/seg_based_scores.tsv",header = T, stringsAsFactors = F)
df_h$Sample <- substr(df_h$Sample,1,12)
df$heterogeneity <- df_h[match(df$Sample_ID,df_h$Sample),"frac_altered"]
# KI67基因表达量，从UCSC Xena下载gene expression RNAseq
# 只找到 Marker Of Proliferation Ki-67（MKI67)
expr <- read.delim("./raw_data/HiSeqV2.gz",header = T, stringsAsFactors = F)
KI67 <-as.numeric(expr[expr$sample == "MKI67",2:1219])
#调整样本名，把"."替换成"-"，并取前12个字符
names(KI67) <- gsub("\\.","-",substr(colnames(expr[2:1219]),1,12))
df$KI67 <- KI67[df$Sample_ID]
# 画图
library(ggpubr)
p1 <- ggboxplot(df, x="group",y="Mutation_burden",
fill = "group",
palette = "jco",
legend = "") +
stat_compare_means(method = "t.test")
p2 <- ggboxplot(df, x="group",y="KI67",
fill = "group",
palette = "jco",
legend = "") +
stat_compare_means(method = "t.test")
p3 <- ggboxplot(df, x="group",y="heterogeneity",
fill = "group",
palette = "jco",
legend = "") +
stat_compare_means(method = "t.test")
ggarrange(p1,p2,p3,ncol = 3,nrow = 1)
KI67基因的表达结果图与文中不一致，可能是因为我们表达矩阵的处理不同吧！
生存分析图
重复别人的生存分析图简直就是噩梦，因为不知道别人的分组究竟是什么样的，我几乎就没有成功过。但这篇文章的方法倒是给了我一个启示：可以在cutoff值的选择上做做文章。
先看看直接按中位数分组的生存分析图：
# 临床数据，从UCSC Xena上下载BRCA的临床数据
phe <- read.delim("./raw_data/BRCA_clinicalMatrix.gz",header = T, stringsAsFactors = F)
phe <- phe[,c("sampleID","OS.time","OS","ER_Status_nature2012",
"PR_Status_nature2012","HER2_Final_Status_nature2012")]
colnames(phe) <- c("Sample_ID","OS.time","OS","ER","PR","HER2")
phe <- phe[substr(phe$Sample_ID,14,15) == "01",] #只保留原发肿瘤样本
phe$Sample_ID <- substr(phe$Sample_ID,1,12)
df <- merge(df,phe)
df$OS.time <- df$OS.time/30
# 全部样本的生存分析
library(survival)
library(survminer)
sfit <- survfit(Surv(OS.time, OS) ~ group, data=df)
# 直接按中位数分组，生存不显著
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
risk.table =TRUE,pval =TRUE,
xlab ="Time in months",
ggtheme =theme_light()
) +
labs(title = "TCGA whole cohort grouped by median MATH")
方案一：本文提到一个网页工具Cutoff finder，你可以上传数据，然后让它帮你选择合适的cutoff值，可以用于做生存分析，ROC曲线。
看起来很好，可是我用网站的example data试了好几次，都是“unkown error”，也不知道是不是我的网络问题，放弃了，感兴趣的小伙伴可以试一下。
方案二：survminer包中自带的surv_cutpoint函数
## 挑选合适的cutoff
cutoff <- surv_cutpoint(df,time="OS.time",event = "OS",variables = "MATH")
cutoff #查看结果
plot(cutoff)
# 按最佳cutoff分组
phe_all <- surv_categorize(cutoff, variables = "MATH", labels = c("low", "high"))
sfit <- survfit(Surv(OS.time, OS) ~ MATH, data=phe_all)
p_sur1 <- ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
risk.table =TRUE,pval =TRUE,
xlab ="Time in months",
ggtheme =theme_light()) +
labs(title = "TCGA whole cohort")
## ER组的生存曲线
table(df$ER)
df_er <- df[df$ER == 'Positive',]
cutoff <- surv_cutpoint(df_er,time="OS.time",event = "OS",variables = "MATH")
phe_er <- surv_categorize(cutoff, variables = "MATH", labels = c("low", "high"))
sfit <- survfit(Surv(OS.time, OS) ~ MATH, data=phe_er)
p_sur2 <- ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
risk.table =TRUE,pval =TRUE,
xlab ="Time in months",
ggtheme =theme_light()) +
labs(title = "ER+")
## PR组的生存曲线
table(df$PR)
df_pr <- df[df$PR == 'Positive',]
cutoff <- surv_cutpoint(df_pr,time="OS.time",event = "OS",variables = "MATH")
phe_pr <- surv_categorize(cutoff, variables = "MATH", labels = c("low", "high"))
sfit <- survfit(Surv(OS.time, OS) ~ MATH, data=phe_pr)
p_sur3 <- ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
risk.table =TRUE,pval =TRUE,
xlab ="Time in months",
ggtheme =theme_light()) +
labs(title = "PR+")
## non-TNBC组的生存曲线
table(df$ER,df$PR,df$HER2)
df_NT <- df[!(df$ER =='Negative' & df$PR =='Negative' & df$HER2 =='Negative'),]
cutoff <- surv_cutpoint(df_NT,time="OS.time",event = "OS",variables = "MATH")
phe_NT <- surv_categorize(cutoff, variables = "MATH", labels = c("low", "high"))
sfit <- survfit(Surv(OS.time, OS) ~ MATH, data=phe_NT)
p_sur4 <- ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
risk.table =TRUE,pval =TRUE,
xlab ="Time in months",
ggtheme =theme_light()) +
labs(title = "Non-TNBC")
splot <- arrange_ggsurvplots(list(p_sur1,p_sur2,p_sur3,p_sur4), print = FALSE, ncol = 2, nrow = 2,
risk.table.height = 0.3)
ggsave("splot.png",plot = splot, width = 25,height = 20, units = "cm")
这下图是好看了，可仔细看看还是有问题的，2组样本数相差有点大啊，这是不是也从侧面说明了我这里的MATH作为生存分析指标是有不合适的呢？至于为什么和作者的结果不一致，我觉得原因很多：
我们的maf文件可能就不一样，得到MATH值也有差异
作者用表达矩阵样本数只有959个，导致我们的样本数也不同
我们用的生存数据可能也是不同的
不同cutoff选择方法的结果也是不同的，作者的意思是他尝试了多个工具
虽然结果不一致，但我觉得cutoff选择的方法还是值得了解一下的，也许你的数据适用呢！
关于生存分析的讨论，在生信技能树出现过几次：
文章剩下的部分是比较不同MATH分组中免疫浸润细胞组成和免疫相关基因的表达，数据下载地址在原文中有提到，代码和Fig1A差不多，小伙伴们可以实战一下啊！
如果你觉得有点困难，那你可能需要复习一下学徒数据挖掘第二期汇总，当然更欢迎直接报名参加Jimmy老师的线下培训班啦 ~
■ ■ ■
全国巡讲约你
第1-10站北上广深杭，西安，郑州， 吉林，武汉和成都（全部结束）
七月份我们不外出，只专注单细胞！
系统学习单细胞分析，报名生信技能树的线下培训，手慢无
全国巡讲第12站-北京（生信技能树爆款入门课）（下一站杭州）