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

全外显子组测序分析中预处理方法和变异识别方法的比较闫瑾, 潘琦, 任红（重庆医科大学附属病毒性肝炎所，重庆 400010）【摘要】目的研究在外显子组数据分析中，使用不同的预处理方法和变异过滤方法对变异识别的影响。方法采用FASTX-Toolkit、Trimmomatic作为预处理方法，修饰后不同的不匹配读长（single-end reads）取舍策略，以及硬过滤（hard filter）和变异质量得分重新校正（VQSR）作为变异过滤方法对两例全外显子组数据进行变异识别，通过数据覆盖深度（DP）、识别变异的数目、Ti/Tv值和基因型一致性等数据进行比较其效果。结果Trimmomatic预处理后的读长的测序覆盖深度与未预处理的原始数据接近，但明显高于FASTX-Toolkit预处理方法。当DP10、基因型质量分数（GQ）20时，经Trimmomatic预处理后识别到的SNV数量比FASTX-Toolkit高，与未预处理组接近。当包含single-end读长时，FASTX-Toolkit组多识别的SNV数量高于 (28%)Trimmomatic组 (5%)。当样本量较少时，在所有试验组中硬过滤方法滤掉的SNV要少于VQSR。结论Trimmomatic修饰（过滤）原始序列更温和，而FASTX-Toolkit可能过度过滤了原始数据。保留single-end读长有利于下游变异识别。硬过滤相较于变异质量得分重新校准表现出更高的容忍度。关键词 全外显子组测序；预处理；变异识别中图法分类号 R857.3；Q344+.12Comparison of methods for pre-processing and variants filtering in analyzing whole exome sequencing dataYan Jin, Pan Qi, Ren Hong（Institute for Viral Hepatitis，Chongqing Medical University， Chongqing，400010，China）【Abstract】 ObjectiveTo investigate how different methods for pre-processing and variants filtering affect variants calling.Method Through the calculation of depth of coverage, number of variants, Ti/Tv ratio and non-reference concordance, we compare the effect of FASTX-Toolkit and Trimmomatic in preprocessing the exome data, the strategies of single-end inclusion and Hard filter and variants quality score recalibration (VQSR) in variants filter by using whole exome sequencing data from two test samples.ResultTrimmomatic pre-processed reads showed similar depth of coverage to reads those without pre-processing, but significantly greater than those by FASTX-Toolkit pre-processed reads. With depth of coverage 10 and genotype quality 20, the number of called SNVs identified by Trimmomatic was greater than FASTX-Toolkit, but similar to those without pre-processing. With the inclusion of single-end reads, the number of variants increased significantly for FASTX-Toolkit pre-processing (28%) than Trimmomatic pre-processing (5%). In the all settings, Hard filtering filtered less SNVs than VQSR filtering in small sample size. ConclusionSequence reads were trimmed and/or filtered moderately by Trimmomatic, whereas it seemed to be over-filtered by FASTX-Toolkit. Keeping the single-end reads is good for variants calling in the downstream analysis. The Hard filtering showed a more favorable tolerability profile than VQSR filtering.Keyword: whole exome sequencing, pre-processing, variants filteringSupported by the General Program of National Natural Science Foundation of China (0318，30930082), National Science and Technology Major Project（2008ZX10002-006,2012ZX10002007）, the Foundation for Sci & Tech Research Project of Chongqing (cstc2012gg-yyjsB10007)Corresponding author: Ren Hong，Tel: E-mail：自全外显子组技术出现以来，研究者们利用该技术不断揭示了众多孟德尔疾病发病的原因1。随着近年来测序技术飞速发展，第二代测序技术应用日趋成熟，费用成本逐渐下降，全外显子组测序被越来越多的实验室和临床检测所应用。虽然第二代测序技术通量大幅提升，测序深度不断提高，在带来更高的碱基识别率的同时，对生物信息学又提出巨大的挑战。由于技术的高速发展以及学界的争议，仍然没有一套公认的、标准的第二代测序数据质量控制方法。是否需要对测序得到的原始序列进行质量控制、以及哪些因素对变异识别造成影响都尚无定论。本研究试图比较不同的预处理方法、预处理后产生的不匹配读长（single-end reads）的取舍策略和变异过滤方法对全外显子组测序数据分析中的测序数据覆盖深度（Depth of Coverage，DP）、识别变异的数目、Ti/Tv值和基因型一致性的影响。1.材料和方法1.1样本我们利用两组1000基因组计划中测序得到的全外显子组数据（NA12878 和NA18967，/swansonw/Swanson_Lab/Data.html）作为样本进行比较。此两个样本均使用NimbleGen SeqCap EZ Exome probes (v1.0)进行外显子捕获，并在Illumina Genome Analyzer IIx测序仪上采用76-bp双末端测序方法（Paired-End Sequencing）进行测序2，共分别生成13.4GB和17.0GB Fastq格式的数据。1.2全外显子组测序分析流程参考GATK建议的第二代测序变异识别和基因分型步骤框架3, 4，我们设计了本实验全外显子组分析流程图（图1），共分三个阶段。图1 本实验室设计及全外显子组数据分析流程图蓝色矩形：输入或输出文件；黄色椭圆：使用的软件名称；红色菱形：条件判断；白色圆框：分析步骤；实线箭头：各阶段内部程序运行方向；虚线箭头：各阶段间程序运行方向。第一阶段为外显子组原始数据处理；第二阶段为变异识别；第三阶段为变异结果统计和质量评价。第一阶段为外显子组原始数据处理。首先，对两例外显子组数据分别用FASTX-Toolkit（v.0.0.13）5和Trimmomatic(v0.20)6进行预处理，并设无预处理对照组(w/o pre-processing)。预处理步骤包括测序接头/引物的剪切（所用引物序列见表1）、依据碱基质量修饰和过滤低质量原始序列（定义质量值小于20为低质量）和滤过人工产物。由于预处理会将原长度一致的配对正反向序列变得读长长度不一，也就是产生single-end读长，由于比对程序无法识别长短不一的读长，故需要把它们与匹配读长（paird-end reads）分开处理。然后，将预处理后生成的数据输入到BWA7中，与参考序列（GRCH7）进行比对。因为BWA无法同时处理Paried-end读长和single-end读长，故需要用“sampe”和“samse”两种算法分别运算，最终生成二进制序列比对/图谱文件（Binary SAM file，BAM）。鉴于经非对称修饰后产生的single-end读长对下游的数据分析能造成多大的影响还不清楚，我们将经过预处理过后的数据分为paired-&single-end组（pse）（使用Picard8“MergeSamFiles”整合SAM文件）和paired-end组（pe）分别分析。 表1 在 FASTX-Toolkit和Trimmomatic中过滤掉的引物（接头）序列。引物（接头）序列Paried_End_Sequencing_Primer_SP1ACACTCTTTCCCTACACGACGCTCTTCCGATCTPaired_End_Sequencing_Primer_SP2CGGTCTCGGCATTCCTACTGAACCGCTCTTCCGATCTPaired_End_PCR_Primer_fAATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTPaired_End_PCR_Primer_rCAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCTPaired_End_Adapter_A2GATCGGAAGAGCGGTTCAGCAGGAATGCCGAGBWA在进行比对时难以避免会产生错误，尤其在插入和缺失变异（Indel）的片段周围。若不进行校正，这些错误很容易被误认为是SNV。使用GATK4（v1.6）RealignerTargetCreator对这些区段进行本地化重比对，Picard工具（MarkDuplicates）移除重复序列，使用CountCovariates和TableRecalibration基于测序循环次数（machine cycle）和周围序列情况（sequence context）对碱基质量分数重校正。最后使用BAM文件使用SAMtools8进行索引编辑。第二阶段为变异识别。以人类基因组（GRCh37）和dbSNP数据库（release 132，build 37）作为参考序列，使用GATK“UnifiedGneotyper”工具处理第一阶段生成的结果，识别两个样本单核苷酸变异（SNV）和插入/缺失变异。“UnifiedGenotyper”的参数“stand_emit_conf ”和“ stand_call_conf”均使用30.0，以忽略所有低质量（Phred-based30）的变异，若有基因型未能识别均标记为“N”。使用GATK“selectVariants”对生成的所有变异结果进行挑选，把SNV(single nucleotide variant)结果挑选出来。选择硬过滤（Hard filtering）和变异质量分数重校正(variant quality score recalibration，VQSR)两种方法来过滤原始SNVs以比较两者过滤效果。将硬过滤工具“VariantFiltration”六个过滤参数设为固定值：“QD (quality by depth) 2.0”、“MQ (mapping quality) 60.0”、“HaplotypeScore (haplotype score) 13.0”、“MQRankSum (mapping quality rank sum test) -12.5”以及“ReadPosRankSum (read position rank sum test) -8.0”。使用固定参数：“maxGaussinas 4”和“percentBadVariants 0.05”运行VQSR“VariantRecalibrator”和“ApplyRecalibration”。第三阶段为变异结果统计和质量评价。依据SNV和Indel上的标记特征和参数，使用VCFtools9进行过滤。针对硬过滤后得到的结果，使用QDFilter, MQFilter, FSFilter, HaplotypeScoreFilter, MQRankSumFilter, 和 ReadPosRankSumFilter过滤器对SNV进行过滤。把VQSR过滤后结果凡是满足“TruthSensitivityTranche99.90to100.00” 和“TruthSensitivityTranche99.00to99.90”的SNV滤掉。最后，运行“VariantEval”工具中的GenotypeConcordance模块将实验中识别出的变异分别与NA12878和NA18967在HapMap3和Omni的变异结果10进行比较。生成的评价指标有SNVs数量、Ti/Tv比值（转换/颠换率）、dbSNP率（即能在dbSNP数据中找到的变异数目的比例）、非参考变异敏感度（Non-reference Sensitivity，NRS）、非参考变异矛盾率（Non-reference Discrepancy，NRD）。1.3覆盖深度分析覆盖深度是测序技术的重要指标之一，当测序长度不变时，提高测序深度可以获得更好的比对率和变异识别率。Nimblegen SeqCap EZ Exome (v1) 探针中包含了25.2MB目标区域（Targeted-region）和34.1MB平铺区域（tiled-region），囊括了CCDS数据库（V20080430）中所公布的11所有无重复的蛋白质编码区域（包括外显子上下游200bp以内），以及约550个miRNAs序列。应用UCSC“liftover”工具，将Nimblegen SeqCap EZ (v1) (hg18) tiled区域的坐标系转换成hg19，共生成了176159个连续的捕获基因组区域（总计34108798 bp），与之前hg18坐标系的数目相近（176,817个连续的捕获基因组区域总计34,108,810 bp）。将GRCh37/hg19人类基因组序列refGene（从UCSC基因组浏览器中下载，2012.6.3）作为参考基因序列，定位于SeqCap Ez Exome（V1）基因组捕获的区间（从Nimblegen公司获得），使用GATK“DepthOfCoverage”工具对预处理后生成的BAM文件进行覆盖深度分析。2.1预处理后的读长使用FASTX-Toolkit、Trimmomatic两个样本后读长剩余数量见表2。FASTX-Toolkit滤过后保留下的读长数明显少于经Trimmomatic处理后的。在保留下的读长中，经FASTX-Toolkit处理后single-end读长所占比例(在NA12878为30.4%，在NA18967中为35.5%)要显著高于Trimmomatic处理后的(在NA12878为3.3%，在NA18967中为8.0%)。若将删除重复序列后的原始序列定义为唯一的读长，由比对后的结果可知，经FASTX-Toolkit处理后比对上参考基因组的读长要显著低于Trimmomatic处理后的以及未预处理的。表2 预处理后、重复序列及比对后读长数量统计处理方式过滤后读长数量重复序列比对后读长d正向序列反向序列总数cNA12878 未预处理375207503752075075041500627071961879633 (90.0%) FASTX psea15901894+752442015901894+638571545713923 (60.9%)1045890235061064 (51.0%) FASTX peb159018941590189431803788 (42.4%)238942429295863 (42.6%) Trimmomatic pse30224628+117319230224628+91925262541700 (83.3%)597008756178578 (81.7%) Trimmomatic pe302246283022462860449256 (80.6%)468217855407250 (80.6%)NA18967 未预处理476280424762804295256084541317765097743 (72.5%) FASTX pse13150803+734082613150803+711766640760098 (42.8%)876877031841000 (35.4%) FASTX pe131508031315080326301606 (27.6%)124936424968823 (27.8%) Trimmomatic pse25761311+221118125761311+225073555984538 (58.8%)487168150297732 (56.0%) Trimmomatic pe257613112576131151522622 (54.1%)241440248339432 (53.8%)apse：paired-end读长和single-end读长。bpe：paired-end读长。c总数：等于正向序列与反向序列数目之和；括号中百分比表示过滤后剩余读长数目占原始读长的比例。d比对后读长：为经过预处理或未经过预处理并去除重复序列后能与参考基因组比对上的读长；括号中的百分比表示比对上的读长占原始读长的比例。2.2覆盖深度每个样本有5个BAM文件被用来检测覆盖深度（图2和表3），包括未预处理的原始读长、FASTX-Toolkit预处理生成的读长（pse组和pe组），以及Trimmomatic预处理后生成的读长（pse组和pe组）。从图2a中可看出未处理的原始序列的覆盖深度最高，NA12878和NA18967均有84%的碱基20，平均覆盖率分别为99和101.88。经Trimmomatic预处理后的读长与原始读长覆盖深度接近，也远高于FASTX-Toolkit。single-end读长的加入对FASTX-Toolkit数据的覆盖深度提升有较显著的影响（14%24%）。以上结果说明使用FASTX-Toolkit做预处理会极大的降低原始数据的碱基位点覆盖深度和平均覆盖深度。 表3 使用不同的预处理方法后tiled区域的碱基覆盖深度分析。 处理方法覆盖深度1 (bp)1 (%)20 (bp)20 (%)平均值 ()NA12878 未预处理3375590999286908198499.00 FASTX psea3323902498192506045637.19 FASTX peb3282708396166210194932.50 Trimmomatic pse3372355099277034788184.05 Trimmomatic pe3371301899275821428183.47NA18967 未预处理33699948992870298484101.88 FASTX pse3306727997174213025131.52 FASTX pe3247846895140695864126.27 Trimmomatic pse3359803699263355807771.10 Trimmomatic pe3356132298259704917669.72apse：paired-end读长和single-end读长bpe：paired-end读长使用refGene（含有23652唯一的基因）作为覆盖率分析的参考序列，总共有15737个（占refGene的66.5%）参考基因被捕获，少于CCDS作为参考序列捕获的蛋白质编码基因(16188个)12。NimbleGen SeqCap EZ (v1.0) 捕获到的基因经过预处理后覆盖深度分布见图2b，当测序DP20时，FASTX-Toolkit处理后的基因数目最多。例在NA18967中，FASTX-Toolkit过滤后有23.3%（pse）和38.1%（pe）的捕获基因DP20；而使用Trimmomatic，仅有4.2% (pse) and 4.7% (pe)的捕获基因DP20。以上结果证明不同的预处理方法会对文库捕获的基因测序覆盖深度造成显著的差别，因此也会影响了识别重要（医学相关）基因的效能（power）。图2 （a）目标区域碱基的测序覆盖深度占所有碱基覆盖深度比例连续曲线；（2）目标基因的平均覆盖深度与基因数量的关系曲线。2.3变异评估对测序数据基因型质量影响最大的两个因素是测序深度和基因型质量分数（Genotype Quality Score，GQ）13。GATK认为外显子组测序数据，比对错误和测序深度之间的关系并不明确，故不建议在外显子组数据处理中直接应用DP过滤器。然而，序列中某一位点的测序深度和该样本识别出的变异数目之间的关系尚不清楚，故分析位点测序深度的阈值与识别变异数目的关系并作图（图3a）。图中可见在不同的方案组中，SNV的数据均随着测序深度的增加而迅速减少，特别是在当测序深度在5的时候。当测序深度10时，识别变异的数量趋于相似。从图3b表现了SNV数目和基因型质量的关系。表4 经不同分析分案生成的测序覆盖深度10、基因型质量20的SNV数量 预处理方法变异过滤方法SNV数Ti/Tv基因型一致性总数已知SNV新SNVdbSNP (%)总数已知SNV新SNVHapMapOmniNRScNRD rated(%)NRSNRD rate(%)NA12878 未预处理Hard 418444102681898.02.672.691.971000.121000.30VQSR390263872030699.22.722.732.121000.111000.28 FASTX_pseaHard 276482724140798.52.792.821.661000.391000.58VQSR273102697733398.82.802.831.471000.381000.56 FASTX_pebHard 239302354338798.42.842.871.601000.401000.64VQSR225502235419699.12.902.921.581000.391000.60 Trimmomatic_pseHard 401463942172598.22.692.701.941000.111000.30VQSR370663680326399.32.742.751.801000.091000.28 Trimmomatic_peHard 398673914672198.22.692.711.981000.111000.30VQSR365083625625299.32.752.751.961000.101000.28NA18967 未预处理Hard 3944537973147296.32.742.752.311000.291000.37VQSR369693602194897.42.792.802.441000.261000.35 FASTX_pseHard 237872299479396.72.862.892.171000.591000.71VQSR234692274772296.92.872.892.171000.581000.72 FASTX_peHard 193861872666096.62.912.942.251000.561000.68VQSR182991781648397.42.952.972.401000.561000.65 Trimmomatic_pseHard 3550634227127996.42.772.782.581000.241000.34VQSR329193206685397.42.822.822.761000.211000.30 Trimmomatic_peHard 3476433513125196.42.772.782.501000.261000.36VQSR320263121980797.52.832.832.771000.231000.33apse: paired-end读长和single-end读长bpe: paired-end读长cNRS ：Non-reference sensitivity，非参考序列敏感度dNRDR：Non-reference discrepancy rate，非参考序列差异度我们以DP1310、GQ20为标准，统计识别到的SNV数目、Ti/Tv比值以及与HapMap和OMNI中已识别的基因型一致性。从表4中可以看到不同的预处理组间变异数目均有明显差别，未经过预处理的和用Trimmomatic预处理能识别出更多的变异。未预处理组中的变异数目约是使用FASTX-Toolkit预处理得到的两倍。当分析中包含single-end读长时，在FASTX-Toolkit组中SNV识别率升高15.5%28.3%，而是在Trimmomatic中仅多识别了0.7%2.8%。而在比较两种基因型过滤方法中，硬过滤则能识别出更多的变异（表5），尤其是在未预处理组和Trimmomatic组。经Trimmomatic预处理后虽然变异数量下降但dbSNP率上升（除NA18967 pse+VQSR组持平），而Fastx-Toolkit变异数量和dbSNP率均下降。dbSNP率在VQSR组中要比硬过滤组中相对要高。表5 硬过滤和VQSR过滤后同组内SNV结果相同的和不同的数目。预处理方法相同不同HardVQSRNA18967 未预处理368582587111 FASTX_psea22855932614 FASTX_peb179861400313 Trimmomatic_pse32869263750 Trimmomatic_pe31991277335NA12878 未预处理38958288668 FASTX_pse266111037699 FASTX_pe221871743363 Trimmomatic_pse37019312747 Trimmomatic_pe36473339435apse: paired-end读长和single-end读长bpe: paired-end读长图3 （a）在不同的测序覆盖深度下变异数量的变化曲线；（b）在不同基因型质量值下变异数量的变化曲线。Ti/Tv值是对于评估SNV质量新变异特异性的关键指标，Ti/Tv越高越准确14。对于在外显子里的SNV，Ti/Tv值约为3.0，其余的约为2.015。在本实验中，在新变异位点中，Trimmomatic预处理后识别的新变异Ti/Tv值均比使用相同single-end取舍策略的FASTX-Toolkit组高。已知变异位点Ti/Tv值则均与之前基于1000基因组数据得到的估计值2.8相近16，但却与使用相同样本不同外显子捕获方法（Agilent exome hybrid capture 28MB）和不同测序方法(Illumina GA （约150，93%碱基20）)的结果要低3。所有变异均能在HapMap和OMNI数据库中找到（NRS=100%），差异性也与之前的研究一致,均小于1%（表4）14。2.4在性染色体上的变异男性性染色体上的发生的变异应是纯合子，而在女性染色体上的变异显然不能出现在Y染色体上。所以我们以此标准检查两个实验对象（NA12878为女性，NA18967为男性）在本流程中发生在性染色体上的错误（图4）。在NA18967中，最多有79个和25个杂合性SNV分别标明在X染色体上和Y染色体上；在NA12878中，最多有40个SNV发生在Y染色体上。图4. 性染色体上SNV分布柱状图。NA18967中的X_hetero和Y_hetero与NA12878中的Y_hetero和Y_homo均为错误识别。Hetero：杂合子；Homo：纯合子。3.1讨论选择NA12878和NA18967这两例样本作为测试对象，不仅是因为他们的数据公开，同时还因为他们已经被作为实验对象反复分析过，有既有的数据可以参考2, 3。FASTX-Toolkit与Trimmomatic都是流行的质控预处理软件，而通过本次实验，我们发现原始数据再经过预处理后，测序覆盖深度下降明显，会直接导致识别变异数目降低，特别是FASTX-Toolkit。Trimmomatic的优势在于能利用配对信息同时处理正反向序列，采用滑动窗口剪切（sliding-window cutting），故比FASTX-Toolkit逐碱基剪切（one-base-by-one-base）方法灵活，同时Trimmomatic能在质量控制的同时降低假阳性率上，使dbSNP率升高，以及有更高的新SNV识别准确率。虽然Trimmomatic在本次试验中体现出了较大的优势，但是还不足以说服对测序原始数据使用预处理软件进行质量控制。因为质量差并不意味着某位点变异就错误的，而过多的引物序列和人工产物的确会降低序列比对的正确率，增加运算时间。因此，仅对原始序列做剪切掉已知污染物序列而不对质量值进行修饰过滤，也不失为一种折中的办法。预处理必定会造成原本配对的双末端序列长短不一，产生single-end读长。 FASTX-Toolkit预处理后产生的single-end reads数目明显高于Trimmomatic，故其在包含single-end reads后测序深度增加了接近5，比对上参考序列的读长升高了19.7%，识别的SNVs也升高了15.5%28.3%SNVs。但在Trimmomatic中此差别并不明显，所以保留single-end reads对FASTX-Toolkit预处理后的数据更加重要。若依据GATK的建议，使用VQSR要达到最佳结果至少需要30个样本。但是在实际操作中会因为诸多限制，如试验对象数量过少、经费有限、计算机硬件资源不足等，而很难达到此标准。若样本量与建议的相差不是太大，除了增加样本量以外，加深测序深度，提高测序覆盖率，效果应该也等同于增加样本数目。而在处理小样本测序数据时，则只能选择使用固定了参数的VQSR（如maximum Gaussians、percentBadVariants）或使用硬过滤。本次研究存在以下局限性：1.使用的样本量较少，提高样本量有助于更进一步验证预处理和变异过滤在全外显子组测序分析中的效果；2.使用于外显子捕获的文库NimbleGen SeqCap EZ exome（v1）发布于2009年，故其捕获效率和捕获基因数量可能会较新版本或者其他外显子组捕获试剂稍差12, 17, 18；3. FASTX-Toolkit方法将原始序列过度过滤的原因有待于进一步研究。最后，虽然出现本次实验中的性染色体错配的数量较小（仅约为所有SNV的0.3%），对结果影响不是太大，但还是有必要完善现有比对程序或开发一种新的软件，以识别并纠正错配到性染色体上的变异。当完成本文时，基因组测序分析小组（Genome Sequencing and Analysis Group, GSA）发布了GATK v4，并更新了变异识别的最佳流程。综上所述，由于不同的算法识别到的变异也差别较大，所以是否需要进行预处理，仍然需要进一步研究。不管是对全基因组分析还是全外显子组分析，建立一套标准化的分析流程仍然是将来的目标。参考文献 1 Gilissen C, Hoischen A, Brunner H G, et al. Unlocking Mendelian disease using exome sequencingJ. Genome Biol, 2011, 12(9): 228. 2 George R D, Mcvicker G, Diederich R, et al. Trans genomic capture and sequencing of primate exomes reveals new targets of positive selectionJ. Genome Res, 2011, 21(10): 1686-1694. 3 Depristo M A, Banks E, Poplin R, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing dataJ. Nat Genet, 2011, 43(5): 491-498. 4 Mckenna A, Hanna M, Banks E, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing dataJ. Genome Res, 2010, 20(9): 1297-1303. 5 Gj H. FASTX-Toolkit: FASTQ/A short-reads pre-processing toolsDB/CD. 6 Lohse M, Bolger A M, Nagel A, et al. RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomicsJ. Nucleic Acids Res, 2012, 40(Web Server issue): W622-W627. 7 Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transformJ. Bioinformatics, 2009, 25(14): 1754-1760. 8 Li H, Handsaker B, Wysoker A, et al. The Sequence Alignment/Map format and