触类旁通！一文掌握泛癌之感兴趣基因集SNV分析

Pancancer一词，在生信文章中出现的频率越来越高，已经成为一个网红分析内容，趁着今天的机会，今天小云想为大家分享一下泛癌分析之感兴趣基因集SNV分析，接下来小云带着大家开始今天的学习之旅。

1. 如何进行泛癌中感兴趣基因集SNV分析？

在进行分析之前，小云想为大家科普一下何为SNV,SNV为单核苷酸多态性，主要包括在coding区域和非coding区域发生突变，编码区的SNV可能影响蛋白质序列，而非编码区的的SNV可能影响基因的表达和剪接过程，这就是小云对SNV的简单介绍哈！

接下来小云为大家介绍一哈咋进行该分析，其实很简单，只需要下载pancer突变数据，以及自己感兴趣的基因集，就可以来统计感兴趣基因集在不同肿瘤组织中发生突变的频率，以及发生的突变类型，小云已经迫不及待想开始实操啦，那就开始吧！

2.需要准备的R包

#安装需要的R包

install.packages("ggplot2")

install.packages("readxl")

install.packages("stringr")

install.packages("dplyr")

install.packages("tidyverse")

install.packages("reshape2")

install.packages("RcolorBrewer")

install.packages("magrittr")

BiocManager::install("maftools")

#导入需要的R包

library(magrittr)

library(readxl)

library(stringr)

library(ggplot2)

library(maftools)

library(dplyr)

library(reshape2)

library(tidyverse)

library(RColorBrewer)

3. 数据准备

在进行感兴趣基因在泛癌中的突变情况时，我们只需要两个文件就可以顺利的完成分析，包括TCGA pancancer突变数据和TCGA pancancer 样本信息文件，小云为大家提供了数据下载网址：https://gdc.cancer.gov/about-data/publications/pancanatlas，下载起来so esay啦!,小云并为大家贴上该网址，方便大家迅速完成数据准备工作。

#TCGA pancancer突变数据下载

mc3.v0.2.8.PUBLIC.maf.gz，泛癌突变数据文件，需要关注基因名，变异类型以及样本名这三列信息。

#泛癌样本信息文件下载

merged_sample_quality_annotations.tsv，只需要关注前三列数据就可以，样本信息和肿瘤类型。

 4. 泛癌中感兴趣基因集SNV分析

如果小伙伴已经下载好了所需要的数据，小云将通过样本信息文件，泛癌突变数据文件和感兴趣基因集，来提取不同肿瘤组织中目标基因集发生的突变文件，统计目标基因发生在不同肿瘤中的突变频率和突变类型，通过该分析小伙伴会掌握很多数据处理的小技巧，小云强烈安利小伙伴认真学习，真的是干活满满。

# 修正TCGA名称

#读取样本信息文件

rawAnno <- read.delim("merged_sample_quality_annotations.tsv",sep = "\t",row.names = NULL,check.names = F,stringsAsFactors = F,header = T)

#提取1-15位字符串

rawAnno$simple_barcode <- substr(rawAnno$aliquot_barcode,1,15)

#去除重复的样本名

samAnno<- rawAnno[!duplicated(rawAnno$simple_barcode),c("cancer type", "simple_barcode")]

#去除cancer type 为空的数据

samAnno <- samAnno[which(samAnno$`cancer type` != ""),]

#保存简化的样本信息文件

write.table(samAnno,"simple_sample_annotation.txt",sep = "\t",row.names = F,col.names = T,quote = F)

# 加载泛癌突变MAF文件，读取文件有点慢，耐心等待奥。

panmaf <- read_tsv("mc3.v0.2.8.PUBLIC.maf.gz", comment = "#")

# 感兴趣基因列表

genelist<- c("PER3","PER2","GPR50","CYP2C19","MTNR1B","CLOCK","CYP1A2","RORA","ARNTL","MTNR1A","ASMT","AANAT")

#codeing区域

vc_nonSyn <- c("Frame_Shift_Del",

               "Frame_Shift_Ins",

               "Splice_Site",

               "Translation_Start_Site",

               "Nonsense_Mutation",

               "Nonstop_Mutation",

               "In_Frame_Del",

               "In_Frame_Ins",

               "Missense_Mutation")

## 非coding区域

vc_Syn <- c("3'UTR",

            "5'UTR",

            "3'Flank",

            "5'Flank",

            "IGR",

            "Intron",

            "Silent")

# 获取不同癌种的样本数

mafsam <- data.frame(samID = panmaf$Tumor_Sample_Barcode,

                     simple_barcode = substr(panmaf$Tumor_Sample_Barcode,1,15),

                     stringsAsFactors = F)

#去除重复的样本

mafsam <- mafsam[!duplicated(mafsam$simple_barcode),]

#合并两个文件

mafsam <- merge(mafsam,samAnno, by = "simple_barcode", all.x = TRUE)

mafsam <- as.data.frame(na.omit(mafsam))

txt <- data.frame(label = paste0(names(table(mafsam$`cancer type`))," (n=", as.numeric(table(mafsam$`cancer type`)),")"),

                  "cancer type" = names(table(mafsam$`cancer type`)),

                  check.names = F) # 将肿瘤和样本数标签合并

#合并两个文件

mafsam <- merge(mafsam, txt, by = "cancer type", all.x = TRUE)

rownames(mafsam) <- mafsam$simple_barcode

# 取感兴趣的突变子集

panmaf <- panmaf[which(panmaf$Hugo_Symbol %in% genelist),]

panmaf$Tumor_Sample_Barcode <- substr(panmaf$Tumor_Sample_Barcode,1,15)

write.table(panmaf, file = "maf_of_genes_of_interest.txt",sep = "\t",row.names = F,col.names = T,quote = F)

# 创建突变矩阵

panmaf <- read_tsv("maf_of_genes_of_interest.txt", comment = "#")

nonSyncount <- matrix(0, nrow = length(genelist), ncol = length(mafsam$simple_barcode),

                   dimnames = list(genelist, mafsam$simple_barcode))

Syncount <- matrix(0, nrow = length(genelist), ncol = length(mafsam$simple_barcode),

                   dimnames = list(genelist, mafsam$simple_barcode))

nonSyncount <- as.data.frame(t(nonSyncount))

Syncount <- as.data.frame(t(Syncount))

#通过循环计算感兴趣基因在每种肿瘤中出现的coding和非coding的数目

for (i in genelist) {

  message(i," starts...")

  submaf <- panmaf[which(panmaf$Hugo_Symbol == i),] # 取出当前基因的突变

  for (j in 1:nrow(submaf)) { # 循环每一行

    tmp <- submaf[j,]

    mutGene <- i # 获取当前突变的基因

    mutSam <- tmp$Tumor_Sample_Barcode # 获取当前突变的样本

    if(is.element(tmp$Variant_Classification, vc_nonSyn)) { # 如果突变类型为coding类

      nonSyncount[mutSam, mutGene] <- 1 # 则coding计数矩阵设为1

    }

    if(is.element(tmp$Variant_Classification, vc_Syn)) { # 如果突变类型为非coding类

      Syncount[mutSam, mutGene] <- 1 # 则非coding计数矩阵设为1

    }

  }

}

# 根据癌种计算累积突变样本量

nonSyncount$tumor <- mafsam[rownames(nonSyncount),"label"]

Syncount$tumor <- mafsam[rownames(Syncount),"label"]

nonSyncount.tumor <- as.data.frame(apply(nonSyncount[,setdiff(colnames(nonSyncount), "tumor")], 2, function(x) tapply(x, INDEX=factor(nonSyncount$tumor), FUN=sum, na.rm=TRUE))) # 对同一种癌种的非沉默突变求和

Syncount.tumor <- as.data.frame(apply(Syncount[,setdiff(colnames(Syncount), "tumor")], 2, function(x) tapply(x, INDEX=factor(Syncount$tumor), FUN=sum, na.rm=TRUE))) # 对同一种癌种的沉默突变求和

mutCount <- matrix(NA, nrow = length(genelist), ncol = length(unique(mafsam$label)),

                   dimnames = list(genelist, unique(mafsam$label)))

mutCount <- as.data.frame(mutCount)

mutFreq <- mutCount

#通过循环计算感兴趣基因在每种肿瘤中出现的coding和非coding的突变频率。

for (i in genelist) {

  for(j in unique(mafsam$label)) {

    tumor <- sapply(strsplit(j," (", fixed = T),"[",1)

    if(nonSyncount.tumor[j,i] > 0) { # 若非沉默突变存在

      mutCount[i,j] <- nonSyncount.tumor[j,i] # 记录该癌种的非沉默突变数

      mutFreq[i,j] <- nonSyncount.tumor[j,i]/as.numeric(table(mafsam$`cancer type`)[tumor]) * 100 # 根据该癌种的样本数计算突变频率（个人认为原文并没有如此计算，原文的颜色映射似乎是按照突变数目来的，而不是突变频率）

    }

    if(nonSyncount.tumor[j,i] == 0 & Syncount.tumor[j,i] > 0) { # 若非沉默突变不存在，而沉默突变存在

      mutCount[i,j] <- 0 # 则记录该癌种的突变数目为0

      mutFreq[i,j] <- 0 # 突变频率为0

    }

    if(nonSyncount.tumor[j,i] == 0 & Syncount.tumor[j,i] == 0) { # 若非沉默突变不存在，沉默突变也不存在

      mutCount[i,j] <- NA # 则此单元为空值

      mutFreq[i,j] <- 0 # 突变频率为0

    }

  }

}

5.绘制热图和瀑布图

#长宽数据转换

ggCount <- melt(as.matrix(mutCount), varnames = c("Cancer", "Gene"), as.is = T)

ggFreq <- melt(as.matrix(mutFreq), varnames = c("Cancer", "Gene"), as.is = T)

#合并两个文件

ggData <- cbind.data.frame(ggCount, Freq = ggFreq$value)

#设置列名

colnames(ggData) <- c("Cancer","Gene","Count","Freq")

#开始绘制热图

p <-

  ggplot(data = ggData) +

  geom_tile(mapping = aes(Gene, Cancer, fill = Freq), col = "grey80") +

  geom_text(mapping = aes(Gene, Cancer, label = Count), size =3) +

  scale_fill_gradient(low = "white","high" = "red") +

  scale_x_discrete(position = "top") +

  labs(fill = "Mutation Frequency (%)") +

  theme_bw() + #theme(axis.text.x = element_text(angle = 45))

  theme(axis.title = element_blank(),

        axis.text.x = element_text(hjust = 0, size = 10, color = "black", angle = 45),

        axis.text.y = element_text(hjust = 1, size = 10, color = "black"),

        axis.ticks = element_line(size=0.2, color="black"),

        axis.ticks.length = unit(0.2, "cm"),

        panel.background = element_blank(),

        panel.grid = element_blank())

ggsave(p, filename = "frequency_heatmap.pdf", width = 10,height = 4)

#绘制瀑布图

# 读取感兴趣的基因子集突变文件

panmaf <- read.maf(maf = "maf of genes of interest.txt",

                   removeDuplicatedVariants = TRUE,

                   isTCGA = FALSE,

                   vc_nonSyn = c("Frame_Shift_Del",

                                 "Frame_Shift_Ins",

                                 "Splice_Site",

                                 "Translation_Start_Site",

                                 "Nonsense_Mutation",

                                 "Nonstop_Mutation",

                                 "In_Frame_Del",

                                 "In_Frame_Ins",

                                 "Missense_Mutation"))

mafAnno<-mafsam[which(mafsam$simple_barcode%in%panmaf@data$Tumor_Sample_Barcode),]

colnames(mafAnno)[1:2] <- c("Cancer","Tumor_Sample_Barcode")

# 设置突变颜色

mutcol <- brewer.pal(n = 10, name = 'Paired')

names(mutcol)<-c('Frame_Shift_Del','Missense_Mutation','Nonsense_Mutation', 'Frame_Shift_Ins','In_Frame_Ins','Splice_Site','In_Frame_Del','Nonstop_Mutation','Translation_Start_Site','Multi_Hit')

cancercolors<-NMF:::ccRamp(brewer.pal(n=12,name='Paired'),length(unique(mafAnno$Cancer)))

names(cancercolors) <- unique(mafAnno$Cancer)

#绘制瀑布图

annocolors = list(Cancer = cancercolors)

pdf(file = "oncoprint.pdf", width = 10, height = 5)

oncoplot(maf = panmaf, # MAF文件

         colors = mutcol, # 突变类型的颜色

         minMut = 0.05 * nrow(mafAnno), # 所显示的突变至少有5%的突变频率

         bgCol = "grey95", # 瀑布图背景色

         borderCol = NA, # 突变的边框（由于样本太多所以不设置边框）

         annotationDat = mafAnno, # 样本注释文件

         annotationColor = annocolors,

         clinicalFeatures = "Cancer", # 从注释文件中显示的信息

         sortByAnnotation = T) # 按照mafAnno的第一列排序，即按照Cancer排序

invisible(dev.off())

6. 结果文件

1. maf_of_genes_of_interest.txt

该文件为提取的感兴趣基因集在泛癌中的突变数据；重点关注基因名，变异类型和样本信息这三列信息。

2. simple_sample_annotation.txt

该文件为简化后的泛癌样本信息文件，第一列样本名，第二列肿瘤类型。

3. frequency_heatmap.pdf

该图为目标基因集在不同肿瘤组织中突变频率热图，横坐标为肿瘤类型以及每种肿瘤所包含的样本数目，纵坐标为目标基因，通过颜色深浅来表示突变频率高低。

4. frequency_heatmap.pdf

该图为绘制的感兴趣基因在不同肿瘤组织中突变情况的瀑布图，利用不同的颜色分别来表示不同肿瘤类型和发生的变异类型。

小云最终绘制了反映感兴趣基因集突变情况的热图和瀑布图，只需要输入自己感兴趣的基因集文件，就可以轻松的完成该分析，泛癌分析的其他分析内容可以尝试使用本公司新开发的云平台生物信息分析小工具，零代码完成分析，云平台网址：http://www.biocloudservice.com/home.html；主要包括泛癌中体细胞拷贝数变异分析(http://www.biocloudservice.com/704/704.php),基因在肿瘤与正常组织中差异分析(http://www.biocloudservice.com/710/710.php)等相关泛癌分析小工具，欢迎大家和小云一起讨论学习奥，下期再见奥。

好了，今天小云的分享就到这里啦，完整版代码上“云生信学生物信息学”公正号回复“代码”领取哟～