人类WGS流程

2023-10-11 21:28:52

目录

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：

  • GATK 的 HaplotypeCaller来检测基因组变异（ SNP  和  InDel）；
  • 利用变异质量值校正（VQSR）对原始的变异检测结果进行过滤；
  • 使用基于深度信号方法的 CNVnator 检测拷贝数变异（ CNV ）；
  • 使用 Breakdancer或者 CREST 检测结构变异（ SV ）；
  • 使用 SnpEff软件对变异结果进行注释及影响预测。

下面分步说明。

3.raw data过滤

原始数据过滤方法比如：

  • 过滤掉 reads 的接头序列；
  • 当SE reads 中低质量碱基数超过该条 reads 的 50% 时，需要去除这一对 reads。
  • 当SE reads 中 N 碱基的含量超过该条 reads 的 10% 时，需要去除这一对 reads。

过滤后得到 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软件对变异结果做注释，注释主要分为两个方面。

  • 基于基因注释：SNPs / InDels  是否导致蛋白编码或氨基酸发生变化。

  • 基于数据库注释：如变异是否出现在 dbSNP数据库中，或者在千人基因组项目中的次等位基因频率 ( MAF ) 是否低于 1%，或者确定编码区非同义突变的  SNPs  的 SIFT 值是否小于 0.05，或者变异的保守性预测值 GERP++ 值是否大于 2，或其它的注释信息等。

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）如何验证变异？

  • SNPs可通过PCR；SNP分型；
  • 小片段InDel可通过PCR；
  • CNVs可通过Real-time PCR，并根据CT值估算不同个体拷贝数变化倍数；
  • 小的SVs可通过PCR和测序，大的SVs需要亚显微方法，如FISH。

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

上一篇：java8判断今天是不是本月最后一天


下一篇：Java 使用SimpleDateFormat和DateTimeFormatter格式化日期时间的方法及示例代码