作者，Evil Genius~~
外显子生信分析流程
基本分析部分主要用于获取样本基因组的突变信息。 首先，对测序原始数据（Raw data）进行质控，得到高质量的Clean data；然后，将Clean data与人参考基因组序列进行比对分析，获得Bam文件； 最后基于Bam进行Somatic Mutation的检出与注释，从而得到疾病机制研究的基本突变信息结果，基本分析版报告主要展示Somatic相关结果。以下为基本分析的流程图：
测序数据过滤
测序获得的原始数据中包含少量带有测序接头或测序质量较低的reads， 为保证数据分析的质量及可靠性，需要对原始数据进行过滤。 通常使用Trimmomatic软件对测序数据进行过滤。
测序错误率分布检查
每个碱基的测序Phred值（Phred score，Qphred）是由测序错误率通过公式一 转化得到的，而测序错误率是在碱基识别（Base Calling）过程中通过一种判别 发生错误概率的模型计算得到的，其数值对应关系如下表所示。
测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。 对于Illumina高通量测序平台，测序错误率分布具有两个特点：
*（1）测序错误率会随着测序的进行而升高，这是由于测序过程中 荧光标记的不完全切割等因素引起荧光信号衰减，因而导致错误率升高；
*（2）每个Read前几个碱基的位置也会有较高的测序错误率， 这是由于边合成边测序过程初始阶段，测序仪荧光感光元件对焦速度较慢， 获取的荧光图像质量较低，导致碱基识别错误率较高。
在该部分分析中，若样品80%的测序序列错误率在0.1%以下即为合格。
GC含量分布
核苷酸序列中鸟嘌呤(G)和胞嘧啶(C)所占的比例称为GC含量。GC含量在物种间存在一定特异性，但由于反转录过程中所使用的6bp随机引物，会引起前几位碱基在核苷酸组成上有一定偏好性，产生正常波动，随后则趋于稳定。对于NEB普通建库方法，由于序列的随机性打断和双链互补等原则，理论上测序读段 在每个位置的GC及AT含量应分别相等，且在整个测序过程基本稳定不变，呈水平线。而对于链特异性建库而言，由于只保留了单链信息，可能会出现AT分离或GC分离现象。
在该部分结果的图中，正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异，因此好的样本四条线应该是平行且接近的。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示有overrepresented sequence的污染。当所有的碱基比例一致地出现bias时，即四条线平行但分开，往往代表文库有bias（建库过程或本身特点），或是测序中的系统误差。当任一位置A/T比例与G/C比例相差超过10%时，为warn； 当任一位置的A/T与G/G比例相差超过20%时，为fail。
测序数据质量分布
样本中超过80%测序序列的测序质量在Q30以上即为合格，就能保证后续分析的正常进行， 根据测序技术的特点，测序片段末端的碱基质量一般会比前端的低。
测序深度、覆盖度统计
有效测序数据通过BWA比对到参考基因组，得到BAM格式的最初的比对结果。 BAM文件再进行标记重复等处理，从而得到BAM格式的最终比对结果， 再利用比对到参考基因组上的有效数据进行覆盖度的统计。 通常，人类样本的测序Reads能达到98%以上的比对率， 测序深度在10×以上Reads覆盖的位点检测出的SNV比较可信。
胚系突变变异检测
胚系突变（Germline mutations）主要是由于生殖细胞（germ cells）突变导致，生殖细胞在男性中为精源细胞，突变发生在睾丸中；生殖细胞在女性中为卵细胞，突变发生在卵巢中。在大多数情况下，Germline mutations 是沉默的，不对亲本产生影响，除非它们影响配子（gametes）的产生。 尽管这些突变会造成负面影响，例如罕见疾病，甚至癌症，但仍可促进人类健康的遗传多样性
Germline Indel检测结果
InDel全称Insertion and Deletion，指小片段的插入和缺失。编码区或剪接位点处发生的插入缺失都可能会改变蛋白的翻译。移码变异，其插入或缺失的碱基串的长度为3的非整数倍，因此可能导致整个阅读框的改变；与非移码变异比较，移码突变对基因功能的影响更大，同时也受到更大的筛选压力。这里利用MutScan模块检测Germline InDel信息， 对得到的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信息
Germline SNV检测结果
SNV全称Single nucleotide variant，即单核苷酸突变， 指在基因组上由单个核苷酸的替换所引起的变异。主要使用GATK muTect2软件 寻找Somatic 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：FORMAT中内容在该样本中对应的数值
体细胞变异检测
体细胞突变（Somatic mutation）是指除胚系细胞以外的体细胞所发生的突变，是发生在正常机体细胞中的突变，如发生在皮肤和器官。体细胞突变既不遗传自亲本，也不会传递给后代，却可以引起当代某些细胞的遗传结构发生改变。体细胞突变，特别是其中的驱动突变（Driver mutation）对解释肿瘤的发生和发展具有非常重要的意义，另一方面肿瘤耐药性的产生也与体细胞突变有关。 因此，关注体细胞突变是肿瘤基因组研究的重心，也是肿瘤基因组研究区别于疾病研究的一个特性。在信息分析时，通过将肿瘤组织与正常组织相比， 过滤掉正常组织中的胚系细胞突变（Germline mutation），只保留下肿瘤细胞携带的体细胞突变。
Somatic SNV检测结果
SNV全称Single nucleotide variant，即单核苷酸突变， 指在基因组上由单个核苷酸的替换所引起的变异。主要使用GATK muTect2软件 寻找Somatic SNV。
设置过滤条件：reads数≥30条，突变丰度＞1%，包括 “Pathogenic” 的变异位点和“Likely_pathogenic” 的变异位点。
Germline/Somatic mutation结果注释
利用ANNOVAR（http://www.openbioinformatics.org/annovar/） 软件对检测出的体细胞突变进行注释。
* 1）使用Refseq注释变异位点的基因结构，基因类型包括mRNA、非编码RNA等；
* 2）变异位点的基因组特征，对于位于基因组重复区段内的突变需谨慎对待；
* 3）通过SIFT、PolyPhen以及MutationTaster等方法全面评估非同义突变 对疾病/肿瘤的影响；
* 4）提供dbSNP、千人基因组SNP数据库、COSMIC已知肿瘤体细胞突变数据库 和esp6500变异数据库等注释，对变异结果可以进行任何组合的筛选；
注：数据库信息见文末
* 该部分注释结果包含200+列数据，此处选择较为重要的列进行说明：
* 第1列：(Chr)染色体信息
* 第2列：(Start)变异起始位点
* 第3列：(End)变异终止位置
* 第4列：(Ref)参考等位
* 第5列：(Alt)突变等位
* 第6列：(Func.refGene)变异位点在基因组的区域
* 第7列：(Gene.refGene)变异位点注释到的基因
* 第8列：(GeneDetail.refGene) 描述UTR、splicing、ncRNA_splicing或intergenic区域的变异情况。当Func列的值为exonic、ncRNA_exonic、intronic、ncRNA_intronic、upstream、downstream、upstream；downstream、ncRNA_UTR3、ncRNA_UTR5时，该列为空；当Func列的值为exonic;splicing时，表示该位点位于某些转录本的exonic区，另一些转录本的splicing区，这种情况下，GeneDetail会给出该位点对于转录本splicing的影响，例如，NM_1524XX:exon3:c.232C>T，表示该变异位于转录本NM_1524XX上，exon3表示第3个外显子，c.232C>T表示cDNA的232bp处发生由C到T的突变；当Func列的值为intergenic时，该列格式为dist=1322;dist=12414，表示该变异位点距离两侧基因的距离。
* 第9列：(ExonicFunc.refGene) 外显子区的SNV or InDel变异类型（SNV的变异类型包括synonymous_SNV, missense_SNV, stopgain,stopgloss和unknown；InDel的变异类型包括frameshift insertion, frameshift deletion, stopgain, stoploss, nonframeshift insertion, nonframeshift deletion和unknown）。
* 第10列：(AAChange.refGene) 氨基酸改变，只有当Func列为exonic或exonic;splicing时，该列才有结果。
第11列：(Interpro_domain) 蛋白序列和蛋白分类数据库interpro中关于结构域的注释
- 第12列：(avsnp150)该变异在 dbSNP中的 ID
- 第13列：(snp138)该变异在dbsnp138中的 ID
- 第14列：(1000g2015aug_all) 千人基因组计划数据（2015年8月公布的版本）的所有人群中，该变异位点上突变碱基的等位基因频率
- 第15列：(ExAC_ALL) 指在所有人群中，该变异位点上突变碱基的等位基因频率
- 第16列：(ExAC_AFR)该变异在ExAC数据库中非洲人群的等位基因频率
- 第17列：(ExAC_AMR)该变异在ExAC数据库中美国人群的等位基因频率
- 第18列：(ExAC_EAS)该变异在ExAC数据库中东亚人群的等位基因频率
- 第19列：(ExAC_FIN) 该变异在ExAC数据库中芬兰人群的等位基因频率
- 第20列：(ExAC_NFE) 该变异在ExAC数据库中非芬兰欧洲人群的等位基因频率
- 第21列：(ExAC_OTH) 该变异在ExAC数据库中除已指定人群之外的人群等位基因频率
第22列：(ExAC_SAS) 该变异在ExAC数据库中南亚人群的等位基因频率
* 第23列：(CLNALLELEID) the ClinVar Allele ID
* 第24列：(CLNDN) ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB
* 第25列：(CLNDISDB) Tag-value pairs of disease database name and identifier, e.g. OMIM:NNNNNN
* 第26列：(CLNREVSTAT) ClinVar review status for the Variation ID
* 第27列：(CLNSIG) Clinical significance for this single variant
* 第28列：(InterVar_automated) InterVar数据库对遗传变异进行临床解释 - 第29-56列：(PVS1-BP7) 美国医学遗传学与基因组学学会（ACMG）和分子病理协会（AMP）在2015年对临床实验室的基因检测进行了指导和规范（PMID: 25741868）。该指导规范主要就是适用于孟德尔遗传病相关基因变异或者是生殖系变异。指导规范推荐记载突变遵循统一的规范-人类基因组变异协会（HGVS），并变异根据人群基因频率（population data）、软件预测（computational data）和功能试验（functional data）等参数分为五个级别：致病性突变（pathogenic）、可能致病性突变（likelypathogenic）、意义不明突变（uncertain significance）、可能良性突变（likely benign）和良性多态性突变（benign）。该规范列出了致病性/可能致病的各种情况的支持证据，证据强度依次包括超强证据（PVS1）、强证据（PS1-4，注意这里的数字不代表证据强度的区别，仅表示同一证据强度的不同的证据情况，下同）、中度证据（PM1-6）、支持性证据（PP1-5），良性多态性/可能良性证据强度依次包括独立证据（BA1）、强证据（BS1-4）、支持性证据（BP1-6）。
- 第57列：(cosmic95) 通过Cosmic 数据库的注释，可知该突变位点在此前文献中是否出现过，出现在什么癌种，出现了几次。
* 第58列：(SIFT_score) SIFT分值，表示该变异对蛋白序列的影响，SIFT 分值越小越“有害”，表明该SNP导致蛋白结构或功能改变的可能性大；
* 第59列：(SIFT_pred) D: Deleterious (sift<=0.05); T: tolerated (sift>0.05)
* 第60列：(Polyphen2_HDIV_score) 利用 PolyPhen2 基于 HumanDiv 数据库预测该变异对蛋白序列的影响，用于复杂疾病，数值越大越“有害”，表明该 SNP 导致蛋白结构或功能改变的可能性大；damaging (0.453<=pp2_hdiv<=0.956); B: benign (pp2_hdiv<=0.452)
* 第61列：(Polyphen2_HDIV_pred) D 或 P 或 B（D: Probably damaging (>=0.957), P: possibly
* 第62列：(Polyphen2_HVAR_score) 利用 PolyPhen2 基于 HumanVar 数据库预测该变异对蛋白序列的影响，用于单基因遗传病。数值越大越“有害”，表明该 SNP 导致蛋白结构或功能改变的可能性大；
* 第63列：(Polyphen2_HVAR_pred) D 或 P 或 B（D: Probably damaging (>=0.909), P: possibly damaging (0.447<=pp2_hvar<=0.909); B: benign (pp2_hvar<=0.446)
* 第64列：(LRT_score) LRT 分值，表示该变异对蛋白序列的影响，值越大越“有害”，表明该 SNP 导致蛋白结构或功能改变的可能性大。
* 第65列：(LRT_pred) D、N 或者 U(D: Deleterious; N: Neutral; U: Unknown)。
* 第66列：(MutationTaster_score) MutationTaster 分值，表示该变异对蛋白序列的影响，值越大越“有害”，表明该 SNP 导致蛋白结构或功能改变的可能性大。("polymorphism_automatic")
* 第67列：(MutationTaster_pred) A ("disease_causing_automatic"); "D" ("disease_causing");"N" ("polymorphism"); "P" (Polymorphism_automatic)。
* 第68列：(MutationAssessor_score) MutationAssessor预测的致病得分
* 第69列：(MutationAssessor_pred) MutationAssessor根据阈值判断得到的预测分类：H为较高可信度的致病位点，M为中等可信的致病位点，L为低可信度的致病位点，N为无害位点
* 第70列：(FATHMM_score) FATHMM软件预测的致病性得分
* 第71列：(FATHMM_pred) FATHMM根据阈值得到的分类：D为较高可信度的致病位点，P为可信度一般的致病位点
* 第72列：(PROVEAN_score) Protein Variation Effect Analyzer，higher values are more deleterious
* 第73列：(PROVEAN_pred) Protein Variation Effect Analyzer，D: Deleterious; N: Neutral - 第74列：(VEST3_score) Variant effect scoring tool；Random forest classifier， higher values are more deleterious
* 第75列：(CADD_raw) CADD raw score
* 第76列：(CADD_phred) CADD phred-like score，higher values are more deleterious - 第77列：(DANN_score) Deleterious Annotation of genetic variants using Neural Networks，higher values are more deleterious
* 第78列：(fathmm-MKL_coding_score) edicting the effects of both coding and non-coding variants using nucleotide-based HMMs, Score >= 0.5: D; Score < 0.5: T - 第79列：(fathmm-MKL_coding_pred) predicting the effects of both coding and non-coding variants using nucleotide-based HMMs, D: Deleterious; T: Tolerated.
* 第80列：(MetaSVM_score) Using Support vector machine. MetaSVM prediction, higher scores are more deleterious. - 第81列：(MetaSVM_pred) Using Support vector machine. MetaSVM prediction, D: Deleterious; T: Tolerated.
- 第82列：(MetaSVM_score) Using Logistic regression. MetaLR score. higher scores are more deleterious.
- 第83列：(MetaLR_pred) Using Logistic regression. MetaLR prediction. D: Deleterious; T: Tolerated.
- 第84列：(integrated_fitCons_score) Integrate functional assays like ChIP-Seq with conservation measure of transcription factor binding sites. Fitness consequences of functional annotation. higher scores are more deleterious.
- 第85列：(integrated_confidence_value) Fitness consequences of functional annotation. confidence level. higher scores are more deleterious.
* 第86列：(GERP++_RS) GREP++ "rejected substitutions" (RS) score，higher scores are more deleterious
* 第87列：(phyloP7way_vertebrate) Phylogenetic p-values for 7 vertebrate species. higher scores are more deleterious
* 第88列：(phyloP20way_mammalian) Phylogenetic p-values for 20 mammalian species. higher scores are more deleterious
* 第89列：(phastCons7way_vertebrate) PhastCons score for 7 vertebrate species. higher scores are more deleterious. - 第90列：(phastCons20way_mammalian) Phylogenetic p-values for 20 mammalian species. higher scores are more deleterious.
- 第91列：(SiPhy_29way_logOdds) SiPhy log odds score for 29 species. higher scores are more deleterious
* 第92-104列：(Otherinfo1-13) 其他注释信息。
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. ) |
网友评论