其他分享
首页 > 其他分享> > 人类WGS流程

人类WGS流程

作者：互联网

目录

1.测序

DNA样品超声随机打断，片段纯化，末端修复，3’端加 A碱基，两端加文库接头。

带接头的 DNA 模板进行 PCR 扩增富集，对DNA 文库浓度定量，质控合格后上机测序。现在国内二代平台一般使用 Illumina HiSeq 或BGIseq/MGIseq系列平台。

测序得到的原始图像数据，经 Illumina 碱基识别软件（如Base Calling）转化为原始序列raw reads（paired-end reads），以 FASTQ 格式存储，即 raw data。

2.流程概述

raw data中包含接头（adapter）序列、测序质量很低的碱基、未测出的碱基N，因此首先需要对raw data进行过滤得到  clean data。

然后，使用比对软件BWA（Burrows-Wheeler Aligner) 将每个样品的  clean data  比对到人的参考基因组GRCh37/hg19，得到  BAM  格式的最初比对结果文件。为了确保变异检测的精准性，按照 GATK（Genome Analysis Toolkit） 官网推荐的最优变异检测分析流程做分析。使用 Picard 工具将比对结果去除重复的reads，使用 GATK做局部重比对（local realignment）和碱基质量值重校正BQSR（base quality recalibration）。

最后，基于比对结果，对每个样品的测序深度、覆盖度、比对率等评价指标进行统计。

分析软件及QC：

下面分步说明。

3.raw data过滤

原始数据过滤方法比如：

过滤后得到 clean data ，并对测序数据进行统计，包括测序reads数量、数据产量、质量值分布等。

4.比对

使用BWA将所有的 clean reads 比对到人参考基因组上。分别对每一条 lane 的测序数据进行比对，并在比对结果中添加 read group ID。

bwa mem -M -R '@RG\tID:GroupID\tSM:SampleID\tPL:illumina\tLB:libraryID' \
ucsc.hg19.fasta read1.fq.gz read2.fq.gz > aligned_reads.sam

对比对结果SAM文件按比对位置排序，并转换为BAM文件：

java -jar picard-tools/SortSam.jar  \
  I=aligned_reads.sam \
  O=aligned_reads.sorted.bam \
  SORT_ORDER=coordinate

由于在建库过程中的  PCR引入很多相同的 DNA 片段，出现duplicate reads 干扰变异检测结果，因此需要标记去除  duplicate reads 。

java -jar picard-tools/MarkDuplicates.jar \
   I=aligned_reads.sorted.bam \
   O=aligned_reads.sorted.dedup.bam \
   METRICS_FILE=metrics.txt
 
java -jar BuildBamIndex.jar \
   I=aligned_reads.sorted.dedup.bam

5.局部重比对

由于 InDels 会影响 BWA 在其附近区域的比对结果的准确性，会出现错误的比对位置，因此需要做局部重比对。它分为两步：一是确定重新比对的区域位置，二是确定最优的一致性序列，并重新调整 reads 的比对信息。

java -jar GenomeAnalysisTK.jar \
   -T RealignerTargetCreator \
   -R gatk_ref/ucsc.hg19.fasta \
   -o indels_religner.intervals \
   -known 1000G_phase1.indels.hg19.vcf \
   -known Mills_and_1000G_gold_standard.indels.hg19.vcf
 
java -jar GenomeAnalysisTK.jar \
   -T IndelRealigner \
   -R ucsc.hg19.fasta \
   -I aligned_reads.sorted.dedup.bam \
   -targetIntervals indels_religner.intervals \
   -known 1000G_phase1.indels.hg19.vcf  \
   -known Mills_and_1000G_gold_standard.indels.hg19.vcf \
   -o aligned_reads.sorted.dedup.realigned.bam

6.碱基质量校正（BQSR）

测序仪器带来的各种系统性误差会引起碱基质量值过高或者过低，因此需要做碱基质量值校正，以得到较精确的碱基质量值，提高变异检测的准确度。

java -jar GenomeAnalysisTK.jar \
   -T BaseRecalibrator \
   -R gatk_ref/ucsc.hg19.fasta \
   -I aligned_reads.sorted.dedup.realigned.bam \
   -knownSites dbsnp_151.hg19.vcf \
   -knownSites Mills_and_1000G_gold_standard.indels.hg19.vcf \
   -knownSites 1000G_phase1.indels.hg19.vcf \
   -o recal.table
 
java -jar GenomeAnalysisTK.jar \
   -T PrintReads \
   -R gatk_ref/ucsc.hg19.fasta \
   -I aligned_reads.sorted.dedup.realigned.bam \
   -BQSR recal.table \
   -o aligned_reads.sorted.dedup.realigned.recal.bam

7.变异检测（SNP和InDel）

使用GATK的 HaplotypeCaller 工具同时检测 SNPs 和 InDels 。其原理是，在有变异信号的区域做单体型的局部 denovo 组装。原始的变异集合以 VCF 格式存储，它包括所有潜在的变异位点。

java -jar GenomeAnalysisTK.jar \
   -T HaplotypeCaller \
   -R gatk_ref/ucsc.hg19.fasta --genotyping_mode DISCOVERY \
   -I aligned_reads.sorted.dedup.realigned.recal.bam \
   -L CallVariantRegion/ex_region.sort.bed \
   -o raw_variants.vcf \
   -stand_call_conf 30 \
   -stand_emit_conf 10 \
   -minPruning 3

8.变异过滤

SNP 和 InDel 的原始变异集过滤可采用基于机器学习算法的变异质量值校正（VQSR）的方法。 GATK VQSR 是使用已知的、高质量的变异集合作为训练集及真实集合，并建立一个预测模型来过滤掉假变异。输出的 VCF 结果中标记为 "PASS" 的 SNP 和 InDel 就是通过了过滤条件的、可信的变异集。

对于 SNP 校正策略，使用以下数据集和特征来训练模型。(a) 训练集：HapMap, Omni2.5M 芯片数据和千人计划得到的高可信的 SNP 位点。(b) 特征：Coverage (DP), Quality/depth (QD), Fisher test on strand bias(FS), Odds ratio for strand bias (SOR), Mapping quality rank sum test (MQRankSum), Read position rank sum test (ReadPosRankSum), RMS mapping quality (MQ)。参数如下：

java -jar GenomeAnalysisTK.jar -T SelectVariants \
   -R gatk_ref/ucsc.hg19.fasta \
   -V raw_variants.vcf -selectType  SNP  \
   -o raw_snps.vcf
 
java -jar GenomeAnalysisTK.jar  -T VariantRecalibrator \
   -R gatk_ref/ucsc.hg19.fasta -input raw_snps.vcf \
   -resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg19.vcf \
   -resource:omni,known=false,training=true,truth=true,prior=12.0 1000G_omni2.5.hg19.vcf \
   -resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg19.vcf \
   -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.hg19.vcf \
   -an DP -an QD -an FS -an SOR -an MQ -an MQRankSum -an ReadPosRankSum \
   -mode  SNP  \
   -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \
   -recalFile recalibrate_SNP.recal \
   -tranchesFile  recalibrate_SNP.tranches \
   -rscriptFile recalibrate_SNP_plots.R 
   
java -jar GenomeAnalysisTK.jar -T ApplyRecalibration \ 
    -R gatk_ref/ucsc.hg19.fasta \ 
    -input raw.snp.vcf \ 
    -mode  SNP  \ 
    --ts_filter_level 99.0 \ 
    -recalFile recalibrate_SNP.recal \ 
    -tranchesFile recalibrate_SNP.tranches \ 
    -o filtered_snp.vcf

对于 InDel 的校正策略，使用以下数据集和特征来训练模型。(a) 训练集：Mills 1000G gold standard indel set. (b) 特征: Coverage (DP), Quality/depth (QD), Fisher test on strand bias(FS), Odds ratio for strand bias (SOR), Mapping quality rank sum test (MQRankSum), Read position rank sum test (ReadPosRankSum)。参数如下：

java -jar GenomeAnalysisTK.jar -T SelectVariants \
   -R gatk_ref/ucsc.hg19.fasta \
   -V raw_variants.vcf -selectType INDEL \
   -o raw_indels.vcf
 
java -jar GenomeAnalysisTK.jar -T VariantRecalibrator \
   -R gatk_ref/ucsc.hg19.fasta -input raw_indels.vcf \
   -resource:mills,known=true,training=true,truth=true,prior=12.0 Mills_and_1000G_gold_standard.indels.hg19.vcf
   -an QD -an DP -an FS -an SOR -an MQRankSum -an ReadPosRankSum -mode INDEL \
   -tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 --maxGaussians 4 \
   -recalFile recalibrate_INDEL.recal \
   -tranchesFile recalibrate_INDEL.tranches  \
   -rscriptFile recalibrate_INDEL_plots.R
   
java -jar GenomeAnalysisTK.jar -T ApplyRecalibration \ 
   -R gatk_ref/ucsc.hg19.fasta \ 
   -input raw_indels.vcf \ 
   -mode INDEL \
   --ts_filter_level 99.0 \
   -recalFile recalibrate_INDEL.recal \
   -tranchesFile recalibrate_INDEL.tranches  \
   -o filtered_indel.vcf

9.CNV检测

可以用基于深度信号算法的 CNVnator来检测拷贝数变异，其原理是，将基因组分割成互不重叠的、相等长度的窗口，再标准化每个窗口中比对上的 reads 数目作为深度信号来检测拷贝数变异。
流程使用标准参数及设置，并选取 100 bp 的窗口长度。包括以下几步：

cnvnator -root out.root -tree sample.bam -unique
cnvnator -root out.root -his 100 -d hg19_chr_fa_dir
cnvnator -root out.root -stat 100
cnvnator -root out.root -partition 100
cnvnator -root out.root -call 100 > sample. CNV

10.SV检测

结构变异有5 种类型：染色体间易位 (CTX)；染色体内易位 (ITX)；倒位 (INV)； 缺失 (DEL)；插入 (INS)。

可使用 Breakdancer检测SV ，采用默认参数设置，主要是基于双末端的 reads 之间分离的距离及比对的方向，使用成对不一致的算法来检测的 。

也可使用 CREST来检测结构变异， 该方法使用 soft-clipping reads 信息，并利用 CAP3 和 BLAT来执行“组装-定位-查找-组装-比对” 这一过程来鉴定结构变异的两个断点位置。

breakdancer_max sample.cfg >sample.out
extractSClip -i sample.bam \
   --ref_genome ucsc.hg19.fasta -o outdir

CREST -f sample.bam.cover -d sample.bam \
   --ref_genome ucsc.hg19.fasta -t ucsc.hg19.fasta.2bit \
   --cap3 /path/to/cap3 --blatclient /path/to/gfClient \
   --blat /path/to/blat -o outdir

11.变异结果注释

得到高可信、高质量的变异集之后，使用 SnpEff软件对变异结果做注释，注释主要分为两个方面。

12.常见问题

1）适合样本？测序深度？

大部分都可，FFPE不推荐，因为DNA存在降解，CNV/SV等结构性变异假阳性高。
一般30X，如果是找较大的结构变异、低丰度突变，需要50X以上；群体重测序可10X。

2）如何找候选变异？

关注非同义突变、剪接突变、移码突变。去除千人基因组数据库中 MAF >=1% 的变异；去除 NHLBI-ESP6500 European American 群体数据库中 MAF >=1% 的变异；去除 NHLBI-ESP6500 African American 群数据库中 MAF >=1%的变异；推测变异的致病性。利用 SIFT/PolyPhen2/Mutation assessor/Condel/FATHMM 进行打分，预测某个变异和氨基酸置换是否影响蛋白功能。如果 score<=0.05 或 PolyPhen2>=0.909 或 MA score>=1.9 或 Condel = deleterious 或 FATHMM=deleterious，就推测该变异可能是有害变异。

尽量避免选择的变异。 InDels 附近的变异；串联重复序列区域的变异；同源序列区域的变异；等位基因不平衡（等位基因频率 <0.25 或 >0.75）的杂合性变异位点。

3）如何验证变异？

Ref：
https://www.haplox.cn/faq.html
--- 整理来自网上，侵删---

标签：vcf,变异,测序,流程,jar,reads,WGS,人类,hg19
来源： https://www.cnblogs.com/jessepeng/p/13062445.html