025 肿瘤突变负荷分析——分析可能被免疫疗法治愈的可能性，突变负荷越大越有可能被免疫疗法治愈

一、生信分析的目的

肿瘤突变负荷也是与免疫治疗相关的一个分析，肿瘤突变负荷越高，免疫治疗的指标就越高。简而言之就是肿瘤突变负荷越大，越有可能被免疫治疗疗法给治愈。

肿瘤突变负荷（TMB）被定义为每百万碱基中被检测出的，体细胞基因编码错误、碱基替换、基因插入或缺失错误的总数。

TMB 是指特定基因组区域内体细胞非同义突变的个数，通常用每兆碱基多少个突变表示（mut/Mb），在早期研究中也直接以突变数量表示。TMB 可以间接反映肿瘤产生新抗原的能力和程度，预测多种肿瘤的免疫治疗疗效。

TMB是指肿瘤细胞基因组中，所评估基因的外显子编码区每兆碱基中发生置换和插入/缺失突变的总数。一方面驱动基因突变可以导致肿瘤的发生，另一方面大量的体细胞突变可以产生新抗原，这些新抗原可以激活CD8+的细胞毒性T细胞，从而发挥T细胞介导的抗肿瘤效应。因此，当基因变异数目累积增多时，就会产生更多的新抗原，进而被免疫系统识别的可能性就越大。TMB最初是在使用Ipilimumab或Tremelimumab治疗晚期黑色素瘤患者时作为预测疗效的生物标记物，结果发现拥有高水平TMB的黑素色瘤和NSCLC患者对PD-1/PD-L1免疫检查点抑制剂疗效往往高于低水平TMB表达的患者。既往研究提示：NSCLC中接受抗PD-1/PD-L1免疫检查点抑制剂治疗的患者，与TMB密切相关的非同义突变数量与肿瘤的客观缓解率（objective response rate, ORR）、持续临床获益时间及无进展生存期（progression-free survival, PFS）的改善均呈正相关

二、具体的代码

1.矩阵提取

去TCGA中下载突变数据
整理好了，就可以得到突变矩阵了：maf_sample_sheet.tsv

remove(list = ls()) ##清空当前环境
setwd("D:/生信代码复现/25.肿瘤突变负荷分析/25.肿瘤突变负荷分析") ##💚💚💚设置路径

  
#这个函数有两个参数，
#@metadata是从TCGA数据下载的sample sheet
#@path是保存maf文件的路径

merge_maf <- function(metadata, path){
  #通过合并path,还有sample sheet前两列得到每一个文件的完整路径
  filenames <- file.path(path, metadata$file_id, metadata$file_name, 
                         fsep = .Platform$file.sep)

  message ('############### Merging maf data ################\n',
           '### This step may take a few minutes ###\n')
    #通过lapply循环去读每一个样本的maf，然后通过rbind合并成矩阵，按行来合并
    #colClasses指定所有列为字符串
    mafMatrix <- do.call("rbind", lapply(filenames, function(fl) 
      read.table(gzfile(fl),header=T,sep="\t",quote="",fill=T,colClasses="character")))
    return (mafMatrix)
}

#定义去除重复样本的函数FilterDuplicate
FilterDuplicate <- function(metadata) {
  filter <- which(duplicated(metadata[,'sample']))
  if (length(filter) != 0) {
    metadata <- metadata[-filter,]
  }
  message (paste('Removed', length(filter), 'samples', sep=' '))
  return (metadata)
}



#读入maf的sample sheet文件
metaMatrix.maf=read.table("maf_sample_sheet.tsv",sep="\t",header=T)
#替换.为下划线，转换成小写，sample_id替换成sample
names(metaMatrix.maf)=gsub("sample_id","sample",gsub("\\.","_",tolower(names(metaMatrix.maf))))
#删掉最后一列sample_type中的空格
metaMatrix.maf$sample_type=gsub(" ","",metaMatrix.maf$sample_type)

#删掉重复的样本
metaMatrix.maf <- FilterDuplicate(metaMatrix.maf)


#调用merge_maf函数合并maf的矩阵
maf_value=merge_maf(metadata=metaMatrix.maf, 
                     path="maf_data"
                     )
#查看前三行前十列
maf_value[1:3,1:10]


#保存合并后的maf文件
write.table(file="combined_maf_value.txt",maf_value,row.names=F,quote=F,sep="\t")

2.突变负荷打分

#安装maftools包
#BiocManager::install("maftools")
library(maftools)

#读取合并的maf文件
laml <- read.maf(maf = "combined_maf_value.txt")   #💛💛💛上面提取出来的突变矩阵

#计算tmb值
tmb_table_wt_log = tmb(maf = laml)
write.table(tmb_table_wt_log,file="TMB_log.txt",sep="\t",row.names=F)   #💜💜💜输出突变负荷打分文件

3.目标基因的肿瘤突变分析

remove(list = ls()) ##清空当前环境
setwd("D:/生信代码复现/25.肿瘤突变负荷分析/25.肿瘤突变负荷分析/3.目标基因肿瘤突变负荷分析") ##💚💚💚设置路径


#引用包
library(limma)
library(ggplot2)
library(ggpubr)
library(ggExtra)

expFile="geneExp.txt"     #💛💛💛我们研究的目标基因的表达数据文件，我们可以通过番外篇：目标基因表达量提取得到这个基因表达文件
tmbFile="TMB.txt"         #💛💛💛肿瘤突变负荷文件，上述的代码分析就可以得到

#读取表达数据文件
rt=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)
gene=colnames(rt)[1]

#删掉正常样品
tumorData=rt[rt$Type=="Tumor",1,drop=F]
tumorData=as.matrix(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3\\", rownames(tumorData))
data=avereps(tumorData)

#读取肿瘤突变负荷的文件
tmb=read.table(tmbFile, header=T, sep="\t", check.names=F, row.names=1)
rownames(tmb)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3\\", rownames(tmb))
tmb=avereps(tmb)
#数据合并并输出结果
sameSample=intersect(row.names(data), row.names(tmb))
data=data[sameSample,,drop=F]
tmb=tmb[sameSample,,drop=F]
rt=cbind(data, tmb)

#相关性分析
x=as.numeric(rt[,gene])
y=log2(as.numeric(rt[,"total_perMB_log"])+1)
df1=as.data.frame(cbind(x,y))
corT=cor.test(x, y, method="spearman")
p1=ggplot(df1, aes(x, y)) + 
			xlab(paste0(gene, " expression"))+ylab("Tumor mutation burden")+
			geom_point()+ geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
			stat_cor(method = 'spearman', aes(x =x, y =y))
p2=ggMarginal(p1, type = "density", xparams = list(fill = "orange"),yparams = list(fill = "blue"))

#输出相关性图形
pdf(file="cor.pdf",width=5,height=5)  #💜💜💜这里最终输出的就是我们研究的基因和肿瘤突变负荷之间的关系
print(p2)
dev.off()

这里就得到了目标基因的表达与肿瘤突变负荷之间的关系，当然了R>0.3和p<0.05才具有意义，这里我是随便调了一个基因给同学们演示。 与肿瘤突变负荷相关性越大则越有理由认为该基因在免疫治疗中起着作用。

三、最终的效果

①R语言中的效果：
image.png

②文件效果：
image.png

四、遇到的问题

出现了部分报错：
image.png

五、参考

①具体代码参考——公众号《叉叉滴同学的生信笔记》