237- 和刘小泽一起跟着官网学maftools
刘小泽写于2021.3.10 官网在：http://bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html
0 简单了解MAF文件
关于MAF文件格式
TCGA测了33种癌症，每个GDC project，每个检测变异的流程都会得到一个MAF。
MAF全称是：Mutation Annotation Format (MAF) ，是tab分隔的文件，用于记录VCF中的 somatic 突变信息。它有两种类型：protected and somatic (or open-access)。
*protected.maf：一般VCF文件中会记录多个转录本变异信息，但是这个maf文件只会记录其中影响最大的那个*somatic.maf：又称Masked Somatic Mutation ，它去除了低质量以及潜在的germline变异，可以公开获取。
如何得到MAF文件
对于VCF文件或者一般的tab分隔的文件，比较简单的办法是利用： vcf2maf ，另外GATK也支持 funcotator可以生成MAF文件
如果使用ANNOVAR注释，可以利用
annovarToMaf进行转换
MAF需要的列
如果看官方说明文档：https://docs.gdc.cancer.gov/Data/File_Formats/MAF_Format/，会发现有上百列的信息，从最基本的基因组坐标，到更复杂的cosmic变异注释。
必须包含的列：Hugo_Symbol, Chromosome, Start_Position, End_Position, Reference_Allele, Tumor_Seq_Allele2, Variant_Classification, Variant_Type and Tumor_Sample_Barcode
推荐添加的其他列： VAF (Variant Allele Frequecy) , amino acid change information.
其中，有这么几列需要看下什么意思：
Reference_Allele： The plus strand reference allele at this position. Includes the deleted sequence for a deletion or “-” for an insertion
Variant_Classification： Translational effect of variant allele
Variant_Type： Type of mutation. TNP (tri-nucleotide polymorphism) is analogous to DNP (di-nucleotide polymorphism) but for three consecutive nucleotides. ONP (oligo-nucleotide polymorphism) is analogous to TNP but for consecutive runs of four or more (SNP, DNP, TNP, ONP, INS, DEL, or Consolidated)
Tumor_Sample_Barcode： Aliquot barcode for the tumor sample
下面👇将以TCGA LAML cohort 为例，来看看maftools的使用
1 准备工作
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("maftools")
这个示意图解释的特别好： 把数据读取成maf对象，然后就可以画图或者分析
输入文件
（必须）MAF文件：可以是gz压缩的
（可选但建议）临床信息，并且与MAF中的sample对应
（可选）拷贝数信息：可以是GISTIC的结果或者自定义的tab文件（包含样本名称、基因名称以及拷贝数变异情况Amp或者Del）
2 读取文件
library(maftools)
suppressMessages(library(dplyr))
#（必须）TCGA LAML 的 MAF 文件
laml.maf = system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools')
#（可选）临床信息
laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools')
laml = read.maf(maf = laml.maf, clinicalData = laml.clin)
## -Reading
## -Validating
## -Silent variants: 475
## -Summarizing
## -Processing clinical data
## -Finished in 0.290s elapsed (0.243s cpu)
看看读取的对象中包含什么
laml
## An object of class MAF
## ID summary Mean Median
## 1: NCBI_Build 37 NA NA
## 2: Center genome.wustl.edu NA NA
## 3: Samples 193 NA NA
## 4: nGenes 1241 NA NA
## 5: Frame_Shift_Del 52 0.269 0
## 6: Frame_Shift_Ins 91 0.472 0
## 7: In_Frame_Del 10 0.052 0
## 8: In_Frame_Ins 42 0.218 0
## 9: Missense_Mutation 1342 6.953 7
## 10: Nonsense_Mutation 103 0.534 0
## 11: Splice_Site 92 0.477 0
## 12: total 1732 8.974 9
#Shows sample summry.
getSampleSummary(laml)
## Tumor_Sample_Barcode Frame_Shift_Del
## 1: TCGA-AB-3009 0
## 2: TCGA-AB-2807 1
## 3: TCGA-AB-2959 0
## 4: TCGA-AB-3002 0
## 5: TCGA-AB-2849 0
## ---
## 189: TCGA-AB-2942 0
## 190: TCGA-AB-2946 0
## 191: TCGA-AB-2954 0
## 192: TCGA-AB-2982 0
## 193: TCGA-AB-2903 0
## Frame_Shift_Ins In_Frame_Del In_Frame_Ins
## 1: 5 0 1
## 2: 0 1 0
## 3: 0 0 0
## 4: 0 0 0
## 5: 1 0 0
## ---
## 189: 0 0 1
## 190: 0 0 0
## 191: 0 0 0
## 192: 0 0 0
## 193: 0 0 0
## Missense_Mutation Nonsense_Mutation Splice_Site
## 1: 25 2 1
## 2: 16 3 4
## 3: 22 0 1
## 4: 15 1 5
## 5: 16 1 2
## ---
## 189: 0 0 0
## 190: 1 0 0
## 191: 1 0 0
## 192: 1 0 0
## 193: 0 0 0
## total
## 1: 34
## 2: 25
## 3: 23
## 4: 21
## 5: 20
## ---
## 189: 1
## 190: 1
## 191: 1
## 192: 1
## 193: 0
#Shows gene summary.
getGeneSummary(laml)
## Hugo_Symbol Frame_Shift_Del Frame_Shift_Ins
## 1: FLT3 0 0
## 2: DNMT3A 4 0
## 3: NPM1 0 33
## 4: IDH2 0 0
## 5: IDH1 0 0
## ---
## 1237: ZNF689 0 0
## 1238: ZNF75D 0 0
## 1239: ZNF827 1 0
## 1240: ZNF99 0 0
## 1241: ZPBP 0 0
## In_Frame_Del In_Frame_Ins Missense_Mutation
## 1: 1 33 15
## 2: 0 0 39
## 3: 0 0 1
## 4: 0 0 20
## 5: 0 0 18
## ---
## 1237: 0 0 1
## 1238: 0 0 1
## 1239: 0 0 0
## 1240: 0 0 1
## 1241: 0 0 1
## Nonsense_Mutation Splice_Site total MutatedSamples
## 1: 0 3 52 52
## 2: 5 6 54 48
## 3: 0 0 34 33
## 4: 0 0 20 20
## 5: 0 0 18 18
## ---
## 1237: 0 0 1 1
## 1238: 0 0 1 1
## 1239: 0 0 1 1
## 1240: 0 0 1 1
## 1241: 0 0 1 1
## AlteredSamples
## 1: 52
## 2: 48
## 3: 33
## 4: 20
## 5: 18
## ---
## 1237: 1
## 1238: 1
## 1239: 1
## 1240: 1
## 1241: 1
#shows clinical data associated with samples
getClinicalData(laml)
## Tumor_Sample_Barcode FAB_classification
## 1: TCGA-AB-2802 M4
## 2: TCGA-AB-2803 M3
## 3: TCGA-AB-2804 M3
## 4: TCGA-AB-2805 M0
## 5: TCGA-AB-2806 M1
## ---
## 189: TCGA-AB-3007 M3
## 190: TCGA-AB-3008 M1
## 191: TCGA-AB-3009 M4
## 192: TCGA-AB-3011 M1
## 193: TCGA-AB-3012 M3
## days_to_last_followup Overall_Survival_Status
## 1: 365 1
## 2: 792 1
## 3: 2557 0
## 4: 577 1
## 5: 945 1
## ---
## 189: 1581 0
## 190: 822 1
## 191: 577 1
## 192: 1885 0
## 193: 1887 0
#Shows all fields in MAF
getFields(laml)
## [1] "Hugo_Symbol" "Entrez_Gene_Id"
## [3] "Center" "NCBI_Build"
## [5] "Chromosome" "Start_Position"
## [7] "End_Position" "Strand"
## [9] "Variant_Classification" "Variant_Type"
## [11] "Reference_Allele" "Tumor_Seq_Allele1"
## [13] "Tumor_Seq_Allele2" "Tumor_Sample_Barcode"
## [15] "Protein_Change" "i_TumorVAF_WU"
## [17] "i_transcript_name"
#Writes maf summary to an output file with basename laml.
write.mafSummary(maf = laml, basename = 'laml')
3 可视化
3.1 MAF概览图 | Summary
plotmafSummary(maf = laml, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
# 下面看看它们各自利用的数据
# 图A
apply(laml@variant.classification.summary[,2:8], 2, sum)
## Frame_Shift_Del Frame_Shift_Ins In_Frame_Del
## 52 91 10
## In_Frame_Ins Missense_Mutation Nonsense_Mutation
## 42 1342 103
## Splice_Site
## 92
# 图B
apply(laml@variant.type.summary[,2:4], 2, sum)
## DEL INS SNP
## 64 141 1527
# 图D
head(laml@variants.per.sample)
## Tumor_Sample_Barcode Variants
## 1: TCGA-AB-3009 34
## 2: TCGA-AB-2807 25
## 3: TCGA-AB-2959 23
## 4: TCGA-AB-3002 21
## 5: TCGA-AB-2923 20
## 6: TCGA-AB-2972 20
# 图E和图A类似，只是把barplot换成了boxplot
# 图F
laml@gene.summary %>% arrange(desc(total))
## Hugo_Symbol Frame_Shift_Del Frame_Shift_Ins
## 1: DNMT3A 4 0
## 2: FLT3 0 0
## 3: NPM1 0 33
## 4: TET2 10 4
## 5: IDH2 0 0
## ---
## 1237: ZNF689 0 0
## 1238: ZNF75D 0 0
## 1239: ZNF827 1 0
## 1240: ZNF99 0 0
## 1241: ZPBP 0 0
## In_Frame_Del In_Frame_Ins Missense_Mutation
## 1: 0 0 39
## 2: 1 33 15
## 3: 0 0 1
## 4: 0 0 4
## 5: 0 0 20
## ---
## 1237: 0 0 1
## 1238: 0 0 1
## 1239: 0 0 0
## 1240: 0 0 1
## 1241: 0 0 1
## Nonsense_Mutation Splice_Site total MutatedSamples
## 1: 5 6 54 48
## 2: 0 3 52 52
## 3: 0 0 34 33
## 4: 8 1 27 17
## 5: 0 0 20 20
## ---
## 1237: 0 0 1 1
## 1238: 0 0 1 1
## 1239: 0 0 1 1
## 1240: 0 0 1 1
## 1241: 0 0 1 1
## AlteredSamples
## 1: 48
## 2: 52
## 3: 33
## 4: 17
## 5: 20
## ---
## 1237: 1
## 1238: 1
## 1239: 1
## 1240: 1
## 1241: 1
至于图3：
# https://rdrr.io/bioc/maftools/src/R/titv.R
maf=laml
mafDat = subsetMaf(maf = maf, query = "Variant_Type == 'SNP'", mafObj = FALSE, includeSyn = T, fields = "Variant_Type")
mafDat$con = paste(mafDat[,Reference_Allele], mafDat[,Tumor_Seq_Allele2], sep = '>')
conv = c("T>C", "T>C", "C>T", "C>T", "T>A", "T>A", "T>G", "T>G", "C>A", "C>A", "C>G", "C>G")
names(conv) = c('A>G', 'T>C', 'C>T', 'G>A', 'A>T', 'T>A', 'A>C', 'T>G', 'C>A', 'G>T', 'C>G', 'G>C')
mafDat = mafDat[con %in% names(conv)]
mafDat$con.class = conv[as.character(mafDat$con)]
table(mafDat$con.class )
##
## C>A C>G C>T T>A T>C T>G
## 235 139 1235 89 236 68
3.2 突变全景图/瀑布图 | Oncoplots
#oncoplot for top ten mutated genes.
oncoplot(maf = laml, top = 10)
- 默认绘制top10基因，当然可以设置参数
genes只展示某几个感兴趣的基因 - 注意其中有一个legend叫做：”Multi Hit“，它指的是某个基因在一个样本中突变次数超过1
给图例换个颜色
#One can use any colors, here in this example color palette from RColorBrewer package is used
vc_cols = RColorBrewer::brewer.pal(n = 8, name = 'Paired')
names(vc_cols) = c(
'Frame_Shift_Del',
'Missense_Mutation',
'Nonsense_Mutation',
'Multi_Hit',
'Frame_Shift_Ins',
'In_Frame_Ins',
'Splice_Site',
'In_Frame_Del'
)
print(vc_cols)
## Frame_Shift_Del Missense_Mutation Nonsense_Mutation
## "#A6CEE3" "#1F78B4" "#B2DF8A"
## Multi_Hit Frame_Shift_Ins In_Frame_Ins
## "#33A02C" "#FB9A99" "#E31A1C"
## Splice_Site In_Frame_Del
## "#FDBF6F" "#FF7F00"
oncoplot(maf = laml, colors = vc_cols, top = 10)
如果要包含拷贝数的数据
使用GISTIC结果
GISTIC会得到许多结果文件，我们需要的是：all_lesions.conf_XX.txt, amp_genes.conf_XX.txt, del_genes.conf_XX.txt, scores.gistic，其中XX指的是可信度（confidence level）
#GISTIC results LAML
all.lesions =
system.file("extdata", "all_lesions.conf_99.txt", package = "maftools")
amp.genes =
system.file("extdata", "amp_genes.conf_99.txt", package = "maftools")
del.genes =
system.file("extdata", "del_genes.conf_99.txt", package = "maftools")
scores.gis =
system.file("extdata", "scores.gistic", package = "maftools")
#Read GISTIC results along with MAF
laml.plus.gistic = read.maf(
maf = laml.maf,
gisticAllLesionsFile = all.lesions,
gisticAmpGenesFile = amp.genes,
gisticDelGenesFile = del.genes,
gisticScoresFile = scores.gis,
isTCGA = TRUE,
verbose = FALSE,
clinicalData = laml.clin
)
oncoplot(maf = laml.plus.gistic, top = 10)
使用自定义的结果
比如这里给DNMT3A基因创建一个20个样本的拷贝数变异数据
set.seed(seed = 1024)
barcodes = as.character(getSampleSummary(x = laml)[,Tumor_Sample_Barcode])
#Random 20 samples
dummy.samples = sample(x = barcodes,
size = 20,
replace = FALSE)
#Genarate random CN status for above samples
cn.status = sample(
x = c('ShallowAmp', 'DeepDel', 'Del', 'Amp'),
size = length(dummy.samples),
replace = TRUE
)
custom.cn.data = data.frame(
Gene = "DNMT3A",
Sample_name = dummy.samples,
CN = cn.status,
stringsAsFactors = FALSE
)
head(custom.cn.data)
## Gene Sample_name CN
## 1 DNMT3A TCGA-AB-2898 ShallowAmp
## 2 DNMT3A TCGA-AB-2879 Del
## 3 DNMT3A TCGA-AB-2920 Amp
## 4 DNMT3A TCGA-AB-2866 Del
## 5 DNMT3A TCGA-AB-2892 Del
## 6 DNMT3A TCGA-AB-2863 ShallowAmp
laml.plus.cn = read.maf(maf = laml.maf,
cnTable = custom.cn.data,
verbose = FALSE)
oncoplot(maf = laml.plus.cn, top = 5)
添加左右两侧的barplot
可以使用leftBarData, rightBarData, topBarData
#Selected AML driver genes
aml_genes = c("TP53", "WT1", "PHF6", "DNMT3A", "DNMT3B", "TET1", "TET2", "IDH1", "IDH2", "FLT3", "KIT", "KRAS", "NRAS", "RUNX1", "CEBPA", "ASXL1", "EZH2", "KDM6A")
#Variant allele frequcnies (Left bar plot)
aml_genes_vaf = subsetMaf(maf = laml, genes = aml_genes, fields = "i_TumorVAF_WU", mafObj = FALSE)[,mean(i_TumorVAF_WU, na.rm = TRUE), Hugo_Symbol]
colnames(aml_genes_vaf)[2] = "VAF"
head(aml_genes_vaf)
## Hugo_Symbol VAF
## 1: ASXL1 37.11250
## 2: CEBPA 22.00235
## 3: DNMT3A 43.51556
## 4: DNMT3B 37.14000
## 5: EZH2 68.88500
## 6: FLT3 34.60294
#MutSig results (Right bar plot)
laml.mutsig = system.file("extdata", "LAML_sig_genes.txt.gz", package = "maftools")
laml.mutsig = data.table::fread(input = laml.mutsig)[,.(gene, q)]
laml.mutsig[,q := -log10(q)] #transoform to log10
head(laml.mutsig)
## gene q
## 1: FLT3 12.64176
## 2: DNMT3A 12.64176
## 3: NPM1 12.64176
## 4: IDH2 12.64176
## 5: IDH1 12.64176
## 6: TET2 12.64176
# 如果这里画图报错，说找不到参数，那么可能就需要更新到R 版本4以上的maftools包了
oncoplot(
maf = laml,
genes = aml_genes,
leftBarData = aml_genes_vaf,
leftBarLims = c(0, 100),
rightBarData = laml.mutsig,
rightBarLims = c(0, 20)
)
添加临床信息
getClinicalData(x = laml)
## Tumor_Sample_Barcode FAB_classification
## 1: TCGA-AB-2802 M4
## 2: TCGA-AB-2803 M3
## 3: TCGA-AB-2804 M3
## 4: TCGA-AB-2805 M0
## 5: TCGA-AB-2806 M1
## ---
## 189: TCGA-AB-3007 M3
## 190: TCGA-AB-3008 M1
## 191: TCGA-AB-3009 M4
## 192: TCGA-AB-3011 M1
## 193: TCGA-AB-3012 M3
## days_to_last_followup Overall_Survival_Status
## 1: 365 1
## 2: 792 1
## 3: 2557 0
## 4: 577 1
## 5: 945 1
## ---
## 189: 1581 0
## 190: 822 1
## 191: 577 1
## 192: 1885 0
## 193: 1887 0
oncoplot(maf = laml, genes = aml_genes, clinicalFeatures = 'FAB_classification')
# 还可以自定义颜色
fabcolors = RColorBrewer::brewer.pal(n = 8,name = 'Spectral')
names(fabcolors) = c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7")
fabcolors = list(FAB_classification = fabcolors)
print(fabcolors)
## $FAB_classification
## M0 M1 M2 M3 M4
## "#D53E4F" "#F46D43" "#FDAE61" "#FEE08B" "#E6F598"
## M5 M6 M7
## "#ABDDA4" "#66C2A5" "#3288BD"
oncoplot(
maf = laml, genes = aml_genes,
clinicalFeatures = 'FAB_classification',
sortByAnnotation = TRUE,
annotationColor = fabcolors
)
突出显示某些样本
比如：突出显示alt allele is C的样本
# 看看有哪些列
getFields(x = laml)
## [1] "Hugo_Symbol" "Entrez_Gene_Id"
## [3] "Center" "NCBI_Build"
## [5] "Chromosome" "Start_Position"
## [7] "End_Position" "Strand"
## [9] "Variant_Classification" "Variant_Type"
## [11] "Reference_Allele" "Tumor_Seq_Allele1"
## [13] "Tumor_Seq_Allele2" "Tumor_Sample_Barcode"
## [15] "Protein_Change" "i_TumorVAF_WU"
## [17] "i_transcript_name"
# 然后使用参数additionalFeature
oncoplot(maf = laml, genes = aml_genes,
additionalFeature = c("Tumor_Seq_Allele2", "C"))
各种组合
oncoplot(
maf = laml.plus.gistic,
draw_titv = TRUE,
clinicalFeatures = c('FAB_classification', 'Overall_Survival_Status'),
sortByAnnotation = TRUE,
additionalFeature = c("Tumor_Seq_Allele2", "C"),
leftBarData = aml_genes_vaf,
leftBarLims = c(0, 100),
rightBarData = laml.mutsig[,.(gene, q)],
)
3.3 转换和颠换 | Transition and Transversions
名词解释：http://www.mun.ca/biology/scarr/Transitions_vs_Transversions.html
使用titv函数，画出箱线图，显示六种不同转换的总体分布,还有堆积条形图显示每个样本中的转换比例。
laml.titv = titv(maf = laml, plot = FALSE, useSyn = TRUE)
plotTiTv(res = laml.titv)
3.4 氨基酸突变热点 | Lollipop plots
许多癌基因，某些位点比其他位点更容易突变，变异频率更高，这些位置被认为是突变热点(hot-spots)。
使用函数lollipopPlot，需要我们在maf中给出氨基酸改变的信息。它默认寻找名字为"AAChange"的列，如果找不到，会给出一个warning信息，并且列出所有的列名。然后通过参数AACol进行指定。另外，它默认使用最长的gene isoform
lollipopPlot(maf = laml, gene = 'DNMT3A', AACol = 'Protein_Change', showMutationRate = TRUE)
## HGNC refseq.ID protein.ID aa.length
## 1: DNMT3A NM_175629 NP_783328 912
## 2: DNMT3A NM_022552 NP_072046 912
## 3: DNMT3A NM_153759 NP_715640 723
可以只画某个位点（如果设置labelPos=='all'，就会画所有的位点）
laml@data$Protein_Change[grepl('882',laml@data$Protein_Change)]
## [1] "p.R882C" "p.R882H" "p.R882C" "p.R882H" "p.R882P"
## [6] "p.R882H" "p.R882H" "p.R882H" "p.R882H" "p.R882H"
## [11] "p.R882H" "p.R882H" "p.R882H" "p.R882H" "p.R882C"
## [16] "p.R882C" "p.R882H" "p.R882C" "p.R882H" "p.R882H"
## [21] "p.R882H" "p.R882H" "p.R882H" "p.R882C" "p.R882H"
## [26] "p.R882C" "p.R882H"
lollipopPlot(maf = laml, gene = 'DNMT3A', showDomainLabel = FALSE, labelPos = 882)
## HGNC refseq.ID protein.ID aa.length
## 1: DNMT3A NM_175629 NP_783328 912
## 2: DNMT3A NM_022552 NP_072046 912
## 3: DNMT3A NM_153759 NP_715640 723
3.5 体细胞超突变 | Rainfall plots (RP)
名词解释：Somatic hypermutation (or SHM) 体细胞超突变, 指让基因组以背景突变率的至少十万倍到百万倍的高频率发生突变的细胞机制，主要形式是点突变（单碱基取代），偶见碱基插入或缺失。免疫系统的B细胞和多种癌细胞都可以使用此机制。
brca <- system.file("extdata", "brca.maf.gz", package = "maftools")
brca = read.maf(maf = brca, verbose = FALSE)
# If detectChangePoints is set to TRUE, rainfall plot also highlights regions where potential changes in inter-event distances are located.
rainfallPlot(maf = brca, detectChangePoints = TRUE, pointSize = 0.4)
## Chromosome Start_Position End_Position nMuts
## 1: 8 98129348 98133560 7
## 2: 8 98398549 98403536 9
## 3: 8 98453076 98456466 9
## 4: 8 124090377 124096810 22
## 5: 12 97436055 97439705 7
## 6: 17 29332072 29336153 8
## Avg_intermutation_dist Size Tumor_Sample_Barcode C>G
## 1: 702.0000 4212 TCGA-A8-A08B 4
## 2: 623.3750 4987 TCGA-A8-A08B 1
## 3: 423.7500 3390 TCGA-A8-A08B 1
## 4: 306.3333 6433 TCGA-A8-A08B 1
## 5: 608.3333 3650 TCGA-A8-A08B 4
## 6: 583.0000 4081 TCGA-A8-A08B 4
## C>T
## 1: 3
## 2: 8
## 3: 8
## 4: 21
## 5: 3
## 6: 4
图片会输出Kataegis的具体信息，并在图中用黑色箭头进行标记。
“Kataegis”are defined as those genomic segments containing six or more consecutive mutations with an average inter-mutation distance of less than or equal to 1,00 bp 即：包含至少6个连续突变的基因组区段，并且平均突变间距离不大于100bp
关于这个图的更多解释，有人写了篇paper：https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-017-1679-8
The plot is mostly used for illustrating the distribution of somatic cancer mutations along a reference genome, typically aiming to identify mutation hotspots. The x-coordinate shows the genomic position of the mutation, while the y-coordinate represents the base pair distance between consecutive mutations on a logarithmic scale.
3.6 样本肿瘤突变负荷 | TMB
TMB在许多癌症与预后相关，maftools也内置了33个TCGA癌症队列的突变数据
# 看源代码https://github.com/PoisonAlien/maftools/blob/master/R/tcgacompare.R，它使用的是mc3.v0.2.8.PUBLIC.maf.gz（https://api.gdc.cancer.gov/data/1c8cfe5f-e52d-41ba-94da-f15ea1337efc, 35.8 MB）
# TCGA的样本
tcga.cohort = system.file('extdata', 'tcga_cohort.txt.gz', package = 'maftools')
tcga.cohort = data.table::fread(file = tcga.cohort, sep = '\t', stringsAsFactors = FALSE)
head(tcga.cohort)
## cohort Tumor_Sample_Barcode total
## 1: ACC TCGA-OR-A5KB-01A-11D-A30A-10 1803
## 2: ACC TCGA-PK-A5HB-01A-11D-A29I-10 962
## 3: ACC TCGA-OR-A5JA-01A-11D-A29I-10 464
## 4: ACC TCGA-OR-A5LJ-01A-11D-A29I-10 395
## 5: ACC TCGA-OR-A5J5-01A-11D-A29I-10 349
## 6: ACC TCGA-OR-A5JB-01A-11D-A29I-10 277
## site
## 1: Adrenal Gland Carcinoma
## 2: Adrenal Gland Carcinoma
## 3: Adrenal Gland Carcinoma
## 4: Adrenal Gland Carcinoma
## 5: Adrenal Gland Carcinoma
## 6: Adrenal Gland Carcinoma
# 我们的样本
dim(laml@variants.per.sample)
## [1] 192 2
head(laml@variants.per.sample)
## Tumor_Sample_Barcode Variants
## 1: TCGA-AB-3009 34
## 2: TCGA-AB-2807 25
## 3: TCGA-AB-2959 23
## 4: TCGA-AB-3002 21
## 5: TCGA-AB-2923 20
## 6: TCGA-AB-2972 20
laml.mutload = tcgaCompare(maf = laml, cohortName = 'Example-LAML', logscale = TRUE, capture_size = 50)
3.7 变异等位基因频率 | Variant Allele Frequencies (VAF)
可以看：https://byteofbio.com/archives/7.html
VAF的全称是Variant Allele Frequency（变异等位基因频率）或Variant Allele Fraction（变异等位基因分数），指的是在基因组某个位点支持alternate/mutant allele的reads覆盖深度占这个位点总reads覆盖深度的比例。
在VCF文件中，DP表示Total Depth，AD表示Allele Depth
这里的VAF=AD/DP
plotVafhelps to estimate clonal status of top mutated genes (clonal genes usually have mean allele frequency around ~50% assuming pure sample)
plotVaf(maf = laml, vafCol = 'i_TumorVAF_WU')
4 处理拷贝数数据
还是利用GISTIC生成的4个文件：all_lesions.conf_XX.txt, amp_genes.conf_XX.txt, del_genes.conf_XX.txt, scores.gistic 从Firehose也能下载，但是版本就会比较老（2016_01_18），而且用的是hg19，而TCGA目前都是基于hg38
如果要从firehose下载
以TCGA-LAML为例： 链接版： http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/LAML/20160128/
界面版：
http://firebrowse.org/?cohort=LAML&download_dialog=true
然后这么操作：
不过还需要注意：基于hg19的firehose的基因组坐标，与官方TCGA的hg38坐标不一致，所以如果CNV数据从firehose下载，那么临床信息与MAF文件也需要从firehose下载
读取数据
# R包自带
all.lesions <- system.file("extdata", "all_lesions.conf_99.txt", package = "maftools")
amp.genes <- system.file("extdata", "amp_genes.conf_99.txt", package = "maftools")
del.genes <- system.file("extdata", "del_genes.conf_99.txt", package = "maftools")
scores.gis <- system.file("extdata", "scores.gistic", package = "maftools")
laml.gistic = readGistic(gisticAllLesionsFile = all.lesions,
gisticAmpGenesFile = amp.genes,
gisticDelGenesFile = del.genes,
gisticScoresFile = scores.gis,
isTCGA = TRUE)
## -Processing Gistic files..
## --Processing amp_genes.conf_99.txt
## --Processing del_genes.conf_99.txt
## --Processing scores.gistic
## --Summarizing by samples
laml.gistic
## An object of class GISTIC
## ID summary
## 1: Samples 191
## 2: nGenes 2622
## 3: cytoBands 16
## 4: Amp 388
## 5: Del 26481
## 6: total 26869
# 自己firehose下载
my.all.lesions <- 'gdac.broadinstitute.org_LAML-TB.CopyNumber_Gistic2.Level_4.2016012800.0.0/all_lesions.conf_99.txt'
my.amp.genes <- 'gdac.broadinstitute.org_LAML-TB.CopyNumber_Gistic2.Level_4.2016012800.0.0/amp_genes.conf_99.txt'
my.del.genes <- 'gdac.broadinstitute.org_LAML-TB.CopyNumber_Gistic2.Level_4.2016012800.0.0/del_genes.conf_99.txt'
my.scores.gis <- 'gdac.broadinstitute.org_LAML-TB.CopyNumber_Gistic2.Level_4.2016012800.0.0/scores.gistic'
my.laml.gistic = readGistic(gisticAllLesionsFile = my.all.lesions,
gisticAmpGenesFile = my.amp.genes,
gisticDelGenesFile = my.del.genes,
gisticScoresFile = my.scores.gis,
isTCGA = TRUE)
可视化
# genome plot
gisticChromPlot(gistic = laml.gistic, markBands = "all")
# Bubble plot
gisticBubblePlot(gistic = laml.gistic)
# oncoplot
# 图中显示：7q deletions 在M4亚型中是缺乏的，但在其他几个亚型中存在
laml.maf = system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools')
laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools')
laml = read.maf(maf = laml.maf, clinicalData = laml.clin)
## -Reading
## -Validating
## -Silent variants: 475
## -Summarizing
## -Processing clinical data
## -Finished in 0.256s elapsed (0.231s cpu)
gisticOncoPlot(gistic = laml.gistic,
clinicalData = getClinicalData(x = laml),
clinicalFeatures = 'FAB_classification',
sortByAnnotation = TRUE, top = 10)
除了输出GISTIC的结果，还可以输出DNAcopy的segment结果
tcga.ab.009.seg <- system.file("extdata", "TCGA.AB.3009.hg19.seg.txt", package = "maftools")
plotCBSsegments(cbsFile = tcga.ab.009.seg)
## NULL
5 可以做的分析
5.1 检测互斥或共存的突变 | somaticInteractions
它是选队列中突变最多的top基因（设定参数top=10），或者指定基因（设定参数genes=c()），进行pair-wise Fisher’s Exact test
检测共突变的意义是区别于单基因突变，治疗方法和后期用药都会存在不同。结果中绿色是倾向共存的突变基因；黄色是倾向互斥，颜色深浅表示显著性
#exclusive/co-occurance event analysis on top 10 mutated genes.
somaticInteractions(maf = laml, top = 25, pvalue = c(0.05, 0.1))
## gene1 gene2 pValue oddsRatio 00 11 01 10
## 1: ASXL1 RUNX1 0.0001541586 55.215541 176 4 12 1
## 2: IDH2 RUNX1 0.0002809928 9.590877 164 7 9 13
## 3: IDH2 ASXL1 0.0004030636 41.077327 172 4 1 16
## 4: FLT3 NPM1 0.0009929836 3.763161 125 17 16 35
## 5: SMC3 DNMT3A 0.0010451985 20.177713 144 6 42 1
## ---
## 296: PLCE1 ASXL1 1.0000000000 0.000000 184 0 5 4
## 297: RAD21 FAM5C 1.0000000000 0.000000 183 0 5 5
## 298: PLCE1 FAM5C 1.0000000000 0.000000 184 0 5 4
## 299: PLCE1 RAD21 1.0000000000 0.000000 184 0 5 4
## 300: EZH2 PLCE1 1.0000000000 0.000000 186 0 4 3
## Event pair event_ratio
## 1: Co_Occurence ASXL1, RUNX1 4/13
## 2: Co_Occurence IDH2, RUNX1 7/22
## 3: Co_Occurence ASXL1, IDH2 4/17
## 4: Co_Occurence FLT3, NPM1 17/51
## 5: Co_Occurence DNMT3A, SMC3 6/43
## ---
## 296: Mutually_Exclusive ASXL1, PLCE1 0/9
## 297: Mutually_Exclusive FAM5C, RAD21 0/10
## 298: Mutually_Exclusive FAM5C, PLCE1 0/9
## 299: Mutually_Exclusive PLCE1, RAD21 0/9
## 300: Mutually_Exclusive EZH2, PLCE1 0/7
5.2 预测癌症驱动基因 | oncodrive
函数oncodrive基于
oncodriveCLUST 算法，从MAF文件中预测驱动基因。它是根据：癌症驱动基因的突变一般富集在某些特定位点（即hot-spots），然后利用这些位点去鉴定基因。
laml.sig = oncodrive(maf = laml, AACol = 'Protein_Change', minMut = 5, pvalMethod = 'zscore')
head(laml.sig)
## Hugo_Symbol Frame_Shift_Del Frame_Shift_Ins
## 1: IDH1 0 0
## 2: IDH2 0 0
## 3: NPM1 0 33
## 4: NRAS 0 0
## 5: U2AF1 0 0
## 6: KIT 1 1
## In_Frame_Del In_Frame_Ins Missense_Mutation
## 1: 0 0 18
## 2: 0 0 20
## 3: 0 0 1
## 4: 0 0 15
## 5: 0 0 8
## 6: 0 1 7
## Nonsense_Mutation Splice_Site total MutatedSamples
## 1: 0 0 18 18
## 2: 0 0 20 20
## 3: 0 0 34 33
## 4: 0 0 15 15
## 5: 0 0 8 8
## 6: 0 0 10 8
## AlteredSamples clusters muts_in_clusters
## 1: 18 1 18
## 2: 20 2 20
## 3: 33 2 32
## 4: 15 2 15
## 5: 8 1 7
## 6: 8 2 9
## clusterScores protLen zscore pval
## 1: 1.0000000 414 5.546154 1.460110e-08
## 2: 1.0000000 452 5.546154 1.460110e-08
## 3: 0.9411765 294 5.093665 1.756034e-07
## 4: 0.9218951 189 4.945347 3.800413e-07
## 5: 0.8750000 240 4.584615 2.274114e-06
## 6: 0.8500000 976 4.392308 5.607691e-06
## fdr fract_muts_in_clusters
## 1: 1.022077e-07 1.0000000
## 2: 1.022077e-07 1.0000000
## 3: 8.194826e-07 0.9411765
## 4: 1.330144e-06 1.0000000
## 5: 6.367520e-06 0.8750000
## 6: 1.308461e-05 0.9000000
plotOncodrive(res = laml.sig, fdrCutOff = 0.1, useFraction = TRUE, labelSize = 0.5,bubbleSize=0.1)
图中点的大小对应基因内突变cluster的数量，x轴是这些cluster中突变数量或比例。上图中，IDH1的18个mutation都存在于1个cluster中，所以显示的基因名后面的括号中显示1
什么是mutation cluster？（来自https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4710516/） An important biological consequence of non-uniform distribution of mutations is that a disproportionately large number of changes can accumulate in a limited section of a genome, forming a mutation cluster
5.3 蛋白结构域注释 | pfamDomains
Pfam 数据库是蛋白结构域注释数据库，它使用基于隐马尔可夫模型的多序列比对方法对蛋白质序列进行家族分类。分为Pfam-A与Pfam-B两个库：
Pfam-A：包含经人工编辑，有完整注释的高质量记录条目
Pfam-B：计算机工具自动编辑，除Pfam-A外的条目，质量较低
这里，可以利用pfamDomains来看看氨基酸改变区域的结构域信息，帮助了解给定的癌症队列中，哪些结构域是经常被影响的。
下图中，圆点表示结构域，x轴是该结构域内存在的突变数量，y轴和圆点的大小是影响的基因数量
# 选择突变最多的10个pfam domain
laml.pfam = pfamDomains(maf = laml, AACol = 'Protein_Change', top = 10)
#Protein summary (Printing first 7 columns for display convenience)
laml.pfam$proteinSummary[,1:7, with = FALSE]
## HGNC AAPos Variant_Classification N total
## 1: DNMT3A 882 Missense_Mutation 27 54
## 2: IDH1 132 Missense_Mutation 18 18
## 3: IDH2 140 Missense_Mutation 17 20
## 4: FLT3 835 Missense_Mutation 14 52
## 5: FLT3 599 In_Frame_Ins 10 52
## ---
## 1512: ZNF646 875 Missense_Mutation 1 1
## 1513: ZNF687 554 Missense_Mutation 1 2
## 1514: ZNF687 363 Missense_Mutation 1 2
## 1515: ZNF75D 5 Missense_Mutation 1 1
## 1516: ZNF827 427 Frame_Shift_Del 1 1
## fraction DomainLabel
## 1: 0.5000000 AdoMet_MTases
## 2: 1.0000000 PTZ00435
## 3: 0.8500000 PTZ00435
## 4: 0.2692308 PKc_like
## 5: 0.1923077 PKc_like
## ---
## 1512: 1.0000000 <NA>
## 1513: 0.5000000 <NA>
## 1514: 0.5000000 <NA>
## 1515: 1.0000000 <NA>
## 1516: 1.0000000 <NA>
#Domain summary (Printing first 3 columns for display convenience)
laml.pfam$domainSummary[,1:3, with = FALSE]
## DomainLabel nMuts nGenes
## 1: PKc_like 55 5
## 2: PTZ00435 38 2
## 3: AdoMet_MTases 33 1
## 4: 7tm_1 24 24
## 5: COG5048 17 17
## ---
## 499: ribokinase 1 1
## 500: rim_protein 1 1
## 501: sigpep_I_bact 1 1
## 502: trp 1 1
## 503: zf-BED 1 1
5.4 生存分析 | Survival analysis
mafSurvival函数可以根据我们自己定义的基因突变情况进行分组，或者自定义样本进行分组，来绘制KM plot。需要输入数据包含 Tumor_Sample_Barcode (要和MAF文件的对应), binary event (1/0), time to event
利用给定基因进行分组，看看突变与不突变的生存差异
#Survival analysis based on grouping of DNMT3A mutation status
# 参数time表示随访时间
# 参数Status表示生存状态
mafSurvival(maf = laml, genes = 'DNMT3A', time = 'days_to_last_followup', Status = 'Overall_Survival_Status', isTCGA = TRUE)
## DNMT3A
## 48
## Group medianTime N
## 1: Mutant 245 45
## 2: WT 396 137
参数isTCGA在这里设定数据是否来自TCGA，看源代码会发现这么一句，表示如果是来自TCGA的数据，临床信息会将前12个字符提取出来，用于和MAF对应：
if (isTCGA) {
clinicalData$Tumor_Sample_Barcode = substr(x = clinicalData$Tumor_Sample_Barcode,
start = 1, stop = 12)
}
与预后较差有关基因集的预测
#Using top 20 mutated genes to identify a set of genes (of size 2) to predict poor prognostic groups
prog_geneset = survGroup(maf = laml, top = 20, geneSetSize = 2, time = "days_to_last_followup", Status = "Overall_Survival_Status", verbose = FALSE)
print(prog_geneset)
## Gene_combination P_value hr WT Mutant
## 1: FLT3_DNMT3A 0.00104 2.510 164 18
## 2: DNMT3A_SMC3 0.04880 2.220 176 6
## 3: DNMT3A_NPM1 0.07190 1.720 166 16
## 4: DNMT3A_TET2 0.19600 1.780 176 6
## 5: FLT3_TET2 0.20700 1.860 177 5
## 6: NPM1_IDH1 0.21900 0.495 176 6
## 7: DNMT3A_IDH1 0.29300 1.500 173 9
## 8: IDH2_RUNX1 0.31800 1.580 176 6
## 9: FLT3_NPM1 0.53600 1.210 165 17
## 10: DNMT3A_IDH2 0.68000 0.747 178 4
## 11: DNMT3A_NRAS 0.99200 0.986 178 4
其中p<0.05表示预后较差，可以选其中一个基因集来看看
mafSurvGroup(maf = laml, geneSet = c("DNMT3A", "FLT3"), time = "days_to_last_followup", Status = "Overall_Survival_Status")
## Group medianTime N
## 1: Mutant 242.5 18
## 2: WT 379.5 164
5.5 比较两个MAF文件 | mafComapre
不同癌症、同种癌症的不同时期之间，或者不同性别之间，它们的突变模式（mutation pattern）是不同的，我们可以比较两个不同的队列（maf文件）来检测这些差异的突变基因。反正只要能明确分组，就可以进行比较
例如：急性早幼粒细胞白血病(Acute Promyelocytic Leukemia)在复发期体现出PML、RARA基因的突变，但是早期没有
这里可以利用mafComapre 对两个不同时期进行fisher检验
#Primary APL MAF
primary.apl = system.file("extdata", "APL_primary.maf.gz", package = "maftools")
primary.apl = read.maf(maf = primary.apl)
## -Reading
## -Validating
## --Non MAF specific values in Variant_Classification column:
## ITD
## -Silent variants: 45
## -Summarizing
## -Processing clinical data
## --Missing clinical data
## -Finished in 0.088s elapsed (0.080s cpu)
#Relapse APL MAF
relapse.apl = system.file("extdata", "APL_relapse.maf.gz", package = "maftools")
relapse.apl = read.maf(maf = relapse.apl)
## -Reading
## -Validating
## --Non MAF specific values in Variant_Classification column:
## ITD
## -Silent variants: 19
## -Summarizing
## -Processing clinical data
## --Missing clinical data
## -Finished in 0.102s elapsed (0.076s cpu)
#Considering only genes which are mutated in at-least in 5 samples in one of the cohort to avoid bias due to genes mutated in single sample.
pt.vs.rt <- mafCompare(m1 = primary.apl, m2 = relapse.apl, m1Name = 'Primary', m2Name = 'Relapse', minMut = 5)
print(pt.vs.rt)
## $results
## Hugo_Symbol Primary Relapse pval or
## 1: PML 1 11 1.529935e-05 0.03537381
## 2: RARA 0 7 2.574810e-04 0.00000000
## 3: RUNX1 1 5 1.310500e-02 0.08740567
## 4: FLT3 26 4 1.812779e-02 3.56086275
## 5: ARID1B 5 8 2.758396e-02 0.26480490
## 6: WT1 20 14 2.229087e-01 0.60619329
## 7: KRAS 6 1 4.334067e-01 2.88486293
## 8: NRAS 15 4 4.353567e-01 1.85209500
## 9: ARID1A 7 4 7.457274e-01 0.80869223
## ci.up ci.low adjPval
## 1: 0.2552937 0.000806034 0.0001376942
## 2: 0.3006159 0.000000000 0.0011586643
## 3: 0.8076265 0.001813280 0.0393149868
## 4: 14.7701728 1.149009169 0.0407875250
## 5: 0.9698686 0.064804160 0.0496511201
## 6: 1.4223101 0.263440988 0.3343630535
## 7: 135.5393108 0.337679367 0.4897762916
## 8: 8.0373994 0.553883512 0.4897762916
## 9: 3.9297309 0.195710173 0.7457273717
##
## $SampleSummary
## Cohort SampleSize
## 1: Primary 124
## 2: Relapse 58
从表中可以看到PML和RARA都是在复发期的突变更多，可以作图再看看：
# 森林图
forestPlot(mafCompareRes = pt.vs.rt, pVal = 0.1, color = c('royalblue', 'maroon'), geneFontSize = 0.5)
# 还可以把这些基因在不同时期的突变画一个瀑布图
genes = c("PML", "RARA", "RUNX1", "ARID1B", "FLT3")
coOncoplot(m1 = primary.apl, m2 = relapse.apl, m1Name = 'PrimaryAPL', m2Name = 'RelapseAPL', genes = genes, removeNonMutated = TRUE)
# 或者条形图
coBarplot(m1 = primary.apl, m2 = relapse.apl, m1Name = "Primary", m2Name = "Relapse")
# 还能看看不同分组之间的突变位点和类型
lollipopPlot2(m1 = primary.apl, m2 = relapse.apl, gene = "PML", AACol1 = "amino_acid_change", AACol2 = "amino_acid_change", m1_name = "Primary", m2_name = "Relapse")
## HGNC refseq.ID protein.ID aa.length
## 1: PML NM_033238 NP_150241 882
## 2: PML NM_002675 NP_002666 633
## 3: PML NM_033249 NP_150252 585
## 4: PML NM_033247 NP_150250 435
## 5: PML NM_033239 NP_150242 829
## 6: PML NM_033250 NP_150253 781
## 7: PML NM_033240 NP_150243 611
## 8: PML NM_033244 NP_150247 560
## 9: PML NM_033246 NP_150249 423
## HGNC refseq.ID protein.ID aa.length
## 1: PML NM_033238 NP_150241 882
## 2: PML NM_002675 NP_002666 633
## 3: PML NM_033249 NP_150252 585
## 4: PML NM_033247 NP_150250 435
## 5: PML NM_033239 NP_150242 829
## 6: PML NM_033250 NP_150253 781
## 7: PML NM_033240 NP_150243 611
## 8: PML NM_033244 NP_150247 560
## 9: PML NM_033246 NP_150249 423
5.6 临床富集分析 | clinicalEnrichment
拿到临床信息中的某个特征，然后看看出现了哪些突变
这个函数clinicalEnrichment()采用Fisher精确检验，会进行两两比较（比如M1和M2，M1和M3），以及分组比较（比如M1和其他各组集合），结果分别存放在pairwise_comparision和groupwise_comparision
fab.ce = clinicalEnrichment(maf = laml, clinicalFeature = 'FAB_classification')
##
## M0 M1 M2 M3 M4 M5 M6 M7
## 19 44 44 21 39 19 3 3
# 两两比较结果
head(fab.ce$pairwise_comparision)
## Hugo_Symbol Feature_1 Feature_2 n_mutated_Feature1
## 1: TP53 M7 M1 3 of 3
## 2: TP53 M7 M3 3 of 3
## 3: TP53 M7 M5 3 of 3
## 4: TP53 M7 M4 3 of 3
## 5: TP53 M7 M2 3 of 3
## 6: DNMT3A M5 M3 10 of 19
## n_mutated_Feature2 fdr
## 1: 2 of 44 0.005844156
## 2: 0 of 21 0.005844156
## 3: 0 of 19 0.005844156
## 4: 2 of 39 0.005879791
## 5: 3 of 44 0.006660500
## 6: 1 of 21 0.028683522
#分组比较结果 Significant associations p-value < 0.05
fab.ce$groupwise_comparision[p_value < 0.05]
## Hugo_Symbol Group1 Group2 n_mutated_group1
## 1: IDH1 M1 Rest 11 of 44
## 2: TP53 M7 Rest 3 of 3
## 3: DNMT3A M5 Rest 10 of 19
## 4: CEBPA M2 Rest 7 of 44
## 5: RUNX1 M0 Rest 5 of 19
## 6: NPM1 M5 Rest 7 of 19
## 7: NPM1 M3 Rest 0 of 21
## 8: DNMT3A M3 Rest 1 of 21
## n_mutated_group2 p_value OR OR_low
## 1: 7 of 149 0.0002597371 6.670592 2.173829026
## 2: 12 of 190 0.0003857187 Inf 5.355415451
## 3: 38 of 174 0.0089427384 3.941207 1.333635173
## 4: 6 of 149 0.0117352110 4.463237 1.204699322
## 5: 11 of 174 0.0117436825 5.216902 1.243812880
## 6: 26 of 174 0.0248582372 3.293201 1.001404899
## 7: 33 of 172 0.0278630823 0.000000 0.000000000
## 8: 47 of 172 0.0294005111 0.133827 0.003146708
## OR_high fdr
## 1: 21.9607250 0.0308575
## 2: Inf 0.0308575
## 3: 11.8455979 0.3757978
## 4: 17.1341278 0.3757978
## 5: 19.4051505 0.3757978
## 6: 10.1210509 0.5880102
## 7: 0.8651972 0.5880102
## 8: 0.8848897 0.5880102
看到IDH1 突变主要发生在M1亚型中 ，同理DNMT3A在M5，TP53在M7。这些信息也能进行可视化
plotEnrichmentResults(enrich_res = fab.ce, pVal = 0.05, geneFontSize = 0.5, annoFontSize = 0.6)
5.7 药物基因互作 | drugInteractions
它可以查询DGIdb数据库(Drug-Gene Interaction database)，用来分析药物-基因作用以及基因的成药性（gene druggability）。不过这里只是用于研究，不能用于临床
DGIdb是一个药物-基因相互作用数据库，其中主要是癌基因，当然也不乏其他健康相关的基因（如心脏病、糖尿病、阿尔兹海默症）。DGIdb contains over 40,000 genes and 10,000 drugs involved in over 100,000 drug-gene interactions or belonging to one of 42 potentially druggable gene categories.
# 显示潜在可成药基因（少于5个基因就全部显示，多于5个就显示前5个）
dgi = drugInteractions(maf = laml, fontSize = 0.75)
# 提供基因名，返回互作的药物
dnmt3a.dgi = drugInteractions(genes = "DNMT3A", drugs = TRUE)
## Number of claimed drugs for given genes:
## Gene N
## 1: DNMT3A 7
#Printing selected columns.
dnmt3a.dgi[,.(Gene, interaction_types, drug_name, drug_claim_name)]
## Gene interaction_types drug_name drug_claim_name
## 1: DNMT3A N/A
## 2: DNMT3A DAUNORUBICIN Daunorubicin
## 3: DNMT3A DECITABINE Decitabine
## 4: DNMT3A IDARUBICIN IDARUBICIN
## 5: DNMT3A DECITABINE DECITABINE
## 6: DNMT3A inhibitor DECITABINE CHEMBL1201129
## 7: DNMT3A inhibitor AZACITIDINE CHEMBL1489
5.8 致癌信号通路 | OncogenicPathways
OncogenicPathways(maf = laml)
## Pathway N n_affected_genes fraction_affected
## 1: PI3K 29 1 0.03448276
## 2: NRF2 3 1 0.33333333
## 3: TP53 6 2 0.33333333
## 4: WNT 68 3 0.04411765
## 5: MYC 13 3 0.23076923
## 6: NOTCH 71 6 0.08450704
## 7: Hippo 38 7 0.18421053
## 8: RTK-RAS 85 18 0.21176471
## Mutated_samples Fraction_mutated_samples
## 1: 1 0.005181347
## 2: 1 0.005181347
## 3: 15 0.077720207
## 4: 4 0.020725389
## 5: 3 0.015544041
## 6: 8 0.041450777
## 7: 7 0.036269430
## 8: 97 0.502590674
# 例如：RTK-RAS这个pathway中所有基因突变情况
PlotOncogenicPathways(maf = laml, pathways = "RTK-RAS")
其中抑癌基因是红色，致癌基因是蓝色
5.9 肿瘤异质性和MATH分数 | inferHeterogeneity
肿瘤是普遍存在异质性的，例如包含多个克隆。这种异质性可以通过等位基因频率（ variant allele frequencies, VAF）的聚类进行推断。
关于等位基因频率（VAF）(来自：t.ly/nDvy)： 如果在某种群中一个等位基因的基因频率为20%，那么在种群的所有成员中，1/5的染色体带有那个等位基因，而其他4/5的染色体带有该等位基因的其他对应变种——可以是一种也可以是很多种 (来自：https://doi.org/10.1186/s13073-019-0664-4) The variant allele frequency (VAF; also known as variant allele fraction) is used to infer whether a variant comes from somatic cells or inherited from parents when a matched normal sample is not provided. A variant is potentially a germline mutation if the VAF is approximately 50% or 100%.
这个聚类操作会使用mclust，函数会寻找带有VAF信息的一列（当然可以用vafCol来指定）
library("mclust")
tcga.ab.2972.het = inferHeterogeneity(maf = laml, tsb = 'TCGA-AB-2972', vafCol = 'i_TumorVAF_WU')
# 各突变基因的VAF数据
print(tcga.ab.2972.het$clusterData)
## Hugo_Symbol Chromosome Start_Position End_Position
## 1: HMCN1 1 185958621 185958621
## 2: TUFT1 1 151546791 151546791
## 3: DNAH5 5 13894834 13894834
## 4: DOCK2 5 169127134 169127134
## 5: RNASEN 5 31409392 31409392
## 6: KIAA0240 6 42797687 42797687
## 7: C10orf118 10 115891901 115891901
## 8: FANCI 15 89847119 89847119
## 9: DNAH3 16 20994057 20994057
## 10: ZC3H18 16 88688690 88688690
## 11: ZNF43 19 21991712 21991712
## 12: MORC3 21 37711165 37711165
## 13: ATP1B4 23 119505057 119505057
## 14: STAG2 23 123197716 123197716
## 15: RIMS1 6 72806854 72806854
## 16: LARP4B 10 871204 871204
## 17: PTPN11 12 112926884 112926884
## 18: SFRS6 20 42088521 42088521
## 19: ASTL 2 96801149 96801149
## Tumor_Sample_Barcode t_vaf cluster MATH
## 1: TCGA-AB-2972 0.4858 2 11.62157
## 2: TCGA-AB-2972 0.4328 2 11.62157
## 3: TCGA-AB-2972 0.4473 2 11.62157
## 4: TCGA-AB-2972 0.4021 2 11.62157
## 5: TCGA-AB-2972 0.4422 2 11.62157
## 6: TCGA-AB-2972 0.4363 2 11.62157
## 7: TCGA-AB-2972 0.4843 2 11.62157
## 8: TCGA-AB-2972 0.4395 2 11.62157
## 9: TCGA-AB-2972 0.4715 2 11.62157
## 10: TCGA-AB-2972 0.4315 2 11.62157
## 11: TCGA-AB-2972 0.4891 2 11.62157
## 12: TCGA-AB-2972 0.4425 2 11.62157
## 13: TCGA-AB-2972 0.4300 2 11.62157
## 14: TCGA-AB-2972 0.4603 2 11.62157
## 15: TCGA-AB-2972 0.1985 1 11.62157
## 16: TCGA-AB-2972 0.2704 1 11.62157
## 17: TCGA-AB-2972 0.2516 1 11.62157
## 18: TCGA-AB-2972 0.2614 1 11.62157
## 19: TCGA-AB-2972 0.3695 outlier 11.62157
## MedianAbsoluteDeviation
## 1: 3.42
## 2: 3.42
## 3: 3.42
## 4: 3.42
## 5: 3.42
## 6: 3.42
## 7: 3.42
## 8: 3.42
## 9: 3.42
## 10: 3.42
## 11: 3.42
## 12: 3.42
## 13: 3.42
## 14: 3.42
## 15: 3.42
## 16: 3.42
## 17: 3.42
## 18: 3.42
## 19: 3.42
# 聚类结果
print(tcga.ab.2972.het$clusterMeans)
## Tumor_Sample_Barcode cluster meanVaf
## 1: TCGA-AB-2972 2 0.4496571
## 2: TCGA-AB-2972 1 0.2454750
## 3: TCGA-AB-2972 outlier 0.3695000
# 可视化
plotClusters(clusters = tcga.ab.2972.het)
上图很明显分成2群，说明可能存在2个克隆。一群的mean variant allele frequencies为45%左右（为主克隆），另一群是25%左右。
虽然这里可能通过数据分布来直观判断，但是没办法具体量化。因此还可以看：Mutant-Allele Tumor Heterogeneity (MATH)
MATH is a simple quantitative measure of intra-tumor heterogeneity, which calculates the width of the vaf distribution.（即量化VAF分布的离散程度） MATH值越高，预后越差 (计算方法来自：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5441963/) The MATH score is calculated as 100 × median absolute deviation (MAD)/median of the variant allele frequencies, and describes the ratio of the width of the data to the center of the distribution
考虑一个特殊情况：拷贝数变异
突变基因的拷贝数变异一般范围会影响到VAF的计算，可能会出现异常高/低，从而影响
拷贝数变异(Copy number variation, CNV)，是由基因组发生重排而导致的,一般指长度1Kb-3Mb 的基因组大片段的拷贝数增加或者减少（多数 CNV 低于 500kb），是基因组结构变异(Structural variation, SV) 的重要组成部分。CNV变异是人类最常见的变异形式。
CNV数据可以通过 DNACopy的Circular Binary Segmentation方法得到Copy Number Segment文件
suppressMessages(library(dplyr))
seg = system.file('extdata', 'TCGA.AB.3009.hg19.seg.txt', package = 'maftools')
read.table(seg) %>% head()
## V1 V2 V3 V4 V5
## 1 Sample Chromosome Start End Num_Probes
## 2 TCGA-AB-3009 1 61735 86341 15
## 3 TCGA-AB-3009 1 98588 3166089 717
## 4 TCGA-AB-3009 1 3169157 9308037 3956
## 5 TCGA-AB-3009 1 9315847 9404980 60
## 6 TCGA-AB-3009 1 9406979 12845863 1777
## V6
## 1 Segment_Mean
## 2 -0.683
## 3 0.3069
## 4 0.0039
## 5 0.2934
## 6 0.0104
tcga.ab.3009.het = inferHeterogeneity(maf = laml, tsb = 'TCGA-AB-3009', segFile = seg, vafCol = 'i_TumorVAF_WU')
## Hugo_Symbol Chromosome Start_Position End_Position
## 1: PHF6 23 133551224 133551224
## Tumor_Sample_Barcode t_vaf Segment_Start
## 1: TCGA-AB-3009 0.9349112 NA
## Segment_End Segment_Mean CN
## 1: NA NA NA
## Hugo_Symbol Chromosome Start_Position End_Position
## 1: NFKBIL2 8 145668658 145668658
## 2: NF1 17 29562981 29562981
## 3: SUZ12 17 30293198 30293198
## Tumor_Sample_Barcode t_vaf Segment_Start
## 1: TCGA-AB-3009 0.4415584 145232496
## 2: TCGA-AB-3009 0.8419000 29054355
## 3: TCGA-AB-3009 0.8958333 29054355
## Segment_End Segment_Mean CN cluster
## 1: 145760746 0.3976 2.634629 CN_altered
## 2: 30363868 -0.9157 1.060173 CN_altered
## 3: 30363868 -0.9157 1.060173 CN_altered
# 根据这个segment文件，函数就会将没有CNV的基因作为离群值，并且标记出有CNV突变的基因
plotClusters(clusters = tcga.ab.3009.het, genes = 'CN_altered', showCNvars = TRUE)
上图中，NF1和SUZ12基因就存在较高的VAF，这是因为出现了CNV（deletion）。为什么说是出现了deletion，是因为上面的结果中：segment mean这一列
来自: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/ The copy number variation (CNV) pipeline uses Affymetrix SNP 6.0 array data to identify genomic regions that are repeated and infer the copy number of these repeats. Diploid regions will have a segment mean of zero, amplified regions will have positive values, and deletions will have negative values.
5.10 突变特征/频谱 | Mutational Signatures
首先，什么是突变特征/频谱？
每种癌症在发展过程中，都会出现一些比较特有的碱基替换，可以认为是一种特征（signature），并且这种特征只是针对于somatic 的mutation，一般一个样本也就几百个somatic 的mutation。somatic突变在全基因组上发生的比例约百万分之一，如果是全基因组肿瘤基因测序， 可能会有3万个左右的somatic突变，如果是全外显子测序，是300个左右。
在http://www.nature.com/nature/journal/v500/n7463/full/nature12477.html这篇文章中，研究了4,938,362 mutations from 7,042 cancers样本；研究了30种癌症，发现21种不同的mutation signature
COSMIC的30个突变特征
来自（https://cancer.sanger.ac.uk/cosmic/signatures）: Somatic mutations are present in all cells of the human body and occur throughout life. Different mutational processes generate unique combinations of mutation types, termed “Mutational Signatures”.
sanger研究所科学家提出来了肿瘤somatic突变的signature概念 ，把96突变频谱的非负矩阵分解后的30个特征。R包deconstructSigs可以把自己的96突变频谱对应到cosmic数据库的30个突变特征，不同的特征有不同的生物学含义
它的30个signature的96突变频谱在：https://cancer.sanger.ac.uk/cancergenome/assets/signatures_probabilities.txt
cosmic=read.table('https://cancer.sanger.ac.uk/cancergenome/assets/signatures_probabilities.txt',
header = T,sep = '\t')[,1:33]
head(cosmic[,1:3])
tmp=cosmic[,2];substr(tmp,2,2) <- '.'
rownames(cosmic)=paste(gsub('>','',cosmic[,1]),
tmp)
96突变频谱是怎么回事？
看这个图：https://en.wikipedia.org/wiki/Mutational_signatures
正链上的C->A突变，对应负链上为G->T，因此某个位点的突变可以划分为6种：
C>A, 表示C>A和G>T两种 C>G, 表示C>G和G>C两种 C>T, 表示C>T和G>A两种 T>A,表示T>A和A>T两种 T>C,表示T>C和A>G两种 T>G,表示T>G和A>C两种
然后在突变位点的上下游各取一个碱基，组成3碱基，所以有4x6x4=96种模式，每种模式的频率分布就是突变频谱。它可以作为一个肿瘤样本的特征，从而进行样本间的比较。比如MutationalPatterns这个R包就可以做这样的分析
我们接下来使用maftools进行计算
#Requires BSgenome object
library(BSgenome.Hsapiens.UCSC.hg19, quietly = TRUE)
laml.tnm = trinucleotideMatrix(maf = laml, prefix = 'chr', add = TRUE, ref_genome = "BSgenome.Hsapiens.UCSC.hg19")
## -Extracting 5' and 3' adjacent bases
## -Extracting +/- 20bp around mutated bases for background C>T estimation
## -Estimating APOBEC enrichment scores
## --Performing one-way Fisher's test for APOBEC enrichment
## ---APOBEC related mutations are enriched in 3.315 % of samples (APOBEC enrichment score > 2 ; 6 of 181 samples)
## -Creating mutation matrix
## --matrix of dimension 188x96
这个计算过程主要分成两部分：
Estimates APOBEC enrichment scores
Prepares a mutational matrix for signature analysis.
APOBEC富集与非富集样本的差异分析
plotApobecDiff函数可以拿trinucleotideMatrix的APOBED富集分数，然后将样本分成两类：APOBEC enriched and non-APOBEC enriched。之后，鉴定两类之间差异基因
# 这里仅仅是代码示例，LAML并没有APOBEC enrichment，所以不适合做这个分析
plotApobecDiff(tnm = laml.tnm, maf = laml, pVal = 0.2)
## $results
## Hugo_Symbol Enriched nonEnriched pval
## 1: TP53 2 13 0.08175632
## 2: TET2 1 16 0.45739351
## 3: FLT3 2 45 0.65523131
## 4: DNMT3A 1 47 1.00000000
## 5: ADAM11 0 2 1.00000000
## ---
## 132: WAC 0 2 1.00000000
## 133: WT1 0 12 1.00000000
## 134: ZBTB33 0 2 1.00000000
## 135: ZC3H18 0 2 1.00000000
## 136: ZNF687 0 2 1.00000000
## or ci.up ci.low adjPval
## 1: 5.9976455 46.608861 0.49875432 1
## 2: 1.9407002 18.983979 0.03882963 1
## 3: 1.4081851 10.211621 0.12341748 1
## 4: 0.5335362 4.949499 0.01101929 1
## 5: 0.0000000 164.191472 0.00000000 1
## ---
## 132: 0.0000000 164.191472 0.00000000 1
## 133: 0.0000000 12.690862 0.00000000 1
## 134: 0.0000000 164.191472 0.00000000 1
## 135: 0.0000000 164.191472 0.00000000 1
## 136: 0.0000000 164.191472 0.00000000 1
##
## $SampleSummary
## Cohort SampleSize Mean Median
## 1: Enriched 6 7.167 6.5
## 2: nonEnriched 172 9.715 9.0
6 小技巧
6.1 MAF取子集
## 常规筛选
#Extract data for samples 'TCGA.AB.3009' and 'TCGA.AB.2933' (Printing just 5 rows for display convenience)
subsetMaf(maf = laml, tsb = c('TCGA-AB-3009', 'TCGA-AB-2933'), mafObj = FALSE)[1:5]
## Hugo_Symbol Entrez_Gene_Id Center
## 1: ABCB11 8647 genome.wustl.edu
## 2: ACSS3 79611 genome.wustl.edu
## 3: ANK3 288 genome.wustl.edu
## 4: AP1G2 8906 genome.wustl.edu
## 5: ARC 23237 genome.wustl.edu
## NCBI_Build Chromosome Start_Position End_Position
## 1: 37 2 169780250 169780250
## 2: 37 12 81536902 81536902
## 3: 37 10 61926581 61926581
## 4: 37 14 24033309 24033309
## 5: 37 8 143694874 143694874
## Strand Variant_Classification Variant_Type
## 1: + Missense_Mutation SNP
## 2: + Missense_Mutation SNP
## 3: + Splice_Site SNP
## 4: + Missense_Mutation SNP
## 5: + Missense_Mutation SNP
## Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2
## 1: G G A
## 2: C C T
## 3: C C A
## 4: C C T
## 5: C C G
## Tumor_Sample_Barcode Protein_Change i_TumorVAF_WU
## 1: TCGA-AB-3009 p.A1283V 46.97218
## 2: TCGA-AB-3009 p.A266V 47.32510
## 3: TCGA-AB-3009 43.99478
## 4: TCGA-AB-3009 p.R346Q 47.08000
## 5: TCGA-AB-3009 p.W253C 42.96435
## i_transcript_name
## 1: NM_003742.2
## 2: NM_024560.2
## 3: NM_020987.2
## 4: NM_003917.2
## 5: NM_015193.3
##Same as above but return output as an MAF object (Default behaviour)
subsetMaf(maf = laml, tsb = c('TCGA-AB-3009', 'TCGA-AB-2933'))
## An object of class MAF
## ID summary Mean Median
## 1: NCBI_Build 37 NA NA
## 2: Center genome.wustl.edu NA NA
## 3: Samples 2 NA NA
## 4: nGenes 34 NA NA
## 5: Frame_Shift_Ins 5 2.5 2.5
## 6: In_Frame_Ins 1 0.5 0.5
## 7: Missense_Mutation 26 13.0 13.0
## 8: Nonsense_Mutation 2 1.0 1.0
## 9: Splice_Site 1 0.5 0.5
## 10: total 35 17.5 17.5
## 指定query筛选
#Select all Splice_Site mutations from DNMT3A and NPM1
subsetMaf(maf = laml, genes = c('DNMT3A', 'NPM1'), mafObj = FALSE,query = "Variant_Classification == 'Splice_Site'")
## Hugo_Symbol Entrez_Gene_Id Center
## 1: DNMT3A 1788 genome.wustl.edu
## 2: DNMT3A 1788 genome.wustl.edu
## 3: DNMT3A 1788 genome.wustl.edu
## 4: DNMT3A 1788 genome.wustl.edu
## 5: DNMT3A 1788 genome.wustl.edu
## 6: DNMT3A 1788 genome.wustl.edu
## NCBI_Build Chromosome Start_Position End_Position
## 1: 37 2 25459874 25459874
## 2: 37 2 25467208 25467208
## 3: 37 2 25467022 25467022
## 4: 37 2 25459803 25459803
## 5: 37 2 25464576 25464576
## 6: 37 2 25469029 25469029
## Strand Variant_Classification Variant_Type
## 1: + Splice_Site SNP
## 2: + Splice_Site SNP
## 3: + Splice_Site SNP
## 4: + Splice_Site SNP
## 5: + Splice_Site SNP
## 6: + Splice_Site SNP
## Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2
## 1: C C G
## 2: C C T
## 3: A A G
## 4: A A C
## 5: C C A
## 6: C C A
## Tumor_Sample_Barcode Protein_Change i_TumorVAF_WU
## 1: TCGA-AB-2818 p.R803S 36.84
## 2: TCGA-AB-2822 62.96
## 3: TCGA-AB-2891 34.78
## 4: TCGA-AB-2898 46.43
## 5: TCGA-AB-2934 p.G646V 37.50
## 6: TCGA-AB-2968 p.E477* 40.00
## i_transcript_name
## 1: NM_022552.3
## 2: NM_022552.3
## 3: NM_022552.3
## 4: NM_022552.3
## 5: NM_022552.3
## 6: NM_022552.3
#Same as above but include only 'i_transcript_name' column in the output
subsetMaf(maf = laml, genes = c('DNMT3A', 'NPM1'), mafObj = FALSE, query = "Variant_Classification == 'Splice_Site'", fields = 'i_transcript_name')
## Hugo_Symbol Chromosome Start_Position End_Position
## 1: DNMT3A 2 25459874 25459874
## 2: DNMT3A 2 25467208 25467208
## 3: DNMT3A 2 25467022 25467022
## 4: DNMT3A 2 25459803 25459803
## 5: DNMT3A 2 25464576 25464576
## 6: DNMT3A 2 25469029 25469029
## Reference_Allele Tumor_Seq_Allele2
## 1: C G
## 2: C T
## 3: A G
## 4: A C
## 5: C A
## 6: C A
## Variant_Classification Variant_Type
## 1: Splice_Site SNP
## 2: Splice_Site SNP
## 3: Splice_Site SNP
## 4: Splice_Site SNP
## 5: Splice_Site SNP
## 6: Splice_Site SNP
## Tumor_Sample_Barcode i_transcript_name
## 1: TCGA-AB-2818 NM_022552.3
## 2: TCGA-AB-2822 NM_022552.3
## 3: TCGA-AB-2891 NM_022552.3
## 4: TCGA-AB-2898 NM_022552.3
## 5: TCGA-AB-2934 NM_022552.3
## 6: TCGA-AB-2968 NM_022552.3
## 按临床信息筛选
#Select all samples with FAB clasification M4 in clinical data
laml_m4 = subsetMaf(maf = laml, clinQuery = "FAB_classification %in% 'M4'")
6.2 拿到已经做好的TCGA MAF
https://github.com/PoisonAlien/TCGAmutations
数据来自(TCGA firehose)[http://firebrowse.org/] and (TCGA MC3)[https://gdc.cancer.gov/about-data/publications/mc3-2017] projects. Every dataset is stored as an MAF object containing somatic mutations along with clinical information.
# 安装（但是我目前安装不上）
# BiocManager::install(pkgs = "PoisonAlien/TCGAmutations")
# 加载 TCGA LUAD cohort. Note that by default MAF from MC3 study is returned
TCGAmutations::tcga_load(study = "LUAD")
# 切换数据源
TCGAmutations::tcga_load(study = "LUAD", source = "Firehose")
