微科盟_人类外显子重测序分析_报告(肿瘤版)


1. 项目背景

1.1 癌症基因组研究简介

肿瘤(Tumor)是机体在各种致癌因素作用下,局部组织的细胞在基因水平上失去对 其生长的正常调控,导致其克隆性异常增生而形成的病变。癌症基因组一般会发生独特且 不可预测的多个点突变、插入缺失、易位、融合等事件,对于这些突变的研究可分为两个 层面:Germline Mutation(胚系突变)和Somatic Mutation(体细胞突变)。 前者侧重于可遗传的、个体背景中所有细胞携带的突变,主要用于遗传易感性及药物基因 组学的相关研究;后者侧重于肿瘤细胞特有、正常细胞没有的一类突变,通过寻找肿瘤与 正常细胞间突变信息的差异,进一步研究癌症发生发展的机制。

外显子组测序(Whole Exome Sequencing,WES)是利用探针杂交富集外显子区域的 DNA序列,然后通过高通量测序,主要识别和研究与癌症相关的编码区相关突变的技术手 段。相较于全基因组测序,由于外显子只占人类基因组约1%(约30Mb),因此更容易做 到高深度测序,检测到更多低频和罕见变异,同时也能降低费用,缩短周期。

1.2 样品信息

对某癌种同一亚型的多个患者成对取样,即取癌组织样本作为 T,取同一个病人的癌 旁组织或血液样本作为N。

表 1-1. 样品信息一览表

2. 建库测序流程

样品检测、建库、上机测序每个环节都会影响测序数据的数量和质量,而数据质量又 会直接影响后续信息分析的结果,因此获得高质量数据是保证生物信息分析正确、 全面、可信的前提。为了从源头上保证测序数据的准确性、可靠性,对样品检测、 建库、测序每一个生产步骤都严格把控,从根本上确保了高质量数据的产出。 流程图如下:

图2-1.外显子组建库测序流程图

2.1 DNA样品检测

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

1)琼脂糖凝胶电泳分析DNA降解程度以及是否存在杂带、RNA及蛋白污染;

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

其中DNA浓度≥20 ng/µL,总量0.6µg以上的DNA样品被用来建库。

2.2 文库构建

将基因组DNA经Covaris破碎仪随机打断成长度为180-280 bp的片段,建库和捕获 实验采用Agilent SureSelect Human All Exon V5/V6 试剂盒,末端修复、 磷酸化以及加A尾后,片段两端分别连接上接头制备 DNA 文库。带有特异 index 的文库与多达543,872 个生物素标记的探针进行液相杂交,再使用带链霉素的磁珠 将20,965个基因的334,378个外显子捕获下来,经PCR线性扩增后进行文库质检, 合格即可进行测序。

图2-2.杂交捕获原理图(Agilent SureSelect)

2.3 库检

文库构建完成后,先使用Qubit 2.0进行初步定量,随后使用Agilent 2100对文库的 Insert size进行检测,Insert size符合预期后,使用Q-PCR方法对文库的有效浓 度(3 nmol/L)进行准确定量,以保证文库质量。

2.4 上机测序

库检合格后,根据文库的有效浓度及数据产出需求进行Illumina HiSeq PE150测序。 PE150(Pair End 150 bp)指高通量双端测序,每端各测150 bp。 在构建的小片段文库中,Insert DNA,即插入片段是高通量测序直接测序的单位。 双端测序是将每条插入片段的两端进行测序的方法,由于插入片段的长度分布已知, 双端测序时不仅可以知道片段两端的序列,也能知道这两段序列之间的长度, 从而便于后续比对分析。

3 生物信息分析流程

基本分析A版部分主要用于获取癌症样本基因组的突变信息。 首先,对测序原始数据(Raw data)进行质控,得到高质量的Clean data; 然后,将Clean data与人参考基因组序列进行比对分析,获得Bam文件; 最后基于Bam进行Somatic Mutation的检出与注释,从而得到癌症机制研究 的基本突变信息结果,基本分析A版报告主要展示Somatic相关结果。

3.1 测序数据质量评估

测序数据的产生经过了DNA提取、建库、测序多个步骤, 这些步骤会产生少许低质量或者无效数据,比如建库阶段会存在建库长度的偏差, 测序阶段会出现测序错误的情况。所以需要对下机得到的原始数据进行质量控制, 保证后续分析的准确性。

3.1.1 原始序列数据

高通量测序(如Illumina测序平台)测序得到的原始图像数据文件经碱基识别 (Base Calling)分析转化为原始测序序列(Sequenced Reads), 我们称之为Raw Data或Raw Reads,结果以FASTQ(简称为fq) 文件格式存储,其中包含测序序列(reads)的序列信息以及其对应的测序质量信息。

FASTQ格式文件中每条Read对应四行内容:

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
   GCTCTTTGCCCTTCTCGTCGAAAATTGTCTCCTCATTCGAAACTTCTCT
   +
   @@CFFFDEHHHHFIJJJ@FHGIIIEHIIJBHHHIJJEGIIJJIGHIGHC

其中第一行以“@”开头,随后为Illumina测序标识符(Sequence Identifiers) 和描述文字(选择性部分);第二行是碱基序列;第三行以“+”开头, 随后为Illumina测序标识符(选择性部分);第四行是对应序列的测序质量。

Illumina测序标识符详细信息如下:

第四行中每个字符对应的ASCII值减去33,即为对应第二行碱基的测序质量值。 如果测序错误率用e表示,Illumina的碱基质量值用Qphred表示,则有下列关系:

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

此公式可说明,质量值越大测序错误率越低,准确性越高。

Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系:

3.1.2 测序数据过滤

测序获得的原始数据中包含少量带有测序接头或测序质量较低的reads, 为保证数据分析的质量及可靠性,需要对原始数据进行过滤。 本分析使用Trimmomatic软件对测序数据进行过滤, 过滤后各部分reads所占比例均在饼图中呈现。 分析结果在文件夹/01dataFilter中呈现,此处以一个样本结果作为示例。

图3.1 原始数据(Raw data)过滤结果

注:
   Both Surviving:该部分为过滤后read1和read2均被保留的数据
   Forward Only Surviving:该部分为仅在read1中被保留的数据
   Reverse Only Surviving:该部分为仅在read2中被保留的数据
   Dropped:该部分为因接头或数据质量等原因被丢弃的数据

3.1.3 测序错误率分布检查

每个碱基的测序Phred值(Phred score,Qphred)是由测序错误率通过公式一 转化得到的,而测序错误率是在碱基识别(Base Calling)过程中通过一种判别 发生错误概率的模型计算得到的。对应关系如下表所显示:

表3.1 Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系:

测序错误率与碱基质量有关,受测序仪本身、测序试剂、样品等多个因素共同影响。 对于Illumina高通量测序平台,测序错误率分布具有两个特点:
   (1)测序错误率会随着测序的进行而升高,这是由于测序过程中 荧光标记的不完全切割等因素引起荧光信号衰减,因而导致错误率升高;
   (2)每个Read前几个碱基的位置也会有较高的测序错误率, 这是由于边合成边测序过程初始阶段,测序仪荧光感光元件对焦速度较慢, 获取的荧光图像质量较低,导致碱基识别错误率较高。

该部分分析结果在文件夹/02fastqcPic中呈现,此处以其中一个样本结果作为示例:

图3.2 测序错误率分布

3.1.4 GC含量分布

核苷酸序列中鸟嘌呤(G)和胞嘧啶(C)所占的比例称为GC含量。 GC含量在物种间存在一定特异性,但由于反转录过程中所使用的6bp随机引物, 会引起前几位碱基在核苷酸组成上有一定偏好性,产生正常波动,随后则趋于稳定。 对于NEB普通建库方法,由于序列的随机性打断和双链互补等原则,理论上测序读段 在每个位置的GC及AT含量应分别相等,且在整个测序过程基本稳定不变,呈水平线。 而对于链特异性建库而言,由于只保留了单链信息,可能会出现AT分离或GC分离现象。

该部分分析结果在文件夹/02fastqcPic中呈现,此处以其中一个样本结果作为示例:

图3.3 样本GC含量分布

3.1.5 测序数据质量分布

测序数据的质量主要分布在Q30(≥80%)以上,这样才能保证后续分析的正常进行, 根据测序技术的特点,测序片段末端的碱基质量一般会比前端的低。

该部分分析结果在文件夹/02fastqcPic中呈现,此处以其中一个样本结果作为示例:

图3.4 样本测序质量分布

3.2 测序深度、覆盖度统计

有效测序数据通过BWA比对到参考基因组,得到BAM格式的最初的比对结果。 BAM文件再进行标记重复等处理,从而得到BAM格式的最终比对结果, 再利用比对到参考基因组上的有效数据进行覆盖度的统计。 通常,人类样本的测序Reads能达到98%以上的比对率, 测序深度在10×以上Reads覆盖的位点检测出的SNV比较可信。

注:该部分结果在文件夹03sequence_depth中。

图3.5 测序深度

该图为测序深度分布图,曲线代表不同测序深度对应的比例。 该部分分析结果在文件夹03sequence_depth中呈现, 此处以一个样本的结果进行示例。

图3.6 累计测序深度

该图为累计测序深度图,即不同测序深度的累计碱基比例。横坐标表示测序深度, 纵坐标表示测序深度超过对应测序深度的碱基比例。 该部分分析结果在文件夹03sequence_depth中呈现, 此处以一个样本的结果进行示例。

图3.7 每条染色体的测序深度(左侧坐标)和覆盖率(右侧坐标)

横坐标表示染色体编号,左侧纵坐标表示平均覆盖深度, 右侧纵坐标表示覆盖率。对每条染色体计算覆盖深度时,计算公式为: 每条染色体的测序数据量/每条染色体上外显子区域的总长度。 计算覆盖率时,公式为:每条染色体被覆盖的总长度/每条染色体上外显子区域的总长度。 该部分结果在文件夹03sequence_depth中呈现, 此处以一个样本的结果进行示例。

4 体细胞变异检测

体细胞突变(Somatic mutation)是指除胚系细胞以外的体细胞所发生的突变, 是发生在正常机体细胞中的突变,如发生在皮肤和器官。 体细胞突变既不遗传自亲本,也不会传递给后代, 却可以引起当代某些细胞的遗传结构发生改变。体细胞突变, 特别是其中的驱动突变(Driver mutation) 对解释肿瘤的发生和发展具有非常重要的意义, 另一方面肿瘤耐药性的产生也与体细胞突变有关。 因此,关注体细胞突变是肿瘤基因组研究的重心, 也是肿瘤基因组研究区别于疾病研究的一个特性。 在信息分析时,通过将肿瘤组织与正常组织相比, 过滤掉正常组织中的胚系细胞突变(Germline mutation), 只保留下肿瘤细胞携带的体细胞突变。

4.1 Somatic SNV检测结果

SNV全称Single nucleotide variant,即单核苷酸突变, 指在基因组上由单个核苷酸的替换所引起的变异。我们主要使用GATK muTect2软件 寻找Somatic SNV。体细胞突变检测结果如下表所示 (具体每个SNV的信息详见结果注释文件):

注:该部分结果在文件夹04SNP中。

表4.1 SNV识别结果(部分展示)

注:
   CHROM:染色体
   POS:SNV在基因组位置
   ID:SNV ID
   REF:参考等位
   ALT:突变等位
   QUAL:质量值
   FILTER:筛选信息
   INFO:SNV信息,该部分信息较多,每项内容参考文件夹中文件的注释信息
   FORMAT:格式信息
      GT:Genotype;
      AD:Allelic depths for the ref and alt alleles in the order listed;
      AF:Allele fractions of alternate alleles in the tumor;
      DP:Approximate read depth; some reads may have been filtered;
      F1R2:Count of reads in F1R2 pair orientation supporting each allele;
      F2R1:Count of reads in F2R1 pair orientation supporting each allele;
      SB:Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.
   sample name:样本sample的GT信息

4.2 Somatic InDel检测结果

InDel全称Insertion and Deletion,指小片段的插入和缺失。 编码区或剪接位点处发生的插入缺失都可能会改变蛋白的翻译。 移码变异,其插入或缺失的碱基串的长度为3的非整数倍,因此可能导致整个阅读框的改变; 与非移码变异比较,移码突变对基因功能的影响更大,同时也受到更大的筛选压力。 我们利用Strelka软件(Saunders CT et al. 2012)检测Somatic InDel信息, 对得到的InDel结果统计如下(具体每个InDel变异的信息详见结果注释文件):

注:该部分结果在05indel中。

表4.2 Indel识别结果(部分展示)

注:
   #CHROM:染色体号
   POS:Indel的基因组位置
   ID:Indel ID
   REF:Indel参考等位
   ALT:突变等位
   QUAL:Indel质量
   FILTER:筛选信息
   INFO:Indel信息,该部分内容较多,可在文件夹下结果文件的注释部分查看具体含义
   FORMAT:格式信息
      DP:Read depth for tier1
      DP2:Read depth for tier2
      TAR:Reads strongly supporting alternate allele for tiers 1,2
      TIR:Reads strongly supporting indel allele for tiers 1,2
      TOR:Other reads (weak support or insufficient indel breakpoint overlap) for tiers 1,2
      DP50:Average tier1 read depth within 50 bases
      FDP50:Average tier1 number of basecalls filtered from original read depth within 50 bases
      SUBDP50:Average number of reads below tier1 mapping quality threshold aligned across sites within 50 bases
      BCN50:Fraction of filtered reads within 50 bases of the indel.
   NORMAL:正常样本的GT信息
   TUMOR:TUMOR的GT信息

4.3 Somatic CNV检测结果

拷贝数变异(Copy number variation, CNV)表现为基因组片段的拷贝数增加或者减少, 是基因组结构变异(Structural variation, SV)的重要组成部分, 可分为缺失(Deletion)和重复(Duplication)两种类型,是一种重要的分子机制。 CNV能够导致包括癌症在内的复杂疾病,对于染色体水平缺失、扩增的研究已经成为肿 瘤研究热点。通常采用Control-FREEC软件 (http://bioinfo-out.curie.fr/projects/freec/) 对Tumor及Normal成对样本检测Somatic CNV,基于比对到参考基因组上Reads的 深度分布情况检测拷贝数的变化,得到如下的统计结果 (具体每个CNV变异的信息详见结果注释文件):

注:该部分结果在06CNV文件夹中。

表4.3 CNV识别结果(部分展示)

注:
   第一列:CNV chromosome
   第二列:start position
   第三列:end position
   第四列:predicted copy number
   第五列:type of alteration
   第六列:genotype

5 Somatic mutation结果注释

利用ANNOVAR(http://www.openbioinformatics.org/annovar/) 软件对检测出的体细胞突变进行注释。
   1)使用Refseq注释变异位点的基因结构,基因类型包括mRNA、非编码RNA等;
   2)变异位点的基因组特征,对于位于基因组重复区段内的突变需谨慎对待;
   3)通过SIFT、PolyPhen以及MutationTaster等方法全面评估非同义突变 对疾病/肿瘤的影响;
   4)我们提供了dbSNP、千人基因组SNP数据库、COSMIC已知肿瘤体细胞突变数据库 和esp6500变异数据库等注释,对变异结果可以进行任何组合的筛选;

注:该部分结果在08variant_annotation文件夹中,该部分结果内容较多, 未在报告中展示。(注释数据库信息见 附表1)

6 附表

附表1 ANNOVAR使用的30个注释数据库相关信息:

refGene FASTA sequences for all annotated transcripts in RefSeq Gene
avsift whole-exome SIFT scores for non-synonymous variants (obselete and should not be uesd any more)
abraom 2.3 million Brazilian genomic variants
cadd13gt20 CADD version 1.3 score>20
clinvar_20210123 Clinvar version 20210123 with separate columns (CLNALLELEID CLNDN CLNDISDB CLNREVSTAT CLNSIG)
cadd13gt20 CADD version 1.3 score>20
dbscsnv11 dbscSNV version 1.1 for splice site prediction by AdaBoost and Random Forest
esp6500siv2_all alternative allele frequency in All subjects in the NHLBI-ESP project with 6500 exomes, including the indel calls and the chrY calls. This is lifted over from hg19 by myself.
esp6500siv2_ea alternative allele frequency in European American subjects in the NHLBI-ESP project with 6500 exomes, including the indel calls and the chrY calls. This is lifted over from hg19 by myself
exac03 ExAC 65000 exome allele frequency data for ALL, AFR (African), AMR (Admixed American), EAS (East Asian), FIN (Finnish), NFE (Non-finnish European), OTH (other), SAS (South Asian)). version 0.3. Left normalization done.
exac03nontcga ExAC on non-TCGA samples (updated header)
gene4denovo201907 gene4denovo database
avsnp150 dbSNP150 with allelic splitting and left-normalization
gerp++elem conserved genomic regions by GERP++
gme Great Middle East allele frequency including NWA (northwest Africa), NEA (northeast Africa), AP (Arabian peninsula), Israel, SD (Syrian desert), TP (Turkish peninsula) and CA (Central Asia)
gwava whole genome GWAVA_region_score and GWAVA_tss_score (GWAVA_unmatched_score has bug in file)
hrcr1 40 million variants from 32K samples in haplotype reference consortium
cadd13gt10 CADD version 1.3 score>10
intervar_20180118 InterVar: clinical interpretation of missense variants (indels not supported)
mcap13 [M-CAP scores v1.3]
mitimpact24 pathogenicity predictions of human mitochondrial missense variants
nci60 NCI-60 human tumor cell line panel exome sequencing allele frequency data
revel REVEL scores for non-synonymous variants
gnomad211_genome gnomAD exome collection (v2.1.1), with "AF AF_popmax AF_male AF_female AF_raw AF_afr AF_sas AF_amr AF_eas AF_nfe AF_fin AF_asj AF_oth non_topmed_AF_popmax non_neuro_AF_popmax non_cancer_AF_popmax controls_AF_popmax" header
gerp++gt2 whole-genome GERP++ scores greater than 2 (RS score threshold of 2 provides high sensitivity while still strongly enriching for truly constrained sites. )
icgc28 International Cancer Genome Consortium version 28
kaviar_20150923 170 million Known VARiants from 13K genomes and 64K exomes in 34 projects
ljb26_all whole-exome SIFT, PolyPhen2 HDIV, PolyPhen2 HVAR, LRT, MutationTaster, MutationAssessor, FATHMM, MetaSVM, MetaLR, VEST, CADD, GERP++, PhyloP and SiPhy scores from dbNSFP version 2.6
eigen whole-genome Eigen scores
popfreq_all_20150413 A database containing all allele frequency from 1000G, ESP6500, ExAC and CG46
regsnpintron prioritize the disease-causing probability of intronic SNVs
gerp++gt2 whole-genome GERP++ scores greater than 2 (RS score threshold of 2 provides high sensitivity while still strongly enriching for truly constrained sites. )

附表2 分析使用软件

分析 软件 参考文献
数据过滤 trimmomatic Trimmomatic: a flexible trimmer for Illumina sequence data.
数据质控 fastqc https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
数据比对 BWA Fast and accurate short read alignment with Burrows–Wheeler transform
SNP识别 GATK4 muTect2 A framework for variation discovery and genotyping using next-generation DNA sequencing data
Indel识别 Strelka Strelka2: fast and accurate calling of germline and somatic variants
CNV识别 Control-FREEC Control-FREEC: a tool for assessing copy number and allelic content using next generation sequencing data
变异注释 ANNOVAR ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data

7. 联系我们

地址:广东省深圳市南山区留仙大道1201号大学城创客小镇

邮编:518055