肿瘤异质性+免疫浸润细胞数据挖掘(可能是最简单的3分文章了)

你好啊，我是Christine，这里是生信菜鸟团周一数据挖掘专栏，欢迎和我们一起学习！

随着肿瘤研究的不断深入，人们逐渐发现肿瘤内部的异质性是个不容忽视的问题，如今大热的单细胞测序技术是研究这个问题的好办法，但是之前那么多bulk测序的数据总不能浪费啊，于是就有了各种推断肿瘤异质性的工具。今天我们要学习的文章Tumor Heterogeneity Correlates with Less Immune Response and Worse Survival in Breast Cancer Patients 就是分析了肿瘤异质性和免疫浸润细胞的关系，今年4月发表在Annals of Surgical Oncology。

全文的分析非常少，而且很容易移植到TCGA的其它癌症，感觉没啥意义，不过数据分析技巧值得学习和训练一下。

   内容提要：

1. 根据体细胞突变数据推断肿瘤的异质性

2. 怎样做出显著的生存分析图

背景介绍

基于拷贝数变异、甲基化、突变等数据有各种推断肿瘤异质性的方法。本文使用的方法是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 <- [e-mail][[e-mail]$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 <- [e-mail]
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")

不出所料，2组的生存差异并不显著。接下来我们借助算法挑选合适的cutoff值：

方案一：本文提到一个网页工具Cutoff finder，你可以上传数据，然后让它帮你选择合适的cutoff值，可以用于做生存分析，ROC曲线。

cutofffinder.png

看起来很好，可是我用网站的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")

splot.png

这下图是好看了，可仔细看看还是有问题的，2组样本数相差有点大啊，这是不是也从侧面说明了我这里的MATH作为生存分析指标是有不合适的呢？至于为什么和作者的结果不一致，我觉得原因很多：

1. 我们的maf文件可能就不一样，得到MATH值也有差异

2. 作者用表达矩阵样本数只有959个，导致我们的样本数也不同

3. 我们用的生存数据可能也是不同的

4. 不同cutoff选择方法的结果也是不同的，作者的意思是他尝试了多个工具

虽然结果不一致，但我觉得cutoff选择的方法还是值得了解一下的，也许你的数据适用呢！

关于生存分析的讨论，在生信技能树出现过几次：

文章剩下的部分是比较不同MATH分组中免疫浸润细胞组成和免疫相关基因的表达，数据下载地址在原文中有提到，代码和Fig1A差不多，小伙伴们可以实战一下啊！

如果你觉得有点困难，那你可能需要复习一下学徒数据挖掘第二期汇总，当然更欢迎直接报名参加Jimmy老师的线下培训班啦 ~

■   ■   ■ 

全国巡讲约你

七月份我们不外出，只专注单细胞！

系统学习单细胞分析，报名生信技能树的线下培训，手慢无。

一年一度的生信技能树单细胞线下培训班火热招生