全基因组(WGS) 思维导图
  1. 首页
  2. 生信组学教程
  3. 全基因组(WGS)
  4. 从Fastq到VCF
切换视频源：
GVCF 模式

从Fastq到VCF

最后发布时间 : 2023-08-05 15:49:18 浏览量 :

生信小木屋

学习资料:

软件安装:

conda create -n bwa bwa
conda activate bwa
conda install samtools
conda install bcftools
conda  create -n gatk4 gatk4==4.4.0.0

生信小木屋

比对

使用bwa建立参考基因组的index，并进行比对

bwa index output/index/ref.fasta

bwa mem -t 10 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:mother' output/index/ref.fasta \
    2-germline/fastq/mother_1.fastq \
    2-germline/fastq/mother_2.fastq > output/mapping/mother.sam

samtools view -bS  -@ 10 output/mapping/mother.sam -o  output/mapping/mother.bam
samtools  sort output/mapping/mother.bam -@ 10 -o   output/mapping/mother.sorted.bam

去除PCR重复

比对得到bam文件后，使用gatk的子命令MarkDuplicates标记PCR重复，并使用samtools排序、建立索引。

gatk MarkDuplicates -I output/mapping/mother.sorted.bam \
    -O output/mapping/mother.sorted.markdup.bam  \
    -M output/mapping/mother.markdup.metrics.txt
	
samtools index output/mapping/mother.sorted.markdup.bam

生信小木屋

碱基质量矫正（BQSR）

samtools faidx output/index/ref.fasta
tabix -p vcf vcf-resource/Homo_sapiens_assembly38.known_indels.vcf.gz
gatk CreateSequenceDictionary -R output/index/ref.fasta -O output/index/ref.dict

vcf文件需要建立索引

gatk BaseRecalibrator \
    -I F11667.new.sort.markdup.bam  \
    -R out/ref/hg38.fa \
    --known-sites Homo_sapiens_assembly38.known_indels.vcf.gz \
    -O F11667.recal_data.table

gatk ApplyBQSR \
   -R out/ref/hg38.fa \
   -I F11667.new.sort.markdup.bam \
   --bqsr-recal-file F11667.recal_data.table \
   -O F11667.new.sort.markdup.BQSR.bam

这里包含三个步骤

  • 第一步：BaseRecalibrator，这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件（F11667.recal_data.table）；
  • 第二步：这一步利用第一步得到的校准表文件（F11667.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。

注意，因为BQSR实际上是为了（尽可能）校正测序过程中的系统性错误，因此，在执行的时候是按照不同的测序lane或者测序文库来进行的，这个时候@RG信息（BWA比对时所设置的）就显得很重要了，算法就是通过@RG中的ID来识别各个独立的测序过程，这也是我开始强调其重要性的原因。

变异检测

生信小木屋

事实上，这是目前所有WGS数据分析流程的一个目标——获得样本准确的变异集合。这里变异检测的内容一般会包括：SNP、Indel，CNV和SV等，这个流程中我们只做其中最主要的两个：SNP和Indel。我们这里使用GATK HaplotypeCaller模块对样本中的变异进行检测，它也是目前最适合用于对二倍体基因组进行变异（SNP+Indel）检测的算法。

HaplotypeCaller和那些直接应用贝叶斯推断的算法有所不同，它会先推断群体的单倍体组合情况，计算各个组合的几率，然后根据这些信息再反推每个样本的基因型组合。因此它不但特别适合应用到群体的变异检测中，而且还能够依据群体的信息更好地计算每个个体的变异数据和它们的基因型组合。

一般来说，在实际的WGS流程中对HaplotypeCaller的应用有两种做法，差别只在于要不要在中间生成一个gVCF：

  • 直接进行HaplotypeCaller，这适合于单样本，或者那种固定样本数量的情况，也就是执行一次HaplotypeCaller之后就老死不相往来了
  • 每个样本先各自生成gVCF，然后再进行群体joint-genotype。gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了

我们首先以第一种直接进行HaplotypeCaller做法为例：
调用gatk子命令CreateSequenceDictonary为参考序列生成一个.dict文件

OK，现在我们就可以进行变异检测了，同样使用GATK 4.0的HaplotypeCaller模块来完成。由于我们只有一个样本，要完成这个工作其实很简单，直接输入比对文件和参考序列就行了。

gatk HaplotypeCaller \
    -R out/ref/hg38.fa \
    -I F11667.new.sort.markdup.BQSR.bam \
    -O F11667.new.sort.markdup.BQSR.vcf  \

    -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \
    -A QualByDepth \ 
    -A RMSMappingQuality \ 
    -A MappingQualityRankSumTest \ 
    -A ReadPosRankSumTest \ 
    -A FisherStrand \ 
    -A StrandOddsRatio \ 
    -A Coverage
  • -D（[--dbsnp](https://gatk.broadinstitute.org/hc/en-us/articles/13832687299739-HaplotypeCaller#--dbsnp)）参数输入的dbSNP同样可以再GATK bundle目录中找到，这份文件汇集的是目前几乎所有的公开人群变异数据集；
  • -A（[--annotation](https://gatk.broadinstitute.org/hc/en-us/articles/13832687299739-HaplotypeCaller#--annotation)）一个或多个特定的注释要添加到变体调用中

另外，由于我们的例子只有一个样本因此只输入一个BAM文件就可以了，如果有多个样本那么可以继续用-I参数输入：

java -jar GenomeAnalysisTK.jar \
    -T HaplotypeCaller \
     -R reference.fasta \
     -I sample1.bam [-I sample2.bam ...] \
     ...

以上的命令是直接对全基因组做变异检测，这个过程会消耗很长的时间，通常需要几十个小时甚至几天。

然而，基因组上各个不同的染色体之间其实是可以理解为相互独立的（结构性变异除外），也就是说，为了提高效率我们可以按照染色体（或者染色体不同区域）一条条来独立执行这个步骤，最后再把结果合并起来就好了，这样的话就能够节省很多的时间。下面我给出一个按照染色体区分的例子：

gatk HaplotypeCaller \
    -R out/ref/hg38.fa \
    -I F11667.new.sort.markdup.BQSR.bam \
    -O F11667.new.sort.markdup.BQSR.vcf  \
    -L 1 \

    -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \
    -A QualByDepth \ 
    -A RMSMappingQuality \ 
    -A MappingQualityRankSumTest \ 
    -A ReadPosRankSumTest \ 
    -A FisherStrand \ 
    -A StrandOddsRatio \ 
    -A Coverage

注意到了吗？其它参数都没任何改变，就只增加了一个 -L 参数，通过这个参数我们可以指定特定的染色体（或者基因组区域-L 20:10,000,000-10,200,000）！我们这里指定的是 1 号染色体，有些地方会写成chr1，具体看human.fasta中如何命名，与其保持一致即可。其他染色体的做法也是如此，就不再举例了。最后合并：

gatk CombineVariants \
 -R /path/to/human.fasta \
 --genotypemergeoption UNSORTED \
 --variant sample_name.HC.1.vcf \
 --variant sample_name.HC.2.vcf \
 ...
 --variant sample_name.HC.MT.vcf \
 -o sample_name.HC.vcf
GVCF 模式