全外显子组测序分析中预处理方法和变异识别方法的比较

基金项目：国家自然科学基金资助项目（编号：0318，30930082)，国家科技重大专项（编号：2008ZX10002-006,2012ZX10002007），重庆市科技项目(编号：cstc2012gg-yyjsB10007作者简介：闫瑾，男，硕士，研究方向：病毒性肝炎相关性肝细胞肝癌的全外显子组分析研究，电话：135940646474，E-mail：[e-mail]通讯作者：任红，电话：[phone]，E-mail：[e-mail]全外显子组测序分析中预处理方法和变异识别方法的比较闫瑾,潘琦,任红（重庆医科大学附属病毒性肝炎所，重庆400010）【摘要】目的研究在外显子组数据分析中，使用不同的预处理方法和变异过滤方法对变异识别的影响。方法采用FASTX-Toolkit、Trimmomatic作为预处理方法，修饰后不同的不匹配读长（single-endreads）取舍策略，以及硬过滤（hardfilter）和变异质量得分重新校正（VQSR）作为变异过滤方法对两例全外显子组数据进行变异识别，通过数据覆盖深度（DP）、识别变异的数目、Ti/Tv值和基因型一致性等数据进行比较其效果。结果Trimmomatic预处理后的读长的测序覆盖深度与未预处理的原始数据接近，但明显高于FASTX-Toolkit预处理方法。当DP≥10×、基因型质量分数（GQ）≥20时，经Trimmomatic预处理后识别到的SNV数量比FASTX-Toolkit高，与未预处理组接近。当包含single-end读长时，FASTX-Toolkit组多识别的SNV数量高于(28%)Trimmomatic组(5%)。当样本量较少时，在所有试验组中硬过滤方法滤掉的SNV要少于VQSR。结论Trimmomatic修饰（过滤）原始序列更温和，而FASTX-Toolkit可能过度过滤了原始数据。保留single-end读长有利于下游变异识别。硬过滤相较于变异质量得分重新校准表现出更高的容忍度。[关键词]全外显子组测序；预处理；变异识别[中图法分类号]R857.3；Q344+.12Comparisonofmethodsforpre-processingandvariantsfilteringinanalyzingwholeexomesequencingdataYanJin,PanQi,RenHong（InstituteforViralHepatitis，ChongqingMedicalUniversity，Chongqing，400010，China）【Abstract】ObjectiveToinvestigatehowdifferentmethodsforpre-processingandvariantsfilteringaffectvariantscalling.MethodThroughthecalculationofdepthofcoverage,numberofvariants,Ti/Tvratioandnon-referenceconcordance,wecomparetheeffectofFASTX-ToolkitandTrimmomaticinpreprocessingtheexomedata,thestrategiesofsingle-endinclusionand‘Hard’filterandvariantsqualityscorerecalibration(VQSR)invariantsfilterbyusingwholeexomesequencingdatafromtwotestsamples.ResultTrimmomaticpre-processedreadsshowedsimilardepthofcoveragetoreadsthosewithoutpre-processing,butsignificantlygreaterthanthosebyFASTX-Toolkitpre-processedreads.Withdepthofcoverage≥10×andgenotypequality≥20,thenumberofcalledSNVsidentifiedbyTrimmomaticwasgreaterthanFASTX-Toolkit,butsimilartothosewithoutpre-processing.Withtheinclusionofsingle-endreads,thenumberofvariantsincreasedsignificantlyforFASTX-Toolkitpre-processing(~28%)thanTrimmomaticpre-processing(~5%).Intheallsettings,‘Hard’filteringfilteredlessSNVsthanVQSRfilteringinsmallsamplesize.ConclusionSequencereadsweretrimmedand/orfilteredmoderatelybyTrimmomatic,whereasitseemedtobeover-filteredbyFASTX-Toolkit.Keepingthesingle-endreadsisgoodforvariantscallinginthedownstreamanalysis.The‘Hard’filteringshowedamorefavorabletolerabilityprofilethan‘VQSR’filtering.Keyword:wholeexomesequencing,pre-processing,variantsfilteringSupportedbytheGeneralProgramofNationalNaturalScienceFoundationofChina(0318，30930082),NationalScienceandTechnologyMajorProject（2008ZX10002-006,2012ZX10002007）,theFoundationforSci&TechResearchProjectofChongqing(cstc2012gg-yyjsB10007)Correspondingauthor:RenHong，Tel:[phone]，E-mail：[e-mail]自全外显子组技术出现以来，研究者们利用该技术不断揭示了众多孟德尔疾病发病的原因[1]。随着近年来测序技术飞速发展，第二代测序技术应用日趋成熟，费用成本逐渐下降，全外显子组测序被越来越多的实验室和临床检测所应用。虽然第二代测序技术通量大幅提升，测序深度不断提高，在带来更高的碱基识别率的同时，对生物信息学又提出巨大的挑战。由于技术的高速发展以及学界的争议，仍然没有一套公认的、标准的第二代测序数据质量控制方法。是否需要对测序得到的原始序列进行质量控制、以及哪些因素对变异识别造成影响都尚无定论。本研究试图比较不同的预处理方法、预处理后产生的不匹配读长（single-endreads）的取舍策略和变异过滤方法对全外显子组测序数据分析中的测序数据覆盖深度（DepthofCoverage，DP）、识别变异的数目、Ti/Tv值和基因型一致性的影响。1.材料和方法1.1样本我们利用两组1000基因组计划中测序得到的全外显子组数据（NA12878和NA18967，）作为样本进行比较。此两个样本均使用NimbleGen®SeqCapEZExomeprobes(v1.0)进行外显子捕获，并在Illumina®GenomeAnalyzerIIx测序仪上采用76-bp双末端测序方法（Paired-EndSequencing）进行测序[2]，共分别生成13.4GB和17.0GBFastq格式的数据。1.2全外显子组测序分析流程参考GATK建议的第二代测序变异识别和基因分型步骤框架[3,4]，我们设计了本实验全外显子组分析流程图（图1），共分三个阶段。图1本实验室设计及全外显子组数据分析流程图蓝色矩形：输入或输出文件；黄色椭圆：使用的软件名称；红色菱形：条件判断；白色圆框：分析步骤；实线箭头：各阶段内部程序运行方向；虚线箭头：各阶段间程序运行方向。第一阶段为外显子组原始数据处理；第二阶段为变异识别；第三阶段为变异结果统计和质量评价。第一阶段为外显子组原始数据处理。首先，对两例外显子组数据分别用FASTX-Toolkit（v.0.0.13）[5]和Trimmomatic(v0.20)[6]进行预处理，并设无预处理对照组(w/opre-processing)。预处理步骤包括测序接头/引物的剪切（所用引物序列见表1）、依据碱基质量修饰和过滤低质量原始序列（定义质量值小于20为低质量）和滤过人工产物。由于预处理会将原长度一致的配对正反向序列变得读长长度不一，也就是产生single-end读长，由于比对程序无法识别长短不一的读长，故需要把它们与匹配读长（paird-endreads）分开处理。然后，将预处理后生成的数据输入到BWA[7]中，与参考序列（GRCH7）进行比对。因为BWA无法同时处理Paried-end读长和single-end读长，故需要用“sampe”和“samse”两种算法分别运算，最终生成二进制序列比对/图谱文件（BinarySAMfile，BAM）。鉴于经非对称修饰后产生的single-end读长对下游的数据分析能造成多大的影响还不清楚，我们将经过预处理过后的数据分为paired-&single-end组（pse）（使用Picard[8]“MergeSamFiles”整合SAM文件）和paired-end组（pe）分别分析。BWA在进行比对时难以避免会产生错误，尤其在插入和缺失变异（Indel）的片段周围。若不进行校正，这些错误很容易被误认为是SNV。使用GATK[4]（v1.6）RealignerTargetCreator对这些区段进行本地化重比对，Picard工具（MarkDuplicates）移除重复序列，使用CountCovariates和TableRecalibration基于测序循环次数（machinecycle）和周围序列情况（sequencecontext）对碱基质量分数重校正。最后使用BAM文件使用SAMtools[8]进行索引编辑。第二阶段为变异识别。以人类基因组（GRCh37）和dbSNP数据库（release132，build37）作为参考序列，使用GATK“UnifiedGneotyper”工具处理第一阶段生成的结果，识别两个样本单核苷酸变异（SNV）和插入/缺失变异。“UnifiedGenotyper”的参数“stand_emit_conf”和“stand_call_conf”均使用30.0，以忽略所有低质量（Phred-based＜30）的变异，若有基因型未能识别均标记为“N”。使用GATK“selectVariants”对生成的所有变异结果进行挑选，把SNV(singlenucleotidevariant)结果挑选出来。选择硬过表1在FASTX-Toolkit和Trimmomatic中过滤掉的引物（接头）序列。引物（接头）序列Paried_End_Sequencing_Primer_SP1ACACTCTTTCCCTACACGACGCTCTTCCGATCTPaired_End_Sequencing_Primer_SP2CGGTCTCGGCATTCCTACTGAACCGCTCTTCCGATCTPaired_End_PCR_Primer_fAATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCC