1.本发明涉及生物信息分析
技术领域：
：，具体来说，涉及一种人类基因全外显检测数据的分析方法。
背景技术：
：：2.新一代技术的迅猛发展在数据通量和成本上都显示出巨大的优势。全外显子组捕获测序技术(wes)针对外显子功能区域进行深度测序，可以更全面的检测编码区域的变异，且能够识别低频的和新的变异位点。全外显子测序包含目标区间的捕获、文库构建和上机测序，以及生物信息学分析三个过程。3.随着高通量测序技术的不断推广，涌现出海量的基因组测序数据，如何方便、系统性地挖掘这些大数据中的信息，以更全面的解析疾病，服务于生物医学，为数据分析团队提出了更高的要求和新的挑战。技术实现要素：4.针对相关技术中的上述技术问题，本发明提出一种人类基因全外显检测数据的分析方法，能够克服现有技术的上述不足。5.为实现上述技术目的，本发明的技术方案是这样实现的：6.一种人类基因全外显检测数据的分析方法，该方法包括以下步骤：7.s1原始数据质控；8.s2数据比对，下载参考基因hg19，为参考基因构建索引，用bwa(burrows-wheeler-alignmenttool)将fastq文件的read比对到参考基因，并生成sam文件，用samtools软件将sam文件转换为bam文件；9.s3排序，对bam文件进行排序；10.s4去重，使用gatk软件对排序后bam文件进行重复序列标记；11.s5数据矫正，使用gatk软件对数据进行重新矫正；12.s6变异检测，进行variantcalling，寻找snp和indel，将比对数据存储在vcf格式的文件中；13.s7变异过滤，变异如果没有预测性的或已知的功能改变，或者被认为是太常见而不能成为特定疾病的原因，或者不符合所观察到的疾病的遗传模式，则排除；14.s8变异注释，通过annovar工具进行基因文件注释。15.进一步地，s1-s8通过shell脚本语言编写自动化运行程序。16.进一步地，s1中利用fastqc软件进行原始数据质控。17.进一步地，所述vcf格式的文件储存了测序数据通过对比参考基因后得到的变异结果。18.本发明的有益效果：本发明的人类基因全外显检测数据的分析方法，易于推广，采用实用的自动化分析；外显子占全基因组的1％，约90％的人类致病基因变异集中在外显子区域，相比于全基因组测序，全外显子测序更加经济、高效；能够方便快速、系统性地挖掘大数据中的信息，以更全面的解析疾病，服务于生物医学和科研，为更多的生物数据分析团队，尤其是不擅长计算机的生物医学和科研工作者提供便利。附图说明19.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。20.图1是根据本发明实施例所述的人类基因全外显检测数据的分析方法的示意图；21.图2是根据本发明实施例所述的摘要列出质控的每个条目的示意图；22.图3是根据本发明实施例所述的重复数据的示意图。具体实施方式23.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。24.如图1所示，根据本发明实施例所述的人类基因全外显检测数据的分析方法，该方法包括以下步骤：25.s1原始数据质控，上机测序的过程中影响测序结果的因素有很多，测序数据的质量好坏会影响下游分析，测序数据的第一步就是要进行测序数据的质量检查，利用fastqc软件进行原始数据质控；26.s2数据比对，下载参考基因hg19，为参考基因构建索引，用bwa将fastq文件的read比对到参考基因，并生成sam文件，用samtools软件将sam文件转换为bam文件；27.s3排序，对bam文件进行排序；28.s4去重，使用gatk软件对排序后bam文件进行重复序列标记；29.s5数据矫正，使用gatk软件对数据进行重新矫正；30.s6变异检测，进行variantcalling，寻找snp(单个碱基上的变异)和indel(插入缺失标记)，将比对数据存储在vcf格式的文件中；31.s7变异过滤，变异如果没有预测性的或已知的功能改变，或者被认为是太常见而不能成为特定疾病的原因，或者不符合所观察到的疾病的遗传模式，则排除；32.s8变异注释，所述vcf格式的文件储存了测序数据通过对比参考基因后得到的变异结果，通过vcf过滤，后续可根据需求对vcf进行对比到已知变异的数据库，得到碱基突变是否致病，通过annovar工具进行基因文件注释。33.以上步骤全是分步执行，过程繁琐，考虑到实用性，本发明通过shell脚本语言编写自动化运行程序。34.其中，fastqc是基于java的软件，可以快速地对测序数据进行质量评估。35.bwa是能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。它有三个不同的算法：36.bwa-backtrack：是用来比对illumina的序列的，reads长度最长能到100bp。37.bwa-sw：用于比对long-read，支持的长度为70bp-1mbp；同时支持剪接性比对(splitalignments)。38.bwa-mem：推荐使用的算法，支持较长的read长度，同时支持剪接性比对，但是bwa-mem是更新的算法，也更快更准确，且bwa-mem对于70bp-100bp的illumina数据来说效果更好。39.使用bwa整个比对过程分为两步，第一步使用索引命令构建参考基因组的索引，第二步使用bwamem进行比对。40.bwa的使用需要两种输入文件：41.referencegenomedata(fasta格式.fa,.fasta,.fna)；42.shortreadsdata(fastaq格式.fastaq,.fq)。43.annovar是一个高效的注释工具，能够利用最新的数据来分析各种基因组中的遗传变异。由perl编写，支持包括vcf在内的多种输入和输出文件格式。主要包含三种不同的注释方法：gene-based，region-based和filter-based。基于基因的注释(gene-basedannotation)揭示variant与已知基因直接的关系以及对其产生的功能性影响，基于区域的注释(region-basedannotation)揭示variant与不同基因组特定段的关系，例如：它是否落在转录因子结合区域等，基于筛选的注释(filter-basedannotation)则分析变异位点是否位于指定的数据库中，比如dbsnp，1000g，esp6500等数据库，计算sift/polyphen/lrt/mutationtaster/mutationassessor等。44.实施例145.通过测序仪测序后得到一对fastq原始测序文件，分析使用的软件为fastqc，bwa，samtools，igv，gatk，annovar，vcftools，picard。46.s1原始数据质控，上机测序的过程中影响测序结果的因素有很多，测序数据的质量好坏会影响下游分析，测序数据的第一步就是要进行测序数据的质量检查，利用fastqc软件进行原始数据质控，qc完生成了两个文件.html和.zip，.html文件是qc的分析结果网页文件，.zip压缩包中是每条结果的资源文件，如图2所示，摘要列出质控的每个条目，此报表显示，数据通过标准，可以进行后续分析；如果有未通过的项目，则需要考虑进行样本过滤，甚至重新测序；47.s2数据比对，下载参考基因hg19ucsc.hg19.fasta.gz，解压ucsc.hg19.fasta.gz为hg19.fasta，为参考基因构建索引：48.[root@qjyzc/home/hg19fa05:27:14]#bwaindexucsc.hg19.fasta[0049]#log[0050][bwa_index]packfasta...29.32sec[0051][bwa_index]constructbwtforthepackedsequence...[0052][bwtinccreate]textlength＝6274322528,availableword＝453484340[0053][bwtincconstructfrompacked]10iterationsdone.100000000charactersprocessed.[0054]...............[0055][bwtincconstructfrompacked]690iterationsdone.6264217232charactersprocessed.[0056][bwt_gen]finishedconstructingbwtin695iterations.[0057][bwa_index]2864.85secondselapse.[0058][bwa_index]updatebwt...36.25sec[0059][bwa_index]packforward-onlyfasta...16.66sec[0060][bwa_index]constructsafrombwtandocc...927.12sec[0061][main]version:0.7.15-r1140[0062][main]cmd:bwaindexhg19.fasta[0063][main]realtime:3958.820sec；cpu:3874.214sec[0064]#用时65分钟[0065]real65m58.887s[0066]user63m46.305s[0067]sys0m47.911s[0068]完成流程后生成如下文件：[0069]├──[3.0g]hg19.fasta[0070]├──[8.4k]hg19.fasta.amb[0071]├──[3.9k]hg19.fasta.ann[0072]├──[2.9g]hg19.fasta.bwt[0073]├──[3.5k]hg19.fasta.fai[0074]├──[748m]hg19.fasta.pac[0075]└──[1.5g]hg19.fasta.sa[0076]用bwa将fastq文件的read比对到参考基因，并生成sam文件；[0077]bwa比对命令用法[0078]usage:bwamem[options]《idxbase》《in1.fq》[in2.fq][0079]其中，[options]是一系列可选的参数，这里的《idxbase》要输入的是参考基因组的bw索引文件，我们上面通过bwaindex构建好的以human.fasta为前缀的文件便是；《in1.fq》和[in2.fq]输入的是质控后的fastq文件。[0080][root@qjyzc/home/wes06:48:34]#cdbamdir/[0081][root@qjyzc/home/wes/bamdir06:48:40]#bwamem\[0082]-t4-r'@rg\tid:foo\tpl:illumina\tlb:library\tsm:bar'\[0083]/home/hg19fa/hg19.fasta\[0084]/home/wes_fastq/s_r1.fq.gz\[0085]/home/wes_fastq/s_r2.fq.gz》s.sam[0086]#log[0087]timebwamem-t4-r'@rg\tid:foo\tpl:illumina\tlb:library\tsm:bar'/home/hg19fa/hg19.fasta/home/wes_fastq/s_r1.fq.gz/home/wes_fastq/s_r2.fq.gz》s.sam[0088]#log[0089][m::bwa_idx_load_from_disk]read0altcontigs[0090][m::process]read266668sequences(40000200bp)...[0091][m::process]read266668sequences(40000200bp)...[0092][m::mem_pestat]#candidateuniquepairsfor(ff,fr,rf,rr):(0,121247,0,0)[0093][m::mem_pestat]skiporientationffastherearenotenoughpairs[0094][m::mem_pestat]analyzinginsertsizedistributionfororientationfr...[0095][m::mem_pestat](25,50,75)percentile:(194,220,250)[0096][m::mem_pestat]lowandhighboundariesforcomputingmeanandstd.dev:(82,362)[0097][m::mem_pestat]meanandstd.dev:(223.41,42.78)[0098][m::mem_pestat]lowandhighboundariesforproperpairs:(26,418)[0099][m::mem_pestat]skiporientationrfastherearenotenoughpairs[0100][m::mem_pestat]skiporientationrrastherearenotenoughpairs[0101][m::mem_process_seqs]processed266668readsin45.116cpusec,12.066realsec[0102][m::process]read266668sequences(40000200bp)...[0103]........[0104][main]version:0.7.15-r1140[0105][main]cmd:bwamem-t4-r@rg\tid:foo\tpl:illumina\tlb:library\tsm:bar/home/hg19fa/hg19.fasta/home/wes_fastq/s_r1.fq.gz/home/wes_fastq/s_r2.fq.gz[0106][main]realtime:2846.717sec；cpu:11552.456sec[0107]#用时47分钟[0108]real47m26.766s[0109]user191m14.005s[0110]sys1m18.471s[0111][root@qjyzc/home/wes/bamdir11:35:33]#tt[0112].[0113]└──[28g]s.sam[0114]比对后生成了一个巨大的sam文件，比fastq大了近10倍，因此，为了有效节省磁盘空间，一般都会用samtools将它转化为bam文件(sam的特殊二进制格式)，而且bam会更加方便于后续的分析。[0115]sam格式转换为bam格式：[0116][root@qjyzc/home/wes_fastq10:09:50]#samtoolsview-s-bs.sam》s.bam[0117][root@qjyzc/home/wes/bamdir11:59:18]#tt[0118].[0119]├──[6.2g]s.bam[0120]└──[28g]s.sam[0121]根据以上比对的命令可以把bwa和samtools结合输出bam文件。[0122][root@qjyzc/home/wes_fastq]#bwamem\[0123]-t4-r'@rg\tid:foo\tpl:illumina\tlb:library\tsm:bar'\[0124]/home/hg19fa/hg19.fasta\[0125]/home/wes_fastq/s_r1.fq.gz\[0126]/home/wes_fastq/s_r2.fq.gz\[0127]|samtoolsview-s-b-》s.bam[0128]我们通过管道(“|”)把比对的输出如同引导水流一样导流给samtools去处理，上面samtoolsview的-b参数指的就是输出为bam文件，这里需要注意的地方是-b后面的'-'，它代表就是上面管道引流过来的数据，经过samtools转换之后我们再重定向为bam文件。[0129]s3排序，比对结束后，使用samtools对s.bam文件进行排序；[0130][root@qjyzc/home/wes_fastq]#samtoolssort-@4-m3g-obam-os.sorted.bams.bam[0131][root@qjyzc/home/wes/bamdir12:32:49]#tt[0132].[0133]├──[6.2g]s.bam[0134]├──[28g]s.sam[0135]└──[3.4g]s.sorted.bam[0136]排序后的bam按染色体的序号进行排序；[0137]s4去重，如图3所示，进行查看发现有部分数据重复，这些由于pcr扩增引起的重复数据将会影响下游分析，使用gtak软件对重复reads进行标记；[0138]s5数据矫正(碱基质量重校)，使用gatk软件对数据进行重新矫正；[0139]s6变异检测，进行variantcalling，寻找snp(单个碱基上的变异)和indel(插入缺失标记)，将比对数据存储在vcf格式的文件中；[0140]s7变异过滤，变异如果没有预测性的或已知的功能改变，或者被认为是太常见而不能成为特定疾病的原因，或者不符合所观察到的疾病的遗传模式，则排除；[0141]s8变异注释，所述vcf格式的文件储存了测序数据通过对比参考基因后得到的变异结果，通过vcf过滤，后续可根据需求对vcf进行对比到已知变异的数据库，得到碱基突变是否致病，通过annovar工具进行基因文件注释。[0142]以上步骤全是分步执行，过程繁琐，考虑到实用性，通过shell脚本语言编写自动化运行程序。[0143]#！/sh/bash[0144]name＝"sample"#修改sample为数据文件名[0145]ref＝"/home/hg19fa/hg19.fasta"#参考基因文件位置[0146]snp＝"/home/hg/dbsnp_138.hg19.vcf"#snp数据库[0147]indel＝"/home/hg/mills_and_1000g_gold_standard.indels.hg19.sites.vcf"#indel数据库[0148]mkdir-p/home/wes/$name/&&cd/home/wes/$name/[0149]mkdirlog/[0150]#第一步比对转bam[0151]bwamem-t4-r'@rg\tid:foo\tpl:illumina\tlb:library\tsm:bar'$ref${name}_r1.fq.gz${name}_r2.fq.gz|samtoolsview-s-b-》${name}.bam1》log/${name}.viewtobam.log2》&1[0152]#第二步排序[0153]samtoolssort-@4-m3g-obam-o${name}.sorted.bam${name}.bam1》log/${name}.sort.log2》&1[0154]#第三步去重复[0155]gatkmarkduplicates-i${name}.sorted.bam-o${name}.sorted.markdup.bam-m${name}_duped.txt1》log/${name}.markduplicates.log2》&1[0156]gatkfixmateinformation-i${name}.sorted.markdup.bam-o${name}.sorted.markdup.fixed.bam-socoordinate1》log/${name}.fixmateinformation.log2》&1[0157]#第四步碱基质量重校[0158]gatkbaserecalibrator-r$ref-i${name}.sorted.markdup.fixed.bam‑‑known-sites$snp‑‑known-sites$indel-o${name}_baserecalibrator.table1》log/${name}.baserecalibrator.log2》&1[0159]gatkapplybqsr-r$ref-i${name}.sorted.markdup.fixed.bam-bqsr${name}_baserecalibrator.table-o${name}.sorted.markdup.fixed.bqsr.bam1》log/${name}.applybqsr.log2》&1[0160]#第五步找变异[0161]gatkhaplotypecaller-r$ref-i${name}.sorted.markdup.fixed.bqsr.bam‑‑dbsnp$snp-o${name}_raw.vcf.gz1》log/${name}.haplotypecaller.log2》&1[0162]综上所述，借助于本发明的上述技术方案，本发明的人类基因全外显检测数据的分析方法，易于推广，采用实用的自动化分析；外显子占全基因组的1％，约90％的人类致病基因变异集中在外显子区域，相比于全基因组测序，全外显子测序更加经济、高效；能够方便快速、系统性地挖掘大数据中的信息，以更全面的解析疾病，服务于生物医学和科研，为更多的生物数据分析团队，尤其是不擅长计算机的生物医学和科研工作者提供便利。[0163]以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页12当前第1页12