2020-05-31 学习maftools : 可视化分析SNV

续昨天的学习，不过是重启了一下，临床富集分析又可以做了，但是Signature分析这一步开始出现了新的错误，后面解决。不知道为什么GitHub都打不开连接了，好多资料都找不到。桑心

9.7 临床富集分析

clinicalEnrichment是另一个函数，它提取与样本相关的任何临床特征并进行富集分析。它执行各种groupwise和pairwise成对比较，以确定临床特征中每个类别的丰富突变。下面是一个识别与FAB_classfication相关的突变的示例。

fab.ce = clinicalEnrichment(maf = laml, clinicalFeature = 'FAB_classification')
#> 
#> M0 M1 M2 M3 M4 M5 M6 M7 
#> 19 44 44 21 39 19  3  3
#Results are returned as a list. Significant associations p-value < 0.05
fab.ce$groupwise_comparision[p_value < 0.05]
#>    Hugo_Symbol Group1 Group2 n_mutated_group1 n_mutated_group2
#> 1:        IDH1     M1   Rest         11 of 44         7 of 149
#> 2:        TP53     M7   Rest           3 of 3        12 of 190
#> 3:      DNMT3A     M5   Rest         10 of 19        38 of 174
#> 4:       CEBPA     M2   Rest          7 of 44         6 of 149
#> 5:       RUNX1     M0   Rest          5 of 19        11 of 174
#> 6:        NPM1     M5   Rest          7 of 19        26 of 174
#> 7:       CEBPA     M1   Rest          6 of 44         7 of 149
#>         p_value OR_low   OR_high       fdr
#> 1: 0.0002597371      0 0.3926994 0.0308575
#> 2: 0.0003857187      0 0.1315271 0.0308575
#> 3: 0.0057610493      0 0.6406007 0.3072560
#> 4: 0.0117352110      0 0.6874270 0.3757978
#> 5: 0.0117436825      0 0.6466787 0.3757978
#> 6: 0.0248582372      0 0.8342897 0.6628863
#> 7: 0.0478737468      0 0.9869971 1.0000000

以上结果表明，与队列的其余部分相比，IDH1突变在M1亚型白血病中富集。同样，DNMT3A在M5中，RUNX1在M0中，依此类推。这些都是众所周知的结果，此函数有效地重新捕获了它们。可以使用任何类型的临床特征来执行这样的分析。还有一个小函数-plotEnrichmentResults，可用于绘制这些结果。

plotEnrichmentResults(enrich_res = fab.ce, pVal = 0.05)

image.png

9.8 药物与基因的相互作用。

从药物基因相互作用数据库汇编而来的药物-基因相互作用和基因和药性信息可用drugInteractions函数查询

dgi = drugInteractions(maf = laml, fontSize = 0.75)

image.png

上图显示了潜在的可用药基因类别，以及与之相关的最多5个基因。人们还可以提取关于药物与基因相互作用的信息。例如，下面是已知/报告的药物与DNMT3A相互作用的结果。

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

如果您认为此函数对你有用 10，请引用DGIdb article。 免责声明：此功能中使用的资源仅用于研究目的。它不应用于紧急情况或医疗或专业建议。。 ### 9.9致癌信号通路。 OncogenicPathways函数查看TCGA队列中已知的致癌信号通路的富集情况 11。

OncogenicPathways(maf = laml)
#>        Pathway  N n_affected_genes fraction_affected
#>  1:    RTK-RAS 85               18        0.21176471
#>  2:      Hippo 38                7        0.18421053
#>  3:      NOTCH 71                6        0.08450704
#>  4:        MYC 13                3        0.23076923
#>  5:        WNT 68                3        0.04411765
#>  6:       TP53  6                2        0.33333333
#>  7:       NRF2  3                1        0.33333333
#>  8:       PI3K 29                1        0.03448276
#>  9: Cell_Cycle 15                0        0.00000000
#> 10:   TGF-Beta  7                0        0.00000000

image.png

也有可能将完整的路径可视化。

PlotOncogenicPathways(maf = laml, pathways = "RTK-RAS")

image.png

抑癌基因用红色表示，癌基因用蓝色字体表示。 ### 9.10 肿瘤异质性与数学成绩 #### 9.10.1 肿瘤样本的异质性肿瘤一般是异质性的，即由多个克隆组成。这种异质性可以通过聚类不同的等位基因频率来推断。inserHetereneity函数使用VAF信息对变量进行聚类(使用mclust)，从而推断出克隆性。默认情况下，inserverHetereneity函数查找包含VAF信息的t_vaf列。但是，如果字段名称与t_vaf不同，我们可以使用参数vafCol手动指定。例如，在本例中，研究vaf存储在字段名i_TumorVAF_WU下。

#Heterogeneity in sample TCGA.AB.2972
tcga.ab.2972.het = inferHeterogeneity(maf = laml, tsb = 'TCGA-AB-2972', vafCol = 'i_TumorVAF_WU')
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
#Visualizing results
plotClusters(clusters = tcga.ab.2972.het)

image.png

上图清楚地显示了两个平均变异等位基因频率为_{45%的克隆(主克隆)和另一个次要克隆的变异等位基因频率为}25%。 虽然变异等位基因频率的聚类使我们对异质性有了一个很好的认识，但也可以用数值来衡量异质性的程度。MATH score(上图中的副标题)是肿瘤内异质性的一种简单的定量测量，它计算出VAF分布的宽度。研究发现，较高的MATH scores 与较差的结果有关。MATH score也可以作为生存分析的代理变量-11.

9.10.2 忽略拷贝数更改区域的变异。

我们可以使用拷贝数信息来忽略位于拷贝数改变区域的变异。拷贝数改变会导致变异等位基因频率异常高/低，这往往会影响聚类。去除这些变异可以改进聚类和密度估计，同时保留有生物学意义的结果。拷贝数信息可以作为从分割程序生成的分割文件来提供，例如来自“DNACopy” Bioconductor package的圆形二进制分割的R包6

seg = system.file('extdata', 'TCGA.AB.3009.hg19.seg.txt', package = 'maftools')
tcga.ab.3009.het = inferHeterogeneity(maf = laml, tsb = 'TCGA-AB-3009', segFile = seg, vafCol = 'i_TumorVAF_WU')
#>    Hugo_Symbol Chromosome Start_Position End_Position Tumor_Sample_Barcode
#> 1:        PHF6         23      133551224    133551224         TCGA-AB-3009
#>        t_vaf Segment_Start Segment_End Segment_Mean CN
#> 1: 0.9349112            NA          NA           NA NA
#>    Hugo_Symbol Chromosome Start_Position End_Position Tumor_Sample_Barcode
#> 1:         NF1         17       29562981     29562981         TCGA-AB-3009
#> 2:       SUZ12         17       30293198     30293198         TCGA-AB-3009
#>        t_vaf Segment_Start Segment_End Segment_Mean       CN    cluster
#> 1: 0.8419000      29054355    30363868      -0.9157 1.060173 CN_altered
#> 2: 0.8958333      29054355    30363868      -0.9157 1.060173 CN_altered
#Visualizing results. Highlighting those variants on copynumber altered variants.
plotClusters(clusters = tcga.ab.3009.het, genes = 'CN_altered', showCNvars = TRUE)

image.png

上图显示了两个VAF高的基因NF1和SUZ12，这是由于拷贝数变异(缺失)所致。这两个基因在分析中被忽略了。 ### 9.11 突变Signatures 每一种癌症，在进展过程中都会留下一个Signatures，其Signatures是特定的核苷酸替换模式。 Alexandrov et.al已经显示出这样的突变signatures，源自5的7,000多个癌症样本。这样的signatures可以通过分解核苷酸替换矩阵来提取，基于突变碱基周围的直接碱基将其分类为96个替换类别。还可以将提取的signatures与那些signatures进行比较validated signatures。 signatures分析的第一步是获得突变碱基周围的相邻碱基，并形成突变矩阵。 >注意：早期版本的maftools需要一个fasta文件作为输入。但是从1.8.0开始，BSgenome对象被用于更快的序列提取。

#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")
#> Warning in trinucleotideMatrix(maf = laml, prefix = "chr", add = TRUE, ref_genome = "BSgenome.Hsapiens.UCSC.hg19"): Chromosome names in MAF must match chromosome names in reference genome.
#> Ignorinig 101 single nucleotide variants from missing chromosomes chr23
#> -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

上述功能执行两个步骤： * 估计APOBEC富集分数 * 为signature分析准备突变矩阵。

9.11.1 APOBEC Enrichment estimation.

APOBEC诱导的突变在实体瘤中更为常见，并且主要与TCW基序中发生的C>T转换事件有关。上述命令中的APOBEC富集分数是使用Roberts et al 13描述的方法估计的。简而言之，将发生在TCW基序内的C>T突变富集在给定样本中的所有C>T突变上，并将其与背景胞嘧啶和发生在突变碱基20bp内的TCW进行比较。

image.png

如Roberts等人的原始研究中所描述的，还执行单侧fishers精确测试以统计地评估富集分数。

9.11.2 APOBEC富集与非富集的区别。

我们还可以分析APOBEC富集型和非APOBEC富集型的突变模式的差异，plotApobecDiff是一个函数，它取trinucleotideMatrix估计的APOBEC富集值，将样本分为APOBEC富集型和非APOBEC富集型。一旦分层，它就会比较这两组人，以确定差异改变的基因。 请注意，没有APOBEC富集的LAML不是此类分析的理想队列，因此下面的图表仅用于演示目的。

plotApobecDiff(tnm = laml.tnm, maf = laml, pVal = 0.2)

image.png

#> $results
#>      Hugo_Symbol Enriched nonEnriched       pval        or      ci.up
#>   1:        TP53        2          13 0.08175632 5.9976455  46.608861
#>   2:        TET2        1          16 0.45739351 1.9407002  18.983979
#>   3:        FLT3        2          45 0.65523131 1.4081851  10.211621
#>   4:      DNMT3A        1          47 1.00000000 0.5335362   4.949499
#>   5:      ADAM11        0           2 1.00000000 0.0000000 164.191472
#>  ---                                                                 
#> 132:         WAC        0           2 1.00000000 0.0000000 164.191472
#> 133:         WT1        0          12 1.00000000 0.0000000  12.690862
#> 134:      ZBTB33        0           2 1.00000000 0.0000000 164.191472
#> 135:      ZC3H18        0           2 1.00000000 0.0000000 164.191472
#> 136:      ZNF687        0           2 1.00000000 0.0000000 164.191472
#>          ci.low adjPval
#>   1: 0.49875432       1
#>   2: 0.03882963       1
#>   3: 0.12341748       1
#>   4: 0.01101929       1
#>   5: 0.00000000       1
#>  ---                   
#> 132: 0.00000000       1
#> 133: 0.00000000       1
#> 134: 0.00000000       1
#> 135: 0.00000000       1
#> 136: 0.00000000       1
#> 
#> $SampleSummary
#>         Cohort SampleSize  Mean Median
#> 1:    Enriched          6 7.167    6.5
#> 2: nonEnriched        172 9.715    9.0

9.11.3 Signature分析。

Signature分析包括以下步骤。 1.estimateSignatures -在一系列值上运行NMF，并衡量拟合的好坏-以时间为单位Cophenetic correlation。 2.plotCophenetic-绘制elblow图，帮助您确定最佳Signature数量。最佳可能的Signature是共生相关性显著下降的值。 3.ExtractSignaturesSignature-使用非负矩阵分解将矩阵分解为三个N个Signature，根据上述两个步骤选择N个Signature。如果您已经对N‘有了很好的估计，您可以跳过以上两步。 4.将上述步骤中提取的SignaturecompareSignatures与signatures[11](http://127.0.0.1:25995/library/maftools/doc/maftools.html#references)数据库中已知的[COSMIC](https://cancer.sanger.ac.uk/cosmic/signatures/SBS/)签名进行比对，并计算余弦相似度，确定最佳匹配。 5.plotSignatures`-plots signatures

image.png

注意：在以前的版本中，上述步骤都是由ExtractSignatures自动完成的。在2.2.0版本之后，Main函数被拆分成不超过5个stpe，以方便用户灵活使用。

library('NMF')
#laml.sign = estimateSignatures(mat = laml.tnm, nTry = 6)

绘制elbow曲线图，根据上述结果可视化并确定最佳signatures数量。

#plotCophenetic(res = laml.sign)

最佳可能值是y轴上的相关值显著下降的值。在这种情况下，它看起来是在n = 3。LAML的突变率较低，不是特征分析的理想例子，但对于突变负担较高的实体肿瘤，只要有足够数量的样本，就可以期待更多的特征。 一旦估计了n，我们就可以运行main函数了。

#laml.sig = extractSignatures(mat = laml.tnm, n = 3)

将检测到的signatures与COSMIC数据库中的已知signatures进行比较。

#Compate against original 30 signatures 
#laml.og30.cosm = compareSignatures(nmfRes = laml.sig, sig_db = "legacy")
#Compate against updated version3 60 signatures 
#laml.v3.cosm = compareSignatures(nmfRes = laml.sig, sig_db = "SBS")

compareSignatures返回COSMIC signatures余弦相似度的完整表，可以进一步分析。下图显示了检测到的signatures与验证过的signatures的相似性比较。

library('pheatmap')
#pheatmap::pheatmap(mat = laml.og30.cosm$cosine_similarities, cluster_rows = FALSE, main = "cosine similarity against validated signatures")

Finally plot signatures

#maftools::plotSignatures(nmfRes = laml.sig, title_size = 0.8)

注： 1.如果您在运行extractSignatures时收到none of the packages are loaded的错误，请手动加载NMF库并重新运行。 2.如果extractSignatures或estimateSignatures在两者之间停止，则可能是因为矩阵中的突变计数较低。在这种情况下，重新运行pConstant参数设置为小正值(例如0.1)的函数。

9.11.4 Signature enrichment analysis

Signature可以进一步赋值给样本，并使用signatureEnrichment函数进行富集分析，该函数识别每个识别出的Signature中富集的突变。

#laml.se = signatureEnrichment(maf = laml, sig_res = laml.sig)

上述结果可进行和临床结果相似的可视化操作。

#plotEnrichmentResults(enrich_res = laml.se, pVal = 0.05)