样本编号(客户):57997-11

样本编号(公司):NovaS02510-BG0D25X096-25X24001-T2691V1hg38-6-CX-20251104cap59-Q1P1C2-25R39916

报告日期:20251113

版本号:TA2024031

T2691V1hg38-单样本分析报告

艾吉泰康生物科技有限公司

一、技术背景

TargetSeq®液相探针杂交捕获技术是由iGeneTech自主研发,基于热力学稳定性算法对目标区域基因组进行探针设计,然后合成出有效的特异性探针,进而与基因组DNA进行液相杂交,对目标区域序列进行捕获并富集后,使用主流的测序平台进行高通量测序。通过对目标区域测序,可以对候选基因或候选位点进行检测,也可以进一步的找到易感位点,结合数据库信息,更好的了解变异之间的关联和疾病的致病机理。

1

二、实验流程

2.1 DNA样本检测

TargetSeq®对 DNA 样品的检测主要包括 2 种方法:

1)琼脂糖凝胶电泳分析 DNA 降解程度以及是否有 RNA 或蛋白等污染。

2)Qubit对DNA浓度进行精准定量。

2.2 建库捕获

将自定义的捕获区域,采用艾吉泰康自主研发的液相捕获技术进行高效、特异富集,富集后在Novaseq 6000、Nextseq 500、DNBSEQ-T7、MGISEQ-2000等主流二代测序平台上进行高通量、高深度测序。其操作流程如下:

1)将基因组DNA经酶切或超声随机打断成150-250 bp的片段;

2)片段化后的DNA进行末端修复、3'端加“A”;

3)接头连接;

4)通过PCR扩增的方法进行样本标记和富集DNA;

5)带有特异index的文库与生物素标记的RNA探针进行液相杂交,利用生物素-亲和素间的结合具有迅速、专一、稳定的特性,再使用链酶亲和素标记的磁珠获取目标基因外显子,然后用PCR扩增进行目标基因的富集。

图1 目标区域建库捕获流程图

2.3 文库质检定量

1)使用 Qubit 3.0 进行文库定量,文库浓度>25 ng/μl 参考为合格文库;

2)使用 Qsep 100 检测,文库主峰要在 220-450 bp 左右,主峰前后无杂峰。

2.4 上机测序

1)取1 μl文库使用 Qubit dsDNA HS Assay Kit 进行定量,记录文库浓度,文库浓度约在1-20 ng/μl;

2)取1 μl样品使用 Qsep 100 进行文库片段长度测定,文库长度约在220-450 bp之间;

3)使用高通量测序平台进行测序。

三、生信分析

获得原始下机数据(一般为Fastq格式)后,基于已知的参考基因组进行生信分析,本产品使用的参考基因组可能为GRCh37/hg19(线粒体版本NC_012920)或者GRCh38/hg38(线粒体版本NC_012920),分析内容主要包括以下几个部分:

1)质量控制及评估:主要包含测序质量评估以及文库质量评估两部分,评估的过程中通常会伴随质量控制,去掉测序接头及部分低质量的序列。评估结果通常体现在一些重要性能指标上,如Q20、Q30、比对率、重复率、覆盖率、捕获率、均一性等指标,可有效评估该数据是否合格,能否进行后续分析。

2)变异检测:将通过质量评估的序列比对到参考基因组上,比较真实数据与参考基因组的不同之处(即变异),该步骤的灵敏度与准确性对本次检测的最终结果影响巨大,同时选择不同版本的参考基因组也会对结果产生些许差异。

3)变异注释:对检出的变异结果进行功能性、人群频率、疾病相关性等数据库的注释,基于注释的结果,挑选出致病或者可能致病的变异。

图2 生信分析流程图

四、质量控制及评估

4.1 原始数据说明

高通量测序(如NovaSeq 6000、NextSeq 500、DNBSEQ-T7、MGISEQ-2000等测序平台)测序得到的原始图象数据文件经碱基识别 (Base Calling)分析转化为原始测序序列 (Sequenced reads),我们称之为Raw data或Raw reads,结果以FASTQ(简称为fq)文件格式存储,其中包含测序序列 (reads)的序列信息以及其对应的测序质量信息。如下图示例:

 @EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
 GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
 +
 !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

1:“@”开头,随后为Illumina 测序标识符 (Sequence Identifiers)和描述文字(选择性部分)

2: 碱基序列

3:“+”开头,随后为Illumina 测序标识符(选择性部分)

4: 对应序列的测序质量

表1 Illumina 测序标识符详细信息表
指标 说明
EAS139 the unique instrument name
136 the run id
FC706VJ the flowcell id
2 flowcell lane
2104 tile number within the flowcell lane
15343 'x'-coordinate of the cluster within the tile
197393 'y'-coordinate of the cluster within the tile
1 the member of a pair, 1or2(paired-end or mate-pair reads only)
Y Y if the read is filtered, N otherwise
18 0 when none of the control bits are on, otherwise it is an even number
ATCACG index sequence

目前大多数平台使用的是phred33测序质量值系统(部分使用phred64质量值系统),在该质量值系统中,每个字符对应的ASCII值减去33,即为对应碱基的测序质量值。如果测序错误率用e表示,碱基质量值用Qphred表示,则有下列关系:

公式:Qphred= -10 log10(e);

此公式表示质量值越大测序错误率越低,准确性越高。质量值为10时准确率为90%,20为99%,30为99.9%,依此类推。

4.2 测序数据过滤

测序得到的原始数据中,会存在少量包含接头(亦称Adapter)、低质量碱基或N碱基的序列,为了保证信息分析质量,需要对Raw reads进行初步过滤,得到Clean reads,后续分析都基于Clean reads进行。 数据过滤的内容及步骤主要如下:

(1)以8 bp滑窗的方式,剪切掉平均碱基质量值小于20的序列;

(2)去掉序列尾部的接头序列;

(3)如果序列首碱基或尾碱基质量值小于20,则会直接剪切掉该碱基;

(4)通常情况下,去除完后若剩余序列长度小于40bp(双端),则丢弃该组序列。

4.3 测序质量统计

原始数据进行初步过滤后,我们会对剩余的数据进行详细的统计(QC.stat.xls文件),主要统计数据量、Q20、Q30、QCrate等指标,一般过滤后的数据,Q20平均比例在90%以上,Q30平均比例在80%以上,属于比较正常的结果:

表2 测序质量统计表
指标 数值 说明
Raw reads(M) 17.54 原始数据总序列数
Raw bases(Mb) 2630.56 原始数据总碱基数
Raw Q20-R1(%) 97.89 原始R1数据中质量值>=20的碱基占比
Raw Q30-R1(%) 94.27 原始R1数据中质量值>=30的碱基占比
Raw Q20-R2(%) 97.18 原始R2数据中质量值>=20的碱基占比
Raw Q30-R2(%) 92.92 原始R2数据中质量值>=30的碱基占比
Raw Q20-all(%) 97.53 原始双端数据中质量值>=20的碱基占比
Raw Q30-all(%) 93.6 原始双端数据中质量值>=30的碱基占比
Clean reads(M) 17.36 原始数据过滤后剩余总序列数
Clean bases(Mb) 2510.96 原始数据过滤后剩余总碱基数
Clean Q20-R1(%) 98.11 过滤后R1数据中质量值>=20的碱基占比
Clean Q30-R1(%) 94.61 过滤后R1数据中质量值>=30的碱基占比
Clean Q20-R2(%) 97.73 过滤后R2数据中质量值>=20的碱基占比
Clean Q30-R2(%) 93.67 过滤后R2数据中质量值>=30的碱基占比
Clean Q20-all(%) 97.92 过滤后双端数据中质量值>=20的碱基占比
Clean Q30-all(%) 94.14 过滤后双端数据中质量值>=30的碱基占比
QCrate(%) 95.45 过滤比率,Clean bases(Mb)/ Raw bases(Mb)

4.4 碱基质量值分布

对于二代测序,测序碱基的质量值通常会随着测序的进行而降低,特别是read2的末端,试剂的性能下降会比较明显。此外每条read的前几个碱基的质量值也可能存在较大波动。我们建议测序数据的质量主要分布至少在Q30(≥80%)以上,这样才能保证后续分析的正常进行。

R1碱基质量值分布
R2碱基质量值分布
图3 碱基质量值分布图(左图为R1,右图为R2)

4.5 四种碱基比例分布

测序过程中,ATCG四种碱基的分布比例会逐渐趋于平衡,只有前几个碱基可能存在较大的波动。另外利用多重PCR技术产生的文库,由于技术特性的原因,四种碱基的分布会呈现出不均衡的状态。

R1碱基比例分布
R2碱基比例分布
图4 碱基比例分布图(左图为R1,右图为R2)

4.6 文库质量统计

经过第一步质量控制Clean data以后,我们采用BWA软件将Clean data与参考基因组做比对,同时进行标记重复、排序等工作,得到bam结果文件。最后我们对bam结果文件进行文库质量统计,统计内容包括文库平均长度、比对率、覆盖率、捕获率、测序深度、均一性等指标,用于评价样本文库质量的好坏,统计的结果存放在statistic.xls中。此外我们还统计了一些不太常用的指标,例如MQ0比例、合适比对率、错误率等,结果存放在statistic_detail.xls中。

表3 文库质量统计表
指标 数值 说明
Raw bases(Mb)2630.56原始数据总碱基数
Clean bases(Mb)2510.96原始数据过滤后剩余总碱基数
QCrate(%)95.45过滤比率,Clean bases(Mb)/ Raw bases(Mb)
Target average read length144目标区域平均序列长度
Target average base quality36.1目标区域平均碱基质量
Target average insert size210.0目标区域平均插入片段长度
Target duplicated bases(Mb)605.44目标区域重复碱基数
Target duplication rate(%)32.08重复率,目标区域重复碱基数/目标区域总碱基数
Total mapped reads(M)17.35总共比对到参考基因组上的序列数
Total reads mapping rate(%)99.94比对率,比对到参考基因组上序列数/质控后总序列数
Target mapped reads(M)13.11比对到目标区域上的序列数
Target reads capture rate(%)75.56捕获率,目标区域序列数/比对到参考基因组上序列数(分子分母均包含重复序列)
Flank mapped reads(M)14.54比对到目标区域及其周围150bp上的序列数
Flank reads capture rate(%)83.85捕获率,目标区域及其周围150bp总序列数/比对到参考基因组上序列数(分子分母均包含重复序列)
Target size513516目标区域大小
Target covered size513516目标区域实际覆盖到的大小
Coverage rate(%)100.0覆盖率,目标区域实际覆盖到的大小/目标区域大小
Target effective bases(Mb)908.86目标区域上的有效碱基数,去除了重复序列碱基,且横跨目标区域内外的序列只截取目标区域内部分
Target effective rate(%)53.06有效碱基率,目标区域上的有效碱基数/比对去重后总的有效碱基数
Target mean depth1769.87目标区域平均深度,目标区域有效碱基数/目标区域大小
T 4X coverage rate(%)100.0深度>=4X的目标区域比率,也叫4X覆盖度
T 10X coverage rate(%)100.0深度>=10X的目标区域比率,也叫10X覆盖度
T 20X coverage rate(%)99.99深度>=20X的目标区域比率,也叫20X覆盖度
T 30X coverage rate(%)99.96深度>=30X的目标区域比率,也叫30X覆盖度
T 10%X coverage rate(%)99.83深度>=平均深度10%X的目标区域比率(假设平均深度200层,就是大于等于20层的比率)
T 20%X coverage rate(%)99.54深度>=平均深度20%X的目标区域比率
T 30%X coverage rate(%)98.77深度>=平均深度30%X的目标区域比率
T 50%X coverage rate(%)92.52深度>=平均深度50%X的目标区域比率
T Mean depth flank 10%X coverage rate(%)21.59平均深度两侧10%深度的目标区域比率(假设平均深度100层,就是90-110层的比率)
Panel nameT2691V1hg38目标区域编号

4.7 文库片段大小分布

文库片段大小对测序的质量存在较大影响,按照PE150的测序模式,文库片段在300bp左右最为合适,但过长的文库可能会导致实验困难,因此我司文库片段的波峰一般在150-250bp之间,呈钟型分布。

文库片段大小分布图
图5 插入片段大小分布图

4.8 测序深度分布

测序深度分布图直观的展示了深度比例分布情况,我们可以很容易地从图上判断出覆盖度、均一性等指标的好坏。

深度积累分布图
深度分布图
图6 测序深度分布图(左图深度积累分布图,右图深度分布图)

注:左图为深度积累分布图,目标区域不同测序深度(均一化处理,平均深度为1)的累积占比,比如深度为0.2时坐标为90,代表90%的目标区域测序深度大于等于平均深度的20%;右图为深度分布图,目标区域不同测序深度(均一化处理,平均深度为1)的占比,比如深度为1时纵坐标为0.5,代表0.5%的目标区域测序深度等于该样本平均深度,该图一般在平均深度周围呈泊松分布。

4.9 单区间质量统计

除了整体的统计结果外,我们还对单个目标区域进行了相关指标的统计,以方便详细查看每个目标区域的覆盖情况,以下是单个目标区域统计格式示例(该表由于结果较多,仅举例说明格式,真实结果请查看regioncount.xls文件):

表4 单区间质量统计表
指标 示例 说明
Chrchr1染色体编号
Start12080区间起始
End12251区间终止
GeneDDX11L1区间覆盖到的基因
Exon_SN1区间覆盖到的外显子编号
TranscriptionNR_046018区间覆盖到的转录本编号
Length(bp)171区间的长度
Target covered size171实际覆盖到的区间大小
Coverage rate(%)100覆盖率
Target bases27314该区间包含的碱基总数
Target mean depth159.73该区间的平均测序深度
Normalized depth1.0524该区间均一化后的平均测序深度
(N)X coverage rate(%)100该区间深度>=(N)X的比率

4.10 单区间深度分布

下图展示了每个区间的深度分布情况,可以直观看出好坏区域的占比。

深度区间分布图
图7 单区间深度分布图

注:图中横坐标为每个区间编号,纵坐标为该区间的深度(经过排序及均一化处理,平均深度值为0)。

4.11 样本性别质控

本司精心选取的5个判定性别位点,可有效判断出样本的测序性别与送样性别是否一致,客户可依据判定文件自行判断。

五、标准分析结果

5.1 SNP检测和注释

SNP(单核苷酸多态性)主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。它是人类可遗传的变异中最常见的一种。占所有已知多态性的90%以上。SNP在人类基因组中广泛存在,平均每500~1000个碱基对中就有1个,估计其总数可达300万个甚至更多。找出SNP以后我们会对结果进行详细的注释,注释内容和说明如下:

特别声明:哺乳动物的X染色体和Y染色体具有共同的进化起源,并保留着高度相似的序列区域。所以就算是女性,也可能存在一些低质量的reads比对到Y染色体,并检测到一些突变(见参考文献10)。

表5 突变(SNP/InDel)结果表
指标 示例 说明
Prioritylow突变的优先级,分为 high | likeyhigh | medium | low 四级,仅供客户参考,其中:high: 位于基因的exonic区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值小于0.01,造成氨基酸改变,且至少一款预测软件预测为有害;likeyhigh: 位于基因的exonic或splicing区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值至少有一个小于0.01,且至少一款预测软件预测为有害;medium: 位于基因的exonic或splicing区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值至少有一个小于0.01;low: 不属于上述3种的结果
Qualityhighquality质量值评分,从高到低分为highquality、mediumquality、lowquality三级,其中highquality代表深度大于等于50层,基因型质量值大于30,频率偏差不超过0.1(杂合0.4-0.6,纯合0.9-1.0)的结果;mediumquality代表深度大约20层,基因型质量值大于20,频率偏差不超过0.2(杂合0.3-0.7,纯合0.8-1.0)的结果,剩余的结果均判定为lowquality。此外位于repeat区域(genomicSuperDups注释结果不为“.”)的结果均会判定为lowquality
Locationontarget突变位置判定,分为ontarget、flank150、offtarget三种,ontarget代表结果位于目标区域内;flank150代表结果位于目标区域周围150bp内;offtarget为目标区域外的结果,该类结果会被直接丢弃,不在结果文件中显示
Mut_typeSNP突变类型,分为SNP和InDel两类
Chrchr1变异所处的染色体
Start11856378变异的起始位置
End11856378变异的终止位置
RefG参考基因组的碱基类型
AltA突变的碱基类型
Vcf_mutG/A突变在vcf中的情况,有些多基因型的位置可以通过该信息看出来
GT0/1基因型,可能会出现0/1(杂合子),1/1(纯合子),0/0(野生型)等情况
AD338/271支持参考基因组的碱基数/支持突变的碱基数
DP609该位置的总深度
GQ255基因型质量值,分值越高越好,一般建议大于30以上
FRE0.445突变频率,支持突变的碱基数/(支持参考基因组的碱基数+支持突变的碱基数)
Func.refGeneexonic参考refGene数据库,对该变异位点所在的区域进行注释 (exonic, splicing, UTR5, UTR3, intronic, ncRNA_exonic, ncRNA_intronic, ncRNA_UTR3, ncRNA_UTR5, ncRNA _splicing, upstream, downstream, intergenic)。说明:1、exonic应该包括coding exonic portion、UTR3和UTR5,但本注释结果中exonic只代表coding exonic portion。2、当一个变异位点位于多个基因或转录本,且功能不同,这些功能按照优先级排序,该列输出优先级最高的功能类型:Exonic = splicing > ncRNA> > UTR5/UTR3 > intron > upstream/downstream > intergenic。当一个变异既位于一个基因的UTR3,又位于另一个基因的UTR5时,该列输出"UTR5,UTR3"。当一个变异既位于一个基因的downstream,又位于另一个基因的upstream时,该列输出" upstream,downstream "
Gene.refGeneMTHFR参考refGene数据库,列出突变所处基因,如果不在基因上,则列出最近的基因名(有可能为两个)
GeneDetail.refGene.参考refGene数据库,描述UTR、splicing、ncRNA_splicing 、 intergenic区域的变异情况,例如:NM_024758:exon3:c.273-5T>C。其中NM_024758表示该变异所在的转录本,exon3表示该变异位于第三个外显子上,c.273-5T>C表示cDNA的273位置上游5bp处发生由T到C的突变;另例:dist=1000;dist=2000。表示该变异位于基因间,距离上下游基因的距离分别为1000 bp和2000 bp
ExonicFunc.refGenenonsynonymous SNV参考refGene数据库,对突变功能类型进行注释,一般位于exonic上的突变才有该信息。其中SNP的突变功能类型包括:nonsynonymous SNV非同义突变、synonymous SNV同义突变、stopgain终止子获得、stoploss终止子丢失;InDel的突变功能类型包括:frameshift insertion移码插入、frameshift deletion移码缺失、nonframeshift insertion非移码插入、nonframeshift deletion非移码缺失
AAChange.refGeneMTHFR:NM_001330358:exon5:c.C788T:p.A222V参考refGene数据库,对氨基酸改变情况进行注释,当该突变位于exonic时才有结果。按照每个转录本进行注释,例如:SAMD11, NM_152486, exon3,c.232C>T, p.H78Y。其中,SAMD11表示该变异所在的基因名称,NM_152486表示该变异所在的转录本ID,exon3表示该变异位于转录本的第三个外显子上,c.232C>T表示cDNA的232位置上发生C到T的突变,p.H78Y表示该蛋白序列的78位置上发生H到Y的氨基酸改变
Func.HGVSmissense_variant另一款软件根据HGVS法则注释出来的突变功能类型结果
AAChange.HGVSMTHFR:NM_005957.4:5/12:c.665C>T:p.(Ala222Val)另一款软件根据HGVS法则注释出来的氨基酸改变情况,格式为 基因:转录本:外显子号/外显子总数:cDNA坐标:氨基酸变化
cytoBand1p36.22突变所处的染色体区段名称
InterVar_automatedBenign突变在InterVar数据库中的有害性注释
ACMG(missense only)PM1,BA1,BS1,BP6突变的ACMG评级结果
wgRna.基于miRBase和snoRNABase,对突变的microRNA和snoRNA进行注释,给出microRNA和snoRNA 的基因名称
genomicSuperDups.检测该变异是否位于重复片段中,位于重复片段中的突变往往是由于比对错误造成的,很有可能是假阳性位点。例如Score=0.994427;Name=chr16:21740424,Score表示两个相似片段的一致性(0-1),Name表示该片段和chr16:21740424位置的片段相似
rmsk.RepeatMasker,检测该变异是否位于散在重复序列区域或者低复杂序列区域中,包括SINE、ALU、LINE、LTR、简单重复序列、卫星序列等多种类型。
gwasCatalogName=Homocysteine levels在以往的GWAS研究中被报道,位于该区间中的结果,表示该区间可能与哪些疾病有关联
avsnp156rs1801133该变异在dbSNP156版本中的ID
1000g2015aug_all0.245407千人基因组计划2015年8月版中该突变的等位基因频率,包含千人基因组中所有人群
1000g2015aug_eas0.2956千人基因组计划2015年8月版中该突变的等位基因频率,仅包含东亚的人群
esp6500siv2_all0.2706美国国家心肺和血液研究所外显子测序计划(NHLBI-ESP),包含6500个外显子,人群包括非洲裔美国人和欧裔美国人
ExAC_ALL0.3037Exome Aggregation Consortium数据库,包括60000多个无亲缘关系个体的数据,包含所有人群
ExAC_EAS0.3052Exome Aggregation Consortium数据库,包括60000多个无亲缘关系个体的数据,仅包含东亚人群
gnomAD_exome_AF_raw0.3143Genome Aggregation Database数据库,是由各国研究者联合发展起来的基因组突变数据库,其中数据集包括123136个全外显子数据和15496个全基因组数据,来源于各种疾病研究项目和大型人群测序项目。本数据库包含所有人群的全外显子数据
gnomAD_exome_AF_eas0.2937Genome Aggregation Database数据库,是由各国研究者联合发展起来的基因组突变数据库,其中数据集包括123136个全外显子数据和15496个全基因组数据,来源于各种疾病研究项目和大型人群测序项目。本数据库包含东亚人群的全外显子数据
gnomAD_genome_AF_raw0.2572Genome Aggregation Database数据库,是由各国研究者联合发展起来的基因组突变数据库,其中数据集包括123136个全外显子数据和15496个全基因组数据,来源于各种疾病研究项目和大型人群测序项目。本数据库包含所有人群的全基因组数据
gnomAD_genome_AF_eas0.2866Genome Aggregation Database数据库,是由各国研究者联合发展起来的基因组突变数据库,其中数据集包括123136个全外显子数据和15496个全基因组数据,来源于各种疾病研究项目和大型人群测序项目。本数据库包含东亚人群的全基因组数据
ICGC_IdMU580088ICGC(International Cancer Genome Consortium)国际癌症基因组联盟打造的癌症数据库,包括上万个癌症数据,本结果显示突变在该库中的id
ICGC_OccurrenceLAML-KR|1|107|0.00935,LUSC-KR|1|66|0.01515ICGC(International Cancer Genome Consortium)国际癌症基因组联盟打造的癌症数据库,包括上万个癌症数据,本结果显示突变在该库中的详细情况。例如LUSC-KR|1|66|0.01515代表了“癌症名称缩写-国家|发生突变的样本数|该项目研究的样本数|突变频率”
iGeneTechDB0.38902743本公司自行整理的全外数据库,包含4600多个正常人全外数据
CLNALLELEID18559该变异在临床疾病数据库clinvar中的注释结果,CLNALLELEID 代表数据库中的编号
CLNDNNeoplasm_of_stomach|Gastrointestinal_stroma_tumor该变异在临床疾病数据库clinvar中的注释结果,CLNDN 代表变异位点相关的疾病名称
CLNDISDBHuman_Phenotype_Ontology:HP:0006753,MeSH:D013274该变异在临床疾病数据库clinvar中的注释结果,CLNDSDB 代表疾病关联信息的数据库来源
CLNREVSTATreviewed_by_expert_panel该变异在临床疾病数据库clinvar中的注释结果
CLNSIGdrug_response该变异在临床疾病数据库clinvar中的注释结果,CLINSIG 代表变异位点在临床意义
cosmic99.Catalogue Of Somatic Mutations In Cancer,世界上研究人体肿瘤体细胞突变最大和最全面的数据库,结果分为id和疾病,用逗号分隔
HGMDCardiovascular_disease_association_with人类基因突变数据库,本结果给出该位点在数据库涉及的疾病名称
dbscSNV_ADA_SCORE.预测突变对可变剪切影响力的数据库,本结果采用AdaBoost算法进行预测,数值表示有影响的概率,例如0.5代表有50%的可能性对可变剪切有影响
dbscSNV_RF_SCORE.预测突变对可变剪切影响力的数据库,本结果采用Random Forest算法进行预测,数值表示有影响的概率,例如0.5代表有50%的可能性对可变剪切有影响
SIFT_score0.002SIFT分值,表示该变异对蛋白质的影响,分值越小越有害
SIFT_predDSIFT预测结果说明,取值为 T 或者 D。当该变异同时影响多个蛋白序列时,对每条蛋白序列有一个 SIFT 值,取其中最小值。D: Deleterious (分值<=0.05), T: tolerated (分值>0.05)
SIFT4G_score0.002SIFT的另一个版本,具有更高的准确度及更快的运行速度,同样表示该变异对蛋白质的影响,分值越小越有害
SIFT4G_predDSIFT4G预测结果说明,取值为 T 或者 D。当该变异同时影响多个蛋白序列时,对每条蛋白序列有一个 SIFT 值,取其中最小值。D: Deleterious (分值<=0.05), T: tolerated (分值>0.05)
Polyphen2_HDIV_score0.998PolyPhen2基于 HumanDiv 数据库预测该变异对蛋白序列影响的分值,适用复杂表型中罕见等位基因位点突变的诊断,数值越大越有害,表明该 SNP 导致蛋白结构或功能改变的可能性越大
Polyphen2_HDIV_predDPolyphen2 HDIV数据库的预测结果说明,取值为D或P或B,代表的意义如下D: Probably damaging (分值>=0.957), P: possibly damaging (0.453<=分值<=0.956), B: benign (分值<=0.452)
Polyphen2_HVAR_score0.941PolyPhen2基于 HumanVar 数据库预测该变异对蛋白序列影响的分值,适用于符合孟德尔遗传定律疾病的诊断,数值越大越有害,表明该 SNP 导致蛋白结构或功能改变的可能性越大
Polyphen2_HVAR_predDPolyphen2 HVAR数据库的预测结果说明,取值为D或P或B,代表的意义如下D: Probably damaging (分值>=0.909), P: possibly damaging (0.447<=分值<=0.909), B: benign (分值<=0.446)
LRT_score0LRT软件的预测结果,分值越小越有害
LRT_predDLRT的预测结果说明,取值为D或N或U,代表的意义如下D: Deleterious, N:Neutral, U:Unknown
MutationTaster_score0MutationTaster软件的预测结果,分值越大结果越可靠
MutationTaster_predPMutationTaster的预测结果说明,取值为A或D或N或P,代表的意义如下A: Disease_causing_automatic(有害), D: Disease_causing(有害), N: Polymorphism(无害), P: Polymorphism_automatic(无害)
MutationAssessor_score0MutationAssessor软件的预测结果,数值越大越有害
MutationAssessor_predNMutationAssessor的预测结果说明,取值为H或M或L或N,代表的意义如下H: High, M: Medium, L: Low, N: Neutral
FATHMM_score-4.03FATHMM软件的预测结果,分值越小越有害
FATHMM_predDFATHMM的预测结果说明,取值为D或T,代表的意义如下D: Deleterious, T: Tolerated
PROVEAN_score-0.82PROVEAN蛋白效应预测软件的预测结果,一般分数小于-2.5就可以认为是有害的
PROVEAN_predNPROVEAN的预测结果说明,取值为D或N,代表的意义如下D: Deleterious, T: Neutral
VEST4_score0.042Variant effect scoring tool分值。该软件采用机器学习随机森林算法,分值越高危害越大
MetaSVM_score-0.934MetaSVM结合SIFT, PolyPhen 和 MutationAssessor 的预测分值,训练SVM模型来预测突变有害性,分值越高危害越大
MetaSVM_predTMetaSVM的预测结果说明,取值为D或T,代表的意义如下D: Deleterious, T: Tolerated
REVEL_score0.05REVEL评分, 分类建议:0.5以下为不致病,0.85以上为致病
ClinPred_score0.002ClinPred主要侧重于对非同义突变的危害性进行预测,分值越高危害越大
ClinPred_predTClinPred的预测结果说明,取值为D或T,代表的意义如下D: Deleterious, T: Tolerated
CADD_raw4.914CADD 预测初始分值。CADD是一种对 SNV、InDel 的有害性进行打分的工具,它整合多种信息来注释变异位点的功能;不仅预测编码区变异(包括同义突变和非同义突变的影响)的功能影响,还预测非编码区变异的功能影响。对于 SNP,仅对 CADD 分值排名在前 10%的SNP 给出分值,"."表示 CADD 分值排名不在前 10%
CADD_phred25CADD 转换后的分值,没有分值,即为"."时,表示 CADD_Phred 值小于10。CADD_Phred 分值中,10 表示 score 排名在前 10%,20 表示前 1%,30 表示前 0.1%,因此,分值要求越低,能保留下来的位点越多。对于 SNP,CADD 作者建议 CADD_Phred 分值>15,文章中通常用 10或 15;InDel 没有建议值
fathmm-MKL_coding_score0.98fathmm-MKL软件的预测结果,分值越小越有害
fathmm-MKL_coding_predDfathmm-MKL的预测结果说明,取值为D或T,代表的意义如下D: Deleterious, T: Tolerated
GERP++_NR5.01GERP++软件NR算法的预测结果,分值越大越有害
GERP++_RS5.08GERP++软件RS算法的预测结果,分值越大越有害
OMIM_InheritanceAutosomal dominant, Autosomal recessive人类孟德尔遗传病数据库注释,本结果给出遗传模式
OMIM_DiseaseHomocystinuria due to MTHFR deficiency, 236250 (3), Autosomal recessive;人类孟德尔遗传病数据库注释,对应关系为基因-疾病,多种疾病以分号隔开,每种疾病的注释格式为:疾病,编号,遗传模式
HGMD_DiseasePlacental vasculopathies, increased risk|Partial hypoplasia of inferior vena cava 人类基因突变数据库,对应关系为基因-疾病,多种疾病以竖线隔开
GOGO:0001666 response to hypoxiaGO注释结果,如果一个基因注释上多个GO,GO间以分号隔开,每个GO的注释格式为:GO:GO_ID空格GO_term
KEGG_pathway00670 One carbon pool by folate,01523 Antifolate resistanceKEGG_pathway注释结果,如果一个基因注释上多个pathway,pathway间以分号隔开,每个pathway的注释格式为:kegg_ID空格kegg_term

5.2 InDel检测和注释

InDel是小片段插入和缺失(Insertion 和 Deletion)的简称,其长度通常在50 bp以下。其中编码区发生的移码突变对蛋白的影响非常大。找出InDel以后我们同样会对结果进行详细的注释,注释内容等同于SNP,见表5突变结果表。

特别声明:哺乳动物的X染色体和Y染色体具有共同的进化起源,并保留着高度相似的序列区域。所以就算是女性,也可能存在一些低质量的reads比对到Y染色体,并检测到一些突变(见参考文献10)。

5.3 突变数量统计

不同样本在同一区间的突变数量级大体接近,且绝大数的突变在dbsnp数据库中有记录,因此我们可以通过突变数量及其在dbsnp中的占比判断突变结果是否正常。统计结果及说明如下,通常注释到dbsnp中的SNP在95%左右,InDel在90%左右:

表6 SNP数量统计表
指标 数值 说明
Total snp(filter offtarget)1914检测出的SNP总数
ontarget506(26.44%)位于目标区域内的结果
flank1501408(73.56%)位于目标区域周围150bp的结果
offtarget0(0.00%)位于目标区域外的结果
dbsnp1914(100.00%)注释到dbsnp中的突变数
1000g_EAS672(35.11%)注释到千人基因组东亚人群上的突变数
ExAC_EAS470(24.56%)注释到ExAC东亚人群上的突变数
gnomAD_exome_EAS1720(89.86%)注释到gnomAD东亚人群全外显子数据上的突变数
gnomAD_genome_EAS1338(69.91%)注释到gnomAD东亚人群全基因组数据上的突变数
1/1(Hom)322(16.82%)纯合突变数目
0/1(Het)1592(83.18%)杂合突变数目
1/2(Het_m)0(0.00%)多基因型杂合突变数目,在InDel发生区域较常见
exonic480(25.08%)注释到外显子区域的突变数
splicing12(0.63%)注释到可变剪切区域的突变数
UTR371(3.71%)注释到UTR3的突变数
UTR566(3.45%)注释到UTR5的突变数
intronic1251(65.36%)注释到内含子区域的突变数
intergenic0(0.00%)注释到基因间的突变数
upstream8(0.42%)注释到基因上游的突变数
downstream1(0.05%)注释到基因下游的突变数
ncRNA_exonic5(0.26%)注释到非编码外显子的突变数
ncRNA_splicing0(0.00%)注释到非编码可变剪切的突变数
ncRNA_intronic20(1.04%)注释到非编码内含子的突变数
synonymous SNV207(10.82%)同义突变数
nonsynonymous SNV266(13.90%)非同义突变数
stopgain7(0.37%)突变导致所在的密码子变为终止密码子
stoploss0(0.00%)突变导致所在的终止密码子变为非终止密码子
unknown0(0.00%)因为数据库中基因结构定义错误等问题导致无法注释的突变
表7 InDel数量统计表
指标 数值 说明
Total indel(filter offtarget)747检测出的InDel总数
ontarget118(15.80%)位于目标区域内的结果
flank150629(84.20%)位于目标区域周围150bp的结果
offtarget0(0.00%)位于目标区域外的结果
dbsnp747(100.00%)注释到dbsnp中的突变数
1000g_EAS124(16.60%)注释到千人基因组东亚人群上的突变数
ExAC_EAS315(42.17%)注释到ExAC东亚人群上的突变数
gnomAD_exome_EAS720(96.39%)注释到gnomAD东亚人群全外显子数据上的突变数
gnomAD_genome_EAS649(86.88%)注释到gnomAD东亚人群全基因组数据上的突变数
1/1(Hom)16(2.14%)纯合突变数目
0/1(Het)731(97.86%)杂合突变数目
1/2(Het_m)0(0.00%)多基因型杂合突变数目,在InDel发生区域较常见
exonic105(14.06%)注释到外显子区域的突变数
splicing6(0.80%)注释到可变剪切区域的突变数
UTR330(4.02%)注释到UTR3的突变数
UTR515(2.01%)注释到UTR5的突变数
intronic584(78.18%)注释到内含子区域的突变数
intergenic0(0.00%)注释到基因间的突变数
upstream2(0.27%)注释到基因上游的突变数
downstream1(0.13%)注释到基因下游的突变数
ncRNA_exonic1(0.13%)注释到非编码外显子的突变数
ncRNA_splicing0(0.00%)注释到非编码可变剪切的突变数
ncRNA_intronic3(0.40%)注释到非编码内含子的突变数
frameshift insertion9(1.20%)移码突变数(插入)
frameshift deletion28(3.75%)移码突变数(缺失)
nonframeshift insertion3(0.40%)非移码突变数(插入)
nonframeshift deletion63(8.43%)非移码突变数(缺失)
stopgain2(0.27%)突变导致所在的密码子变为终止密码子
stoploss0(0.00%)突变导致所在的终止密码子变为非终止密码子
unknown0(0.00%)因为数据库中基因结构定义错误等问题导致无法注释的突变

5.4 突变类型分布

同一目标区域中,突变类型分布理论上是接近的,如SNP的转换和颠换的总数,InDel的长度分布趋势等。我们通过图片直观地展示分布类型。

SNP突变类型分布图
图8 SNP突变类型分布图
InDel突变类型分布图
图9 InDel突变长度分布图

5.5 CNV检测与注释(需基线库)

CNV(Copy number variations)是基因拷贝数变异,异常的DNA拷贝数变化是许多人类疾病(如癌症、遗传性疾病、心血管疾病)的一种重要分子机制。作为疾病的一项生物标志,染色体水平的缺失、扩增等变化已成为许多疾病研究的热点。针对单样本全外显子,我们采用公司自行构建的数据库做为正常参照进行CNV检测。找出结果以后我们对CNV进行注释(没有构建基线库则无此结果),注释内容及说明如下:

表8 CNV结果表
指标 示例 说明
chromosomechr9染色体号
start84543421CNV区间起始位置
end84562625CNV区间终止位置
geneSPATA31D4,SPATA31D3区间覆盖到的基因,覆盖到多个基因则用逗号隔开
log2-1.27188平均深度的log值,经过了与数据库的均一化处理及异常数据筛除
cn0CNV判断的结果,N代表N个拷贝(人类中正常为2),判断过程受log2值影响,大致的范围为0(log2值小于-1.1),1(log2值-1.1至-0.4),2(log2值-0.4至0.3),3(log2值0.3至0.7)
depth87.1687区间平均深度
probes21该区间包含的bin数量(默认300bp为一个bin)
weight12.3587衡量该区间可靠性的值,数值越大越可靠,一般与区间大小成正比
length19204区间长度
Func.refGeneexonic参考refGene数据库,对CNV所在的区域进行注释,如果该区域覆盖到了exonic区域,则会优先注释为exonic
Gene.refGeneSPATA31D4参考refGene数据库,列出CNV所交集的基因,如果不在基因上,则列出最近的基因名(有可能为两个)
cytoBand9q21.32CNV覆盖到的染色体区段名称,多个区段用横线连接
wgRna.基于miRBase和snoRNABase,对突变的microRNA和snoRNA进行注释,给出microRNA和snoRNA 的基因名称
genomicSuperDupsScore=0.986134;Name=chr9:84536419,chr9:84551491检测区间是否包含与基因组其它位置相似片段
gwasCatalog.在以往的GWAS研究中被报道,位于该区间中的结果,表示该区间可能与哪些疾病有关联
CNVD.CNVD数据库,包含和人类疾病相关的212277条数据。与该库中overlap区域超过60%的CNV会被注释出来

5.6 SV检测与注释

SV(Structural variation)指的是结构变异,主要类型有缺失、重复、倒位、易位。人类的许多遗传病是由染色体结构改变引起的,在癌症中也已经证实部分癌症的产生和结构变异导致的基因融合有关。在单样本分析中我们采用检测断点周围序列特征的方式来检测结构变异。检测结果说明如下:

表9 SV结果表
指标 示例 说明
Chr1chrM断点1染色体号
Breakpoint110断点1位置
Gene1.断点1涉及到的基因,格式为基因名:转录本号:外显子/内含子编号,没有在基因上则以点号表示
Chr2chrM断点2染色体号
Breakpoint216569断点2位置
Gene2.断点2涉及到的基因,格式同Gene1
Length16559两断点之间的距离,如果断点在不同的染色体则以点号表示
Sv_typeDUP结构变异类型,分为DUP(片段重复),DEL(片段缺失),INV(片段倒位),BND(片段易位)四种
Vcf_altDUPVCF结果文件中显示的变异类型,其中DUP、DEL、INV显示结果与Sv_type一致,BND显示结果为另一侧断点的位置及断点组合顺序
Gene_type1/1该位置的基因型
Ref/Alt0/23支持参考基因组的深度/支持突变的深度
Depth23该位置的总深度
Quality649.93质量值,值越高代表可信度越高
Fre1该位置的突变频率
Vcf_infoSVTYPE=DUP;STRANDS=-+:373;SVLEN=16559;END=16569...结构变异在VCF中的详细信息

六、高级分析结果

6.1 突变位点意义评价

在标准分析中会检测到较多的突变位点(通常几千到几万),其中大部分突变位点是良性位点或者与疾病无关的位点,为了方便用户从海量位点出筛选出有意义的位点,我们按照以下判定条件对突变的意义进行了分级评价。优先级分为high、likeyhigh、medium、low 四级,我们将优先级不为low的结果存为一个新的文件(significant.snp_indel.xls),初步缩小了挑选范围。

high: 位于基因的exonic区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值小于0.01,造成氨基酸改变,且至少一款预测软件预测为有害;

likeyhigh: 位于基因的exonic或splicing区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值至少有一个小于0.01,且至少一款预测软件预测为有害;

medium: 位于基因的exonic或splicing区域,在东亚人群中(1000g,ExAC,gnomAD)的MAF值至少有一个小于0.01;

low: 不属于上述3种的其它结果。

6.2 突变位点质量值评价

标准分析检测到的结果,根据其深度、突变频率、质量值等不同,结果有真有假,通常我们可以采用IGV软件查看bam文件或者一代验证的方法确认其真实性。然而大量的突变位点结果会造成时间上的浪费,为了节约时间,我们按照以下判定条件对突变的质量值进行了评价。评价级别分为highquality、mediumquality、lowquality三级,建议用户挑选中高质量的位点进行后续一代验证。

highquality:深度大于等于50层,基因型质量值大于30,频率偏差不超过0.1(杂合0.4-0.6,纯合0.9-1.0)的结果;

mediumquality:深度大约20层,基因型质量值大于20,频率偏差不超过0.2(杂合0.3-0.7,纯合0.8-1.0)的结果;

lowquality:不属于上述2种的其它结果,以及位于repeat区域的结果(genomicSuperDups数据库注释结果不为“.”)。

6.3 突变位置筛选

目标区域捕获除了区域内的结果外,也会产出部分非目标区域的结果。为了保证结果的准确性,我们只保留目标区域内及目标区域周围150bp的突变结果,如无特殊需要,建议用户筛选出目标区域内的结果用于后续分析。

ontarget:结果位于目标区域内;

flank150:结果位于目标区域周围150bp内;

offtarget:目标区域外的结果,该类结果会被直接丢弃。

6.4 突变位点有害性分类

2015年,美国权威机构——美国医学遗传学与基因组学学会(ACMG)编写和发布了《ACMG遗传变异分类标准与指南》。该指南将变异位点的致病、良性证据列为具体的28条评判标准。首先将证据按类型分类(如人群数据、计算预测数据、功能数据等),并将证据的支持度分为几类(支持,中等,强,非常强以及独立);然后使用“标准组合”的形式来评估致病性。不同组合将产生五个类别的致病性分类:致病,可能致病,临床意义不明,可能良性,良性,是数据解读中的金标准,其判断致病或者良性的标准如下:

图10 ACMG标准与指南

我们基于ANNOVAR整理的intervar数据库对突变位点有害性进行分类,该数据库目前只针对错义突变进行评级,且随着新的证据出现,其评级可能产生变化。

6.5 突变位点有害性预测

对于尚未报道或者缺乏有效证据支持的SNP结果,我们提供了多款软件预测突变的有害性或者保守性,例如SIFT、Polyphen2、LRT、MutationTaster、CADD等,建议选择至少60%以上软件预测为有害的突变结果用于后续分析。

附1、数据库列表

数据库 版本或下载日期
参考基因组 hg19/GRCh38 2009/2013
突变注释数据库,
大多数为ANNOVAR整理
refGene 202110
cytoBand202204
genomicSuperDups202204
rmsk202204
wgRna202204
gwasCatalog202204
dbnsfp42a202107
intervar201801
cosmic99202403
clinvar202403
dbscsnv11201512
icgc28202101
gnomad211_genome201904
gnomad211_exome201904
gnomad40_genome202311
gnomad40_exome202311
1000g201508
ExAC201511
esp6500201412
avsnp156202404
HGMD201412
OMIM202403
iGeneTechDB201901
GO201801
KEGG_pathway201801

附2、软件及工具列表

功能 工具
质量控制及评估Trimmomatic-0.38
比对Bwa-0.7.12
去重samblaster 0.1.22
矫正GATK-3.8-0
SNP,InDel检测GATK-3.8-0、Samtools-1.9、Varscan-v2.4.3
CNV检测Cnvkit-0.9.9、ADTEx-v.2.0
SV检测Lumpy-0.2.13
突变注释Annovar-201707、snpEff-4.3i

附3、参考文献

Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754-60 (2009).
Faust, G.G. & Hall, I.M. SAMBLASTER: fast duplicate marking and structural variant read extraction. Bioinformatics 30, 2503-5 (2014).
McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20, 1297-303 (2010).
Koboldt, D.C. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res 22, 568-76 (2012).
Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078-9 (2009).
Talevich, E., Shain, A.H., Botton, T. & Bastian, B.C. CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing. PLoS Comput Biol 12, e1004873 (2016).
Amarasinghe, K.C. et al. Inferring copy number and genotype in tumour exome data. BMC Genomics 15, 732 (2014).
Layer, R.M., Chiang, C., Quinlan, A.R. & Hall, I.M. LUMPY: a probabilistic framework for structural variant discovery. Genome Biol 15, R84 (2014).
Wang, K., Li, M. & Hakonarson, H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38, e164 (2010).
Webster, T.H. et al. Identifying, understanding, and correcting technical artifacts on the sex chromosomes in next-generation sequencing data. GigaScience 8 (2019).