DNA 4. SCI 文章中基因组的突变信号（maftools）

点击关注，桓峰基因

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

87篇原创内容

公众号

突变信号（Mutational Signatures）首次2013年在《nature》进行报道。并做了相关的定义：细胞在成长过程中，基因组不断受到内源性和外源性DNA损伤的威胁，正是由于这些威胁，使得细胞基因组不断发生变化，并最终发生一些突变的积累。每一个突变过程都会留下一个不同的基因组标记，也就称为突变信号。

有需要上述分析内容的老师可以加我微信：

前言

研究人员对TCGA数据的深度挖掘，从而提出的一个统计学概念。文章研究了30种癌症，发现21种不同的mutation signature。如果理解了，就会发现这个其实蛮简单的，他们并不重新测序，只是拿已经有了的TCGA数据进行分析，而且居然是发表在nature上面！

研究了4,938,362 mutations from 7,042 cancers样本，突变频谱的概念只是针对于somatic 的mutation。一般是对癌症病人的肿瘤组织和癌旁组织配对测序，过滤得到的somatic mutation，一般一个样本也就几百个somatic 的mutation。

基本概念

Signature 概念

突变信号（Mutational Signatures）是细胞在成长过程中，基因组不断受到内源性和外源性DNA损伤的威胁，正是由于这些威胁，使得细胞基因组不断发生变化，并最终发生一些突变的积累。每一个突变过程都会留下一个不同的基因组标记，也就称为突变信号。今天突变信号分析已成为基因组研究的标准组成部分，因为它们可以揭示每个肿瘤中突变的环境和内源性来源。事实上，这一新生领域正变得越来越突出，并朝着临床意义上的应用方向发展。虽然这些都是积极的趋势，但这一领域的大幅扩大是否存在局限性？随着越来越多的不同突变类别的信号被报道，它们与各种因素之间的相关性正在被解释。然而，许多信号的起源仍然是神秘的。

APOBEC 富集

我们先来介绍一下APOBEC（apolipoprotein B mRNA editing enzyme catalytic polypeptide like）富集的意义。APOBEC家族(载脂蛋白BmRNA 编辑酶催化多肽家族)是一种具有抗病毒活性的蛋白分子。先前的研究认为APOBEC家族的亚型APOBEC3G具有抗艾滋病病毒的能力。

家族成员:

老大:APOBEC1

老二: APOBEC2

老三:APOBEC3A-G

老四:AID（活化诱导脱氨酶)

APOBEC1 、AID以串联重复形式排列于第12号染色体; APOBEC2位于第6号染色体; APOBEC3A-G以串联重复形式排列于第22号染色体。APOBEC1:最早发现的 APOBEC 家族成员。参与载脂蛋白 B mRNA 的编辑，是 RNA编辑典型范例之一。RNA 编辑:转录后mRNA发生特异性位点的核苷酸改变,结果产生无义突变或错义突变而导致mRNA编码容量的变化。

主要表达于肠道和某些哺乳动物的肝脏内，能导致apoB(载脂蛋白B)mRNA特异性位点单个核苷酸C (C6666)转变为U，使得编码谷氨酰胺的密码子CAA转变为终止密码子UAA ，导致翻译提前终止，生成 apoB100蛋白的截短形式apoB48 。在大肠杆菌的胞苷脱氨酶中发现APOBEC1有保守的活性位点功能域。此外,它有一个独特的功能域含有两个苯丙氨酸残基和插入的4个氨基酸残基跨越活性位点功能域。该功能域在APOBEC家族包括激活诱导的胞苷脱氨酶（AID ) ，APOBEC2，APOBEC3A和APOBEC3G在内的成员中均存在。APOBEC3 : 主要表达于淋巴细胞和骨髓细胞谱系。APOBEC3G 基因位于人的第22号染色体长臂22q13.1-q13.2，含有8个外显子和 7个内含子，编码384 个氨基酸，其中第128 - 194位和第320-380位氨基酸残基为两个重要的活性功能区:锌离子结合功能区和水解酶功能区除APOBEC3G外，APOBEC3 家族其它成员APOBEC3B 、 APOBEC3F等均有很强的抗逆转录病毒活性 AID(活化诱导脱氨酶) 它在启动抗体类型转换重组，体细胞高突变和基因转化中至关重要。通过将免疫球蛋白位点上的胞喀啶脱氨基化来启动这些过程，具有潜在的基因突变活性在大肠杆菌中表达AID，发现细菌对抗抗生素的能力上升了，就是说基因突变率增加了，突变的倾向也是C-T和G-A，这说明AID确实具有对脱氧胞苷的脱氨能力。AID/APOBEC蛋白家族的许多成员都能够将单链多核苷酸上的胞喀啶脱氨基化为尿嚰啶，从而实现各种各样的生理功能。该家族的酶的过度表达会导致癌症，提示编码APOBEC家族的蛋白的基因是原癌基因，与细胞增殖密切相关。它们能够在DNA或RNA水平上改变病毒的遗传信息,这一称为编辑(editing)的修饰和加工过程,可以在多种病毒的基因组或其逆转录产物中引入高频突变，进而诱导其降解、干扰其复制或者严重影响病毒蛋白的生物学功能。研究发现人体细胞内的特定酵素可增强蛋白质APOBEC抑制HIV病毒繁殖的能力。APOBEC通过使HIV的遗传基因的变异来阻碍其繁殖。通常情况下HIVC会在分泌病毒感染因子蛋白(Vif )，分解APOBEC，所以APOBEC不能发挥足够的防御力。将酵素注射进感染了HIV病毒的细胞中，可抑制Vif分解APOBE，促进 APOBEC发挥抑制HIV增殖的作用。APOBEC3G会在HIV进行反转录成cDNA的阶段时，造成其cDNA上碱基的突变（胞嚰啶(cytidine)转变成尿嚰啶(uraci1 ))，使得HIV的感染力下降。APOBEC3G及其家族成员在除HIV 、 HBV外的其它病毒的影响: APOBEC1家族成员之一的hA1 ( ratAPOBEC1，)对DNA 、RNA均有编辑作用，提示该因子对包括RNA病毒在内的其它病毒也有潜在的抑制作用;除胞喀啶脱氨酶作用外， APOBEC3G在免疫系统中是否具有其它作用机制;APOBEC3G 家族在人类肿瘤性疾病中的作用; APOBEC3G与机体其它细胞因子相互作用等问题是目前研究的热点。对上述机制的研究对于阐明体内固有的对病毒的限制因素至关重要并可在将来进一步对病毒进行药物干预研究产生重要影响。

分析步骤

estimatsignatures——它在一系列值上运行NMF，并度量拟合度——就Cophenetic相关性而言;
plotCophenetic 绘制一个elblow plot，并帮助您决定Signatures的最佳数量。最好的可能标志是Cophenetic相关性显著下降的值;
extractSignatures—使用非负矩阵分解将矩阵分解为n个签名。根据以上两步选择N。如果你已经对n有一个很好的估计，你可以跳过以上两步;
compareSignatures -从以上步骤提取的Signatures可以与COSMIC数据库中的已知签名11进行比较，并计算余弦相似度以确定最佳匹配;
plotSignatures—绘制Signatures。

实例讲解

1. 软件安装

在安装这个软件maftools时，需要先安装BioManager,然后在安装maftools，BSgenome.Hsapiens.UCSC.hg19，NMF，barplot3d，如下：

if (!require("BiocManager")) install.packages("BiocManager")if (!require("maftools")) BiocManager::install("maftools")if (!require("BSgenome.Hsapiens.UCSC.hg19")) BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")if (!require("NMF")) BiocManager::install("NMF")if (!require("barplot3d")) BiocManager::install("barplot3d")if (!require("pheatmap")) install.packages("pheatmap")library("maftools")
library("BSgenome.Hsapiens.UCSC.hg19", quietly = TRUE)
library("NMF")
library("pheatmap")
library("barplot3d")

2. 数据读取

maftools工具需要读入两个文件，如下：

1.MAF文件-可以是gz压缩。必需的；

2.与MAF中每个样本/肿瘤样本条码相关的可选但推荐的临床数据；

3.一个可选的拷贝数数据：可以是GISTIC输出或自定义表。

# path to TCGA LAML MAF file
laml.maf = system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
# clinical information containing survival information and histology. This is
# optional
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 5.440s elapsed (1.550s 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

3. 实例操作

第一步：Signature分析的第一步是获取突变基周围的相邻基，形成突变矩阵。注意:maftools的早期版本需要一个fasta文件作为输入。但是从1.8.0开始，BSgenome对象被用于更快的序列提取。包括两部：

Estimates APOBEC enrichment scores;
Prepares a mutational matrix for signature analysis.

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

APOBEC 富集分析

APOBEC诱导的突变在实体肿瘤中更为常见，主要与TCW motif中发生的C>T转换事件有关。使用Roberts等人所描述的方法计算上述命令中的APOBEC富集分数。简单地说，在一个给定的样本中，将TCW motif中发生的C>T突变与所有C>T突变的富集情况的比例与背景胞密啶和发生在突变碱基20bp内的TCWs进行比较。

我们还可以分析APOBEC富集和非APOBEC富集的样品在突变模式上的差异。plotApobecDiff 函数采用trinucleotideMatrix计算APOBEC富集分数，将样本分为APOBEC富集和非APOBEC富集。分组后，比较这两组，以确定改变的基因的差异。

plotApobecDiff(tnm = laml.tnm, maf = laml, pVal = 0.5)
## -Processing clinical data
## -Processing clinical data

## $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

APOBEC 富集分析

NMF 算法

非负矩阵分解（NMF）是一种无监督学习算法，其目的在于提取有用的特征。它的工作原理类似于PCA，也可以用于降维。与PCA相同，我们试图将每个数据点写成一些分量的加权求和。但在PCA中，我们想要的是正负分量，并且能够解释尽可能多的数据方差；而在NMF中，我们希望分量和系数均为负，也就是说，我们希望分量和系数都大于或等于0。因此，NMF只能应用于每个特征都是非负的数据，因为非负分量的非负求和不可能变为负值。

将数据分解成非负加权求和的这个过程，对由多个独立源相加（或叠加）创建而成的数据特别有用，比如多人说话的音轨或包含很多乐器的音乐。在这种情况下，NMF可以识别出组合成数据的原始分量。总的来说，与PCA相比，NMF得到的分量更容易解释，因为负的分量和系数可能会导致难以解释的抵消效应。

从图上可以看到，在3的时候是最优的选择，如下：

library("NMF")
laml.sign = estimateSignatures(mat = laml.tnm, nTry = 6, pConstant = 0.1)
## -Running NMF for 6 ranks
## Compute NMF rank= 2  ... + measures ... OK
## Compute NMF rank= 3  ... + measures ... OK
## Compute NMF rank= 4  ... + measures ... OK
## Compute NMF rank= 5  ... + measures ... OK
## Compute NMF rank= 6  ... + measures ... OK
## -Finished in 39.1s elapsed (10.8s cpu)
plotCophenetic(res = laml.sign)

提取n=3的最优数据，然后与COSMIC数据库进行比较，如下：

laml.sig = extractSignatures(mat = laml.tnm, n = 3, pConstant = 0.1)
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")

4. 结果展示

1. 绘制热图

热图显示了检测 Signatures 与验证 Signatures 的相似性比较,如下：

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

2.绘制 Signatures

plotSignatures(nmfRes = laml.sig, title_size = 1.2, sig_db = "SBS")

3. 绘制3D效果图

library("barplot3d")
# Visualize first signature
sig1 = laml.sig$signatures[, 1]
barplot3d::legoplot3d(contextdata = sig1, labels = FALSE, scalexy = 0.01, sixcolors = "sanger",alpha = 0.5)

结果解读

注意事项：

如果您在运行extractSignatures时收到一个错误，抱怨没有加载任何包，请手动加载NMF库并重新运行。
如果extractSignatures或estimatessignatures在中间停止，可能是因为矩阵中的突变计数很低。在这种情况下，重新运行将pConstant参数设置为小正值(例如0.1)的函数。

References:

Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. 2018. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Resarch.
Revathidevi S, Murugan AK, Nakaoka H, Inoue I, Munirajan AK. APOBEC: A molecular driver in cervical cancer pathogenesis. Cancer Lett. 2021;496:104-116. doi:10.1016/j.canlet.2020.10.004
Alexandrov LB, Nik-Zainal S, Wedge DC, et al. Signatures of mutational processes in human cancer [published correction appears in Nature. 2013 Oct 10;502(7470):258. Imielinsk, Marcin [corrected to Imielinski, Marcin]]. Nature. 2013;500(7463):415-421. doi:10.1038/nature12477