使用ICGC数据库进行肿瘤组织突变分析，绘制瀑布图等

admin管理员组

文章数量:1531792

癌症组织突变分析

---

文章目录

• 癌症组织突变分析
• • ICGC 数据库
  • 下载数据
  • • ICGC下载突变数据
    • genecode网站下载基因注释文件
  • 数据预处理
  • • 读入突变数据
    • 对data文件进行基因注释
  • 瀑布图
  • • R包GenVisR
    • 开始绘画
  • 制作突变矩阵
  • • 读入需要的数据
    • 生成突变矩阵
  • 突变位置可视化
  • 计算基因的突变次数
  • • 注释位置信息
  • 根据突变矩阵做生存分析
  • • 处理donor.KIRC-US.tsv文件
    • 将生存信息匹配到突变矩阵中
    • 生存分析图
  • 计算肿瘤突变负荷
  • • 选择肿瘤样本
    • 读入数据
    • 计算TMB
  • 根据TBM中值的生存分析

ICGC 数据库

• 全称：国际肿瘤基因组协作组
• 主要作用：全面阐明导致全球人类疾病负担的多种癌症中存在的基因组变化。ICGC收集了50种不同癌症类型（或亚型）的肿瘤数据，其中包括基因异常表达，体细胞突变，表观遗传修饰，临床数据等；且数据来自亚洲，澳大利亚，欧洲，北美和南美等
• LIHC数据库一般用作验证
• ICGC数据库
• ICGC收集了50种不同癌症类型（或亚型）的肿瘤数据，其中包括基因异常表达，体细胞突变，表观遗传修饰，临床数据等。
  ICGC包括亚洲、澳大利亚、欧洲、北美和南美17个行政区的89项目，包括25000个癌症基因组

下载数据

ICGC下载突变数据

首先点击DCC Data Releases

选择current，下载最新的数据，当然有特殊需求可以下载之前的数据

然后再点击project就可进入到下图的界面：

中括号内分别是癌型和地区的信息，大家可以根据自己的研究方向进行选择并点击进去

点进去之后可以看到很多不同的数据。不同的癌型不同地区里面的数据类型可能是不一样的。

一般来说会有拷贝数信息，基因变异信息，临床信息，样本信息等等

选择自己需要的文件点击即可下载，这里我选择了一下三个文件：

genecode网站下载基因注释文件

本来是想使用R包biomaRt完成对基因的注释的（ensembl gene ID --> gene SYMBOL）,但是因为这个因为实在是太慢了，太慢了，最后还是决定自己去genecode上下载基因注释文件来完成对基因的注释,并且这个注释文件还有基因的位置信息，在画突变图的时候可以用到。

这样子就下载好了我们需要的基因注释文件啦！

数据预处理

读入突变数据

options(stringsAsFactors = F)
data <- read.delim("E:/daiMa/R/puBu/simple_somatic_mutation.open.LIHC-US.tsv", stringsAsFactors=FALSE)
data <- as.data.frame(data)

## 筛选掉不导致氨基酸改变的突变
data <- data[data$consequence_type != "synonymous_variant",]

对data文件进行基因注释

## 制作id，制作突变矩阵的时候要用
Id <- paste(data$icgc_specimen_id,data$icgc_donor_id,sep = "-")
data[,43] <- Id
colnames(data)[43] <- "ID"




## 读入gtf基因注释文件
GeneAnnotationFile <- read.delim("E:/daiMa/R/puBu/gencode.v34.chr_patch_hapl_scaff.annotation.gtf", header=FALSE, comment.char="#")

library(stringr)
## 取出基因注释文件汇总的基因在染色体位置，起始和终止位点
## 这个在画基因突变位点的时候要用
GeneSite <- GeneAnnotationFile[,c(1,4,5)]
colnames(GeneSite) <- c("chr","start","end")

EnsgToName <- str_split(GeneAnnotationFile[,9],";",simplify = T)[,c(1,3)]
colnames(EnsgToName) <- c("ENSG","NAME")


a <- cbind(EnsgToName,GeneSite)
a <- unique.data.frame(a)
## 因为注释文件中基因有许多重复的信息，这里只要带有gene_name的行
a <- a[grep("gene_name",a[,2]),]


a[,1] <- substring(a[,1],9,23)
a[,2] <- str_split(a[,2]," ",simplify =T)[,3]


## 取出突变数据的ENSG编号
ENSG <- data$gene_affected 
loc <- match(ENSG,a$ENSG)
na_loc <- which(is.na(loc))
data <- data[-na_loc,]
loc <- na.omit(loc)
data[,44:47] <- a[loc,-1]
write.csv(data,"new_data.csv",quote = F,row.names = F)

注释完成后的data文件：

瀑布图

什么是瀑布图呢？请看下图：

在SNP的瀑布图中，横轴表示的是不同的样本，纵轴是基因，填充则代表该基因发生突变，不同的颜色代表不同的突变情况。上面的柱状图是对于每个样本突变情况的统计。

从图中我们可以清楚的看到每个样本中基因具体的突变情况，和所有样本整体的突变情况。

那么如何画出瀑布图呢？

R包GenVisR

GenVisR用来绘制瀑布图，首先支持三种通用的突变数据格式

GenVisR是通过从提供的数据集中提取这三列信息进行绘制瀑布图，所以一定要遵守格式相关的列名要求给列名命名。这三列信息分别是：样本名，基因symbol和突变类型

1.MAF
必须包括列名为"Tumor_Sample_Barcode",“Hugo_Symbol”,“Variant_Classification"的信息
2.MGI
必须包括列名为"sample”,“gene_name”,“trv_type"的信息
3.Custom
必须包括列名为"sample”, “gene”, "variant_class"的信息

TCGA上下载的MAF文件可以直接使用，还有下载的MGI也行，但是要注意这两个版本的输入文件列名一定要符合规定，突变类型也要符合这两种文件规定的类型，否则就会报错（默认格式为MAF）:

Detected an invalid mutation type, valid values for MAF are: Nonsense_Mutation, Frame_Shift_Ins, Frame_Shift_Del, Translation_Start_Site, Splice_Site, Nonstop_Mutation, In_Frame_Ins, In_Frame_Del, Missense_Mutation, 5’Flank, 3’Flank, 5’UTR, 3’UTR, RNA, Intron, IGR, Silent, Targeted_Region, NA

这时候只有两种方法：要么把突变类型转化成要求的类型，但有时候我们下载的数据可能出现上面要求的类型之外的突变类型，那么就只能使用另外一种方法了，那就是使用Custom格式，这是一个自定义的格式，对基因突变类型没有要求

首先将filetype定义为Custom，这个参数允许用户输入自定义的文件。但是列名不同。

waterfall(inputMatrix,fileType = "Custom",variant_class_order = levels(factor(inputMatrix$variant_class)))

其中的variant_class_order参数的意思是规定突变类型在图中图例排序的顺序，一般使用levels(factor(inputMatrix$variant_class))即可。在MAF和MGI中这个顺序是默认的，但是在Custom中，需要设置这个参数，否则会报错

Detected NULL in variant_class_order, this parameter is required if file_type is set to “Custom”

我还遇到了一个神奇的错误：

错误: Insufficient values in manual scale. 21 needed but only 20 provided.

查了一下，错误来自于ggplot2，意思是需要21中颜色，但是只定义了20种颜色，waterfall函数中有一个参数mainPalette可以定义颜色，只要加上这个参数便可

mainPalette = rainbow(21)

顺带一提另一个参数：

plotGenes = gene[1:20]

这个参数是一个用来规定画出哪些基因的字符型向量，再用来规定自己感兴趣的基因的时候非常有用。

参考：GenVisR包：突变类型不在waterfall函数指定范围？

开始绘画

## 制作瀑布图输入文件
inputMatrix <- data.frame(sample=data$icgc_donor_id,
                          variant_class=data$consequence_type,gene=data$NAME)

inputMatrix <- na.omit(inputMatrix)

## 选取自己关注的基因，这里随机选取了20个基因作为自己关注的基因
ForceGene <- sample(inputMatrix$gene,20)
library(GenVisR)
pdf(file="waterfall.pdf",height=6,width=10)
waterfall(inputMatrix,fileType = "Custom",variant_class_order = levels(factor(inputMatrix$variant_class)),mainPalette = rainbow(21),plotGenes = ForceGene)
dev.off()

制作突变矩阵

先看一下已经生成好的突变矩阵，可以看出来行是基因，列是样本，样本名是由“样本ID-病人ID”组成的（就是上面我们注释data文件时生成的IP）。

这个突变矩阵行是所有的基因，列是所有的样本，如果出现在我们刚刚下载的突变数据（data）中的地方我们就标记为“Mutation”，如果没出现在突变数据中，那就是“Wild”。即突变矩阵中标注为Mutation的地方代表着这个基因在这个样本中突变，如果为Wild，则代表没有突变

知道了需要什么数据以后，我们就开始写代码生成该矩阵吧！

读入需要的数据

首先制作突变矩阵需要用到data文件的样本id和gene，这样子我们就能知道哪个样本的哪个基因发生了突变了。

mutation <- data[,c(43,44)]
mutation <- unique.data.frame(mutation)
colnames(mutation) <- c("sampleId","gene")

生成突变矩阵

sampleId <- unique(mutation$sampleId)
gene <- unique(mutation$gene)
mutationMatriix <- matrix(data = "Wild",ncol = length(sampleId),nrow = length(gene))
colnames(mutationMatriix) <- sampleId
row.names(mutationMatriix) <- gene


loc1 <- match(mutation$gene,gene)
loc2 <- match(mutation$sampleId,sampleId)

for(i in 1:nrow(mutation)){
  mutationMatriix[loc1[i],loc2[i]] <- "Mutation"
}

write.csv(mutationMatriix,"突变矩阵.csv",quote = F,row.names = T)

突变位置可视化

计算基因的突变次数

首先需要用到我们刚刚生成的突变矩阵

NumMutationMatrix <-data.frame(Name=row.names(mutationMatriix),
                                Num=unname(rowSums(mutationMatriix == "Mutation")))

生成的数据框包含了每个基因突变的次数（当然这个也可以用来筛选自己关注的基因）：

注释位置信息

GeneSiteLoc <- match(NumMutationMatrix$Name,data$NAME)
NumMutationMatrix <- cbind(NumMutationMatrix,data[GeneSiteLoc,45:47])

library(karyoploteR)

## 为了让画出来的突变位点不会太小，所以添加一个比较大的数值
NumMutationMatrix$end=NumMutationMatrix$end+1000000

## 筛选出突变样本大于5的基因，可以自己定义
NumMutationMatrix=NumMutationMatrix[NumMutationMatrix$Num>5,]
dd=makeGRangesFromDataFrame(NumMutationMatrix)
pdf(file="mutationPos.pdf",width=10,height=7)

kp <- plotKaryotype("hg19", plot.type=1)
y1=log2(NumMutationMatrix$Num+1)
kpHeatmap(kp, dd, y=y1, colors = c("green", "black", "red"), r0=0, r1=0.65)
kpAddBaseNumbers(kp, tick.dist=10000000, minor.tick.dist=1000000)
dev.off()

结果图：

根据突变矩阵做生存分析

处理donor.KIRC-US.tsv文件

将下载好的文件donor.KIRC-US.tsv只留下这四列

将alive设置为0，deceased设置为1

并且我们可以看到活人的都没有生存时间，只有随访时间，我们随后的分析可以把活人的随访时间直接当做生存时间

time <- read.delim("E:/daiMa/R/puBu/donor.LIHC-US.tsv")


## 活人的随访时间就当做是生存时间
na_loc <- is.na(time$donor_survival_time)
time[na_loc,3] <- time[na_loc,4]
time<- time[,-4]
colnames(time) <- c("ID","status","surTime")

处理好的time数据框：

将生存信息匹配到突变矩阵中

mutationMatriix <- read.csv("E:/daiMa/R/puBu/mutationMatriix.csv", row.names=1, stringsAsFactors=FALSE)

## 将突变数量少于10个样本的基因去掉
num <- 10
mutationMatriix <- mutationMatriix[rowSums(mutationMatriix == "Mutation") > num,]


library(stringr)
## 取出突变矩阵中的样本名称
matrixSampleName <- str_split(colnames(mutationMatriix),"\\.",simplify = T)[,2]
## 将生存信息匹配到突变矩阵
sampleLoc <- match(matrixSampleName,time$ID)

na_loc <- which(is.na(sampleLoc))
if(length(na_loc)){
  mutationMatriix <- mutationMatriix[,-na_loc]
  sampleLoc <- na.omit(sampleLoc)
}  
mutationMatriix <- t(mutationMatriix)
mutationMatriix <- cbind(time[sampleLoc,],mutationMatriix)

write.csv(mutationMatriix,"surInput.csv",quote = F,row.names = F)

上面的代码中我们将突变样本小于10的基因删除，因为所有样本有365个样本，其中突变样本却还不到10个，这样子分析起来并不具有意义。

结果图：

有了这个文件我们就能根据基因的突变情况，将样本分为突变和未突变两类样本，对每个基因进行生存分析。

生存分析图

library(survival)
## 生存时间的单位是年
mutationMatriix$surTime <- mutationMatriix$surTime/365


## 用来存放基因生存分析的p值，方便查阅
outTab=data.frame()

rt <- mutationMatriix

for(gene in colnames(rt[,4:ncol(rt)])){
  a=rt[,gene]
  diff=survdiff(Surv(surTime, status) ~a,data = rt)
  pValue=1-pchisq(diff$chisq,df=1)
  outTab=rbind(outTab,cbind(gene=gene,pvalue=pValue))
  if(pValue<0.001){
    pValue=format(pValue, scientific = TRUE)
  }else{
    pValue=sprintf("%.3f",pValue)
  }
  
  fit <- survfit(Surv(surTime, status) ~ a, data = rt)
  summary(fit)
  
  pdf(file=paste(getwd(),"/SurPdf/",gene,".survival.pdf",sep=""),
      width = 5.5,           
      height =5,             
  )
  plot(fit, 
       lwd=2,
       col=c("red","blue"),
       xlab="Time (year)",
       mark.time=T,
       ylim=c(0,1.05),
       ylab="Survival rate",
       main=paste(gene,"(p=", pValue ,")",sep="") )
  legend("topright", 
         c("Mutation","Wild"), 
         lwd=2, 
         col=c("red","blue"))
  dev.off()
}
write.table(outTab,file=paste(getwd(),"/SurPdf/","survival.xls",sep = ""),sep="\t",row.names=F,quote=F)

这串代码将会在SurPdf目录下生成所有基因的生存分析图以及一个包含所有基因p值的csv文件，csv文件将方便大家选择比较显著的基因查看。

计算肿瘤突变负荷

肿瘤突变负荷（TMB），即肿瘤组织在每兆（百万）碱基中突变的总数。通俗来讲就是肿瘤组织的突变密度。肿瘤突变密度越高代表着肿瘤组织和正常组织的差异越大，也就越容易被免疫系统发现，相应的免疫治疗的效果就越好。

通常会检测肿瘤组织的TMB作为免疫治疗的一个生物标志，来判断病人是否适合免疫治疗。

计算TMB的主要思想就是根据我们下载的突变数据计算每个样本出现的突变次数，除以整个人类外显子的碱基长度38兆（这个是我在网上查到的，可能不准确，大家可以自己查找一下）

选择肿瘤样本

因为TMB是肿瘤突变负荷，所以首先我们要先搞清楚哪些样本是肿瘤组织。

首先用WPS打开我们下载的specimen.LIHC-US.tsv文件，并利用筛选功能选出所有肿瘤的样本ID

然后把这些样本名称复制到新建的文件中：tumorSampleName.txt

读入数据

## 读入肿瘤组织文件名
tumorSampleName <- read.table("E:/daiMa/R/puBu/tumorSampleName.txt", quote="\"", comment.char="")
tumorSampleName <- tumorSampleName[,1]
data <- read.delim("E:/daiMa/R/puBu/simple_somatic_mutation.open.LIHC-US.tsv", stringsAsFactors=FALSE)
data <- as.data.frame(data)
## 筛选掉不导致氨基酸改变的突变
data <- data[data$consequence_type != "synonymous_variant",]

计算TMB

TBM <- data.frame()
for (i in tumorSampleName){
  a <- length(which(data$icgc_donor_id == i))/38    ##文章查询到人体全外显子是38M碱基
  TBM <- rbind(TBM,c(i,a))
}

colnames(TBM) <- c("ID","TBM")
write.csv(TBM,"TBM.csv",row.names = F,quote = F)

结果：

​

根据TBM中值的生存分析

## 读入生存数据
time <- read.delim("E:/daiMa/R/puBu/donor.LIHC-US.tsv")
## 活人的随访时间就当做是生存时间
na_loc <- is.na(time$donor_survival_time)
time[na_loc,3] <- time[na_loc,4]
time<- time[,-4]
colnames(time) <- c("ID","status","surTime")
time$surTime <- time$surTime/365


TBMLoc <- match(TBM$ID,time$ID)
na_loc <- which(is.na(TBMLoc))
TBMSur <- TBM
if(length(na_loc)){
  TBMSur <- TBMSur[-na_loc,]
  TBMLoc <- na.omit(TBMLoc)
}
TBMSur <- cbind(time[TBMLoc,],TBMSur[,-1])
colnames(TBMSur) <- c("ID","Status","surTime","TBM")

## 以均值中位数做标签
MeanTBM <- mean(as.numeric(TBMSur$TBM))

label <- as.numeric(TBMSur$TBM > MeanTBM)

library(survival)

diff=survdiff(Surv(surTime,Status) ~label,data = TBMSur)
pValue=1-pchisq(diff$chisq,df=1)
pValue=sprintf("%0.3f",pValue)
fit <- survfit(Surv(surTime, Status) ~ label, data = TBMSur)
summary(fit)


pdf(file="TBM-survival.pdf",
    width=5,
    height=4.5)
plot(fit, 
     lwd=2,
     col=c("red","blue"),
     xlab="Time (year)",
     mark.time=T,
     ylim=c(0,1.05),
     ylab="Survival rate",
     main=paste("TBM","(p=", pValue ,")",sep="") )
legend("topright", 
       c("High TMB","Low TMB"), 
       lwd=2, 
       col=c("red","blue"))
dev.off()

本文标签： 突变瀑布肿瘤组织数据库