maftools
变异数据处理的总体思路如下：
image.png
思维导图
1.数据下载
TCGA的突变数据有4个软件得到的不同版本：
image.png
这个可以在gdc的官网上找到，case选择KIRC，文件类型选择maf即可获得。
image.png
image.png
选择mutect，就一个文件，直接点进去，download就行，下载下来只有一个tar.gz文件，解压放在工作目录下。
tar -xzvf file.tar.gz 解压tar.gz，即可得到一个maf.gz文件。
同样的筛选条件，参考//www.greatytc.com/p/559d9604fcdf下载临床信息数据并整理。
2.数据读取
使用R包maftools读取。
rm(list=ls())
options(stringsAsFactors = F)
library(maftools)
library(dplyr)
project='TCGA_KIRC'
laml = read.maf(maf = 'TCGA.KIRC.mutect.2a8f2c83-8b5e-4987-8dbf-01f7ee24dc26.DR-10.0.somatic.maf.gz')
#> -Reading
#> -Validating
#> -Silent variants: 8383
#> -Summarizing
#> --Mutiple centers found
#> BI;BCM--Possible FLAGS among top ten genes:
#> TTN
#> MUC16
#> HMCN1
#> -Processing clinical data
#> --Missing clinical data
#> -Finished in 2.474s elapsed (2.281s cpu)
laml
#> An object of class MAF
#> ID summary Mean Median
#> 1: NCBI_Build GRCh38 NA NA
#> 2: Center BI;BCM NA NA
#> 3: Samples 336 NA NA
#> 4: nGenes 9444 NA NA
#> 5: Frame_Shift_Del 1732 5.155 4
#> 6: Frame_Shift_Ins 1201 3.574 1
#> 7: In_Frame_Del 238 0.708 0
#> 8: In_Frame_Ins 350 1.042 0
#> 9: Missense_Mutation 12997 38.682 36
#> 10: Nonsense_Mutation 1259 3.747 2
#> 11: Nonstop_Mutation 18 0.054 0
#> 12: Splice_Site 490 1.458 1
#> 13: Translation_Start_Site 25 0.074 0
#> 14: total 18310 54.494 47
maf_df = laml@data
save(laml,maf_df,file = "maf.Rdata")
length(unique(maf_df$Tumor_Sample_Barcode))#统计样本个数
#> [1] 336
length(unique(maf_df$Hugo_Symbol))#统计基因个数
#> [1] 9444
因此，有336个病人,9444个突变基因信息。
3.突变数据的可视化
3.1 plotmafSummary
maftools 自带可视化函数plotmafSummary，可以比较直观的统计maf文件的数据。
#if (as.numeric(dev.cur()) != 1) graphics.off()
plotmafSummary(maf = laml, rmOutlier = TRUE,
#showBarcodes = FALSE,
addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
image.png
就是将maf_df 数据框做了统计，用barplot和boxplot做了可视化。
3.2 突变频谱图
选择突变数量前30的基因
oncoplot(maf = laml, top = 30, fontSize = 0.7)
image.png
下面展开一下这个图的解读
主体热图
一行是一个基因，总共是9444个基因，从中截取了top30；一列是一个样本，总共是336个样本。不同颜色代表不同类型的突变。
右侧条形图
右侧的条形图是每个基因的突变样本数、突变类型和比例
验证一下突变样本数
count(maf_df,Hugo_Symbol,sort = T)
#> Hugo_Symbol n
#> 1: VHL 169
#> 2: PBRM1 148
#> 3: TTN 77
#> 4: SETD2 46
#> 5: BAP1 37
#> ---
#> 9440: ZWINT 1
#> 9441: ZXDA 1
#> 9442: ZXDB 1
#> 9443: ZXDC 1
#> 9444: ZYX 1
结果显示VHL在169样本中突变，样本总数336，所以是49%，以此类推
条形图的颜色是突变类型，以VHL基因为例，他的突变类型分别是：
maf_df %>% filter(Hugo_Symbol=="VHL") %>%
count(Variant_Classification,sort = T)
#> Variant_Classification n
#> 1: Missense_Mutation 60
#> 2: Frame_Shift_Del 41
#> 3: Nonsense_Mutation 27
#> 4: Frame_Shift_Ins 22
#> 5: Splice_Site 16
#> 6: In_Frame_Del 2
#> 7: Nonstop_Mutation 1
顶部条形图
显示每个样本里突变的基因个数，可以看到最高的是那个一枝独秀的1600多。
laml@variants.per.sample %>% head()
#> Tumor_Sample_Barcode Variants
#> 1: TCGA-B8-4143-01A-01D-1806-10 1611
#> 2: TCGA-B0-5098-01A-01D-1421-08 550
#> 3: TCGA-A3-A8OV-01A-11D-A36X-10 120
#> 4: TCGA-CJ-4920-01A-01D-1429-08 117
#> 5: TCGA-CZ-5468-01A-01D-1501-10 102
#> 6: TCGA-B0-5713-01A-11D-1669-08 97
3.2 可以给oncoplot加上一些临床信息
pd 是临床信息数据框，第一列是Tumor_Sample_Barcode，后面几列是各种分组，如gender，age,stage等等。
load("clinical.Rdata")
pd = clinical[,c("bcr_patient_barcode",
"stage_event",
"gender",
"vital_status")]
library(stringr)
pd$stage_event = sapply(str_extract_all(pd$stage_event,"I|V"), paste,collapse = "")
pd[pd==""] = NA
#调整pd的病人id顺序和laml里的样本id顺序一致
sam_id = unique(laml@data$Tumor_Sample_Barcode)
library(stringr)
pd = pd[match(str_sub(sam_id,1,12),
pd$bcr_patient_barcode),]#根据前十二位匹配
pd$Tumor_Sample_Barcode = sam_id
head(pd)
#> bcr_patient_barcode stage_event gender vital_status
#> 482 TCGA-A3-3383 I MALE Alive
#> 438 TCGA-CJ-4920 I FEMALE Dead
#> 142 TCGA-CJ-4908 I MALE Alive
#> 311 TCGA-CZ-5469 II MALE Dead
#> 234 TCGA-CZ-5986 I MALE Alive
#> 202 TCGA-A3-3365 I MALE Alive
#> Tumor_Sample_Barcode
#> 482 TCGA-A3-3383-01A-01D-0966-08
#> 438 TCGA-CJ-4920-01A-01D-1429-08
#> 142 TCGA-CJ-4908-01A-01D-1429-08
#> 311 TCGA-CZ-5469-01A-01D-1501-10
#> 234 TCGA-CZ-5986-01A-11D-1669-08
#> 202 TCGA-A3-3365-01A-01D-0966-08
pd = select(pd,Tumor_Sample_Barcode,everything())
pd = pd[,-2]
head(pd)
#> Tumor_Sample_Barcode stage_event gender vital_status
#> 482 TCGA-A3-3383-01A-01D-0966-08 I MALE Alive
#> 438 TCGA-CJ-4920-01A-01D-1429-08 I FEMALE Dead
#> 142 TCGA-CJ-4908-01A-01D-1429-08 I MALE Alive
#> 311 TCGA-CZ-5469-01A-01D-1501-10 II MALE Dead
#> 234 TCGA-CZ-5986-01A-11D-1669-08 I MALE Alive
#> 202 TCGA-A3-3365-01A-01D-0966-08 I MALE Alive
save(pd,file = "KIRC_pd.Rdata")
laml = read.maf(maf = 'TCGA.KIRC.mutect.2a8f2c83-8b5e-4987-8dbf-01f7ee24dc26.DR-10.0.somatic.maf.gz',
clinicalData = pd)
#> -Reading
#> -Validating
#> -Silent variants: 8383
#> -Summarizing
#> --Mutiple centers found
#> BI;BCM--Possible FLAGS among top ten genes:
#> TTN
#> MUC16
#> HMCN1
#> -Processing clinical data
#> -Finished in 2.022s elapsed (1.889s cpu)
oncoplot(maf = laml,
sortByAnnotation = TRUE,
clinicalFeatures = c("stage_event",
'vital_status',
'gender')
)
#> [1] "stage"
image.png
自定义颜色
col = RColorBrewer::brewer.pal(n = 8,name = 'Set3')
stagecolors = col[1:4]
names(stagecolors) = na.omit(unique(pd$stage_event))
vscolors = col[5:6]
names(vscolors) = unique(pd$vital_status)
gendercolors = col[7:8]
names(gendercolors) = unique(pd$gender)
ancolors = list(stage_event = stagecolors,
vital_status = vscolors,
gender = gendercolors)
oncoplot(maf = laml,
sortByAnnotation = TRUE,
clinicalFeatures = c("stage_event",
'vital_status',
'gender'),
annotationColor = ancolors,
)
image.png
如果瀑布图想要画自己想要的基因，使用帮助文档中的genes参数，将要画的基因组织成向量的形式传递给这个参数
browseVignettes("maftools")#可查看怎么自定义瀑布图的帮助文档html
*生信技能树课程笔记
