maftools
變異數(shù)據(jù)處理的總體思路如下：
image.png
思維導(dǎo)圖
1.數(shù)據(jù)下載
TCGA的突變數(shù)據(jù)有4個(gè)軟件得到的不同版本：
image.png
這個(gè)可以在gdc的官網(wǎng)上找到，case選擇KIRC，文件類型選擇maf即可獲得。
image.png
image.png
選擇mutect，就一個(gè)文件，直接點(diǎn)進(jìn)去，download就行，下載下來只有一個(gè)tar.gz文件，解壓放在工作目錄下。
tar -xzvf file.tar.gz 解壓tar.gz，即可得到一個(gè)maf.gz文件。
同樣的篩選條件，參考http://m.itdecent.cn/p/559d9604fcdf下載臨床信息數(shù)據(jù)并整理。
2.數(shù)據(jù)讀取
使用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))#統(tǒng)計(jì)樣本個(gè)數(shù)
#> [1] 336
length(unique(maf_df$Hugo_Symbol))#統(tǒng)計(jì)基因個(gè)數(shù)
#> [1] 9444
因此，有336個(gè)病人,9444個(gè)突變基因信息。
3.突變數(shù)據(jù)的可視化
3.1 plotmafSummary
maftools 自帶可視化函數(shù)plotmafSummary，可以比較直觀的統(tǒng)計(jì)maf文件的數(shù)據(jù)。
#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 數(shù)據(jù)框做了統(tǒng)計(jì)，用barplot和boxplot做了可視化。
3.2 突變頻譜圖
選擇突變數(shù)量前30的基因
oncoplot(maf = laml, top = 30, fontSize = 0.7)
image.png
下面展開一下這個(gè)圖的解讀
主體熱圖
一行是一個(gè)基因，總共是9444個(gè)基因，從中截取了top30；一列是一個(gè)樣本，總共是336個(gè)樣本。不同顏色代表不同類型的突變。
右側(cè)條形圖
右側(cè)的條形圖是每個(gè)基因的突變樣本數(shù)、突變類型和比例
驗(yàn)證一下突變樣本數(shù)
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
結(jié)果顯示VHL在169樣本中突變，樣本總數(shù)336，所以是49%，以此類推
條形圖的顏色是突變類型，以VHL基因?yàn)槔?，他的突變類型分別是：
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
頂部條形圖
顯示每個(gè)樣本里突變的基因個(gè)數(shù)，可以看到最高的是那個(gè)一枝獨(dú)秀的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 是臨床信息數(shù)據(jù)框，第一列是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
#調(diào)整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),]#根據(jù)前十二位匹配
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參數(shù)，將要畫的基因組織成向量的形式傳遞給這個(gè)參數(shù)
browseVignettes("maftools")#可查看怎么自定義瀑布圖的幫助文檔html
*生信技能樹課程筆記
