[更新中] NGS常见软件和数据库（肿瘤/遗传病）

本文主要介绍软件和数据库在基因检测中的应用，不追求对软件尽可能完整的介绍。
本文第1-3节介绍基因检测流程的主要环节；第4节介绍了流程搭建的常见框架。

0. 基础环境

• 操作系统：ubuntu 18.04 LTS
• （推荐）包和环境管理：miniconda（安装miniconda2）

1. 质控

trimmomatic

fastqc

2. 比对 & 比对后处理

比对：bwa

基于BWT算法将reads比对到基因组。注：IonTorrent 平台数据使用TMAP进行短序列比对。

1）编译安装

# download from https://sourceforge.net/projects/bio-bwa/
tar jxvf bwa-0.7.17.tar.bz2
cd bwa-0.7.17
make
# 添加到系统路径

2）使用方法

# 建立索引
bwa index -a bwtsw -p out_prefix /ref_path/to/reference.fasta &

#将reads比对到参考基因组
bwa mem \
    -t 4 \
    -T 0 \
    -R '@RG\tID:SP00001\tPL:illumina\tLB:WGS\tSM:SP00001' \
    /ref_path/to/out_prefix \
    S00001_L1_R1.clean.fq S00001_L1_R2.clean.fq \
    1>S00001_L1.sam 2>/dev/null &

比对后：samtools/picard/gatk

方法一：首先，使用samtools view将生成的sam格式文件转换为bam格式文件。该步骤一般通过管道接收bwa的输出，避免重复读写。实践中推荐使用set -o pipefail以便更加有效的捕获异常。其次，使用samtools sort对bam文件进行排序。
方法二：直接使用picard sort或gatk sort完成格式转换和排序。

# picard
java -jar picard.jar SortSam \
    CREATE_INDEX=true \
    INPUT=<input.bam> \
    OUTPUT=<output.bam> \
    SORT_ORDER=coordinate \
    VALIDATION_STRINGENCY=STRICT
# gatk
gatk --java-options "-Xmx12G -Djava.io.tmpdir=/path/tmp" SortSam \
    -I S00001_L1.sam \
    -O S00001_L1.sorted.bam \
    -SO coordinate \
    --CREATE_INDEX true

然后使用gatk MarkDuplicates标记PCR重复？

gatk --java-options "-Xmx12G -Djava.io.tmpdir=/path/tmp" MarkDuplicates \
    -I S00001_L1.sorted.bam \
    -O S00001_L1.sorted.marked.bam \
    -M S00001_L1.metrics

3. 变异识别（SNP/InDel）

GATK4最佳实践中不再需要Indel local realignment步骤。

GATK4最佳实践_1：碱基质量分数校正（BQSR）

基于系统误差调整碱基质量分数，可以增加下游变异识别的准确性。
碱基在reads中的位置、碱基的上下文环境、碱基原始的质量值。根据这3这个因素，首先计算出原始碱基质量中错误的分布模型，然后利用这个模型对碱基质量校正，生成新的碱基质量值。

#
gatk --java-options "-Xmx12G -Djava.io.tmpdir=/path1/tmp" BaseRecalibrator \
-R /ref_path/to/reference.fasta \
-I S00001_L1.sorted.marked.bam \
-O recal_data.table \
--known-sites /path2/bundle/dbsnp_146.hg38.vcf.gz \
--known-sites /path2/bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
--known-sites /path2/bundle/1000G_phase1.snps.high_confidence.hg38.vcf.gz

#gatk GatherBQSRReports -I ... -O ...

#
gatk --java-options "-Xmx12G -Djava.io.tmpdir=/path1/tmp" ApplyBQSR \
-R /ref_path/to/reference.fasta \
-I S00001_L1.sorted.marked.bam \
--bqsr-recal-file recal_data.table \
--add-output-sam-program-record \
-O S00001_L1.sorted.marked.BQSR.bam

GATK4最佳实践_2：Mutect2检测体细胞突变

1）

Mutect2可以进行配对和TumorOnly的分析。配对分析可以去除胚系变异和

# 检测体细胞突变
gatk --java-options "-Xmx2g" Mutect2 \
-R hg38/Homo_sapiens_assembly38.fasta \
-I tumor.bam -I normal.bam \
-tumor HCC1143_tumor -normal HCC1143_normal \  # read group id
-pon resources/chr17_pon.vcf.gz \  # 
--germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \  #人群germline
--af-of-alleles-not-in-resource 0.0000025 \    #1/(2*200k), 200k为上述vcf文件的WES样本数
-L chr17.interval_list \
-O 1_somatic_m2.vcf.gz
#--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
#-bamout 2_tumor_normal_m2.bam 

PoN（正常样本callset）定义了一个预过滤变异位点集。PoN不仅呈现了通常的germline变异位点，也呈现了测序数据引入的噪音，如测序bias或比对bias。默认情况下，出现在PoN中的变异位点不认为是tumor样本的somatic变异也不会进行重组装。参数–genotype-pon-sites可以得到PoN位点的基因型信息如果tumor位点和PoN位点变异不是精确匹配，Mutect2会重新组装此区域，在变异文件中输出此位点，并在INFO列标记为“IN_PON”

等位基因频率：在评估一个变异为体细胞突变的概率时，人群等位基因频率和参数af-of-alleles-not-in-resource都将考虑在内。

GATK最佳实践3_VQSR或hard-filtering

VQSR综合考虑了多个指标，避免了hard-filtering单一维度的缺点，因而当所分析物种存在高质量变异集、而且变异数量较多时，推荐优先使用VQSR方法完成变异过滤。

上述指标包括：Coverage, QualByDepth, FisherStrand, StrandOddsRatio, RMSMappingQuality, MappingQualityRankSumTest, ReadPosRankSumTest。

上述人类高质量SNP变异资源包括：
||

• HapMap，hapmap_3.3.hg38.vcf.gz
• Omni，1000G_omni2.5.hg38.vcf.gz，truth=true，training=false（文档中写着是true，参数建议中写着的是false。。。我就按照参数上的来了），Q12 (93.69%)
  1000G，1000G_phase1.snps.high_confidence.hg38.vcf.gz，truth=false表示VQSR考虑到在1000G数据集中的不仅包含了true variants还有false positives，training=true，Q10 (90%)
  dbSNP，dbsnp_146.hg38.vcf.gz，truth=false表示VQSR未将dbSNP数据集中的位点作为可信数据集，training=false表示不用于训练数据集，known=true表示stratify output metrics such as Ti/Tv ratio by whether variants are present in dbsnp or not，Q2 (36.90%)
  INDEL的VQSR过滤，选用的resource datasets为：

Mills，Mills_and_1000G_gold_standard.indels.hg38.vcf.gz，truth=true，training=true，Q12 (93.69%)
dbSNP，dbsnp_146.hg38.vcf.gz，truth=false，training=false，known=true，Q2 (36.90%)

4. 变异注释

ANNOVAR

1）VCF预处理（推荐）

在ANNOVAR最初开发时，很多变异检测工具（samtools,CASAVA,SOAPsnp等）输出文件的格式都不相同。因而ANNOVAR定义了一种avinput格式（chr start end ref alt op）作为后续输入格式，同时提供了convert2annovar.pl用于格式转换。
后来，VCF（Variant Call Format）逐渐成为描述变异的主流格式。convert2annovar.pl也提供了从VCF到avinput的格式转换功能。但VCF本身具有以下特点，因而推荐进行预处理。

1）从技术角度讲，VCF是一种描述基因座（locus），而非变异或基因型的格式。
2）由于一个基因座可以有多个变异，因而VCF一行原则上可以描述多个变异或基因型。
3）一行多个变异可能出现indel劫持同行其它变异的情况，见示例。
4) 理想情况下一个变异在给定参考基因组中应该有且只有一种描述方式，但如何唯一地描述目前还没有统一共识。很多用户倾向于左归一化（left-normalization），即变异起始位置应该尽可能的向左移动。然而HGVS明确指出要在cDNA坐标体系上进行左归一化，这意味着人类基因组一半基因需要右归一化。

# T->A由于indel“劫持”基因座，而被记作CTTT->CTTA
1 112240038 . CTTT CTTTT,CTTTTT,CTTA,CTT,CT,C . PASS AC=986,3,1249,3,127,3;AF=0.196885,0.000599042,0.249401,0.000599042,0.0253594,0.000599042

面对以上问题，ANNOVAR作者给出了以下VCF预处理建议，以获得更准确的注释结果。

• 直接使用左归一化
• 一行只描述一个变异，进而避免indel影响SNP的描述（尽管table_annovar.pl支持在单行注释多个变异信息）

作者基于以上两条规则对1kg, esp6000si, dbSNP等文件进行了预处理，更新后的数据库（avsnp138, avsnp142, 1000g2014oct, esp6500siv2, exac03, clinvar_20150330）已在2014年12月提供。

此外，作者推荐用户进行以下操作：

# 分割vcf行，使得每一行包含一个且只有一个变异
bcftools norm -m-both -o ex1.step1.vcf ex1.vcf.gz

# 对vcf进行左归一化
bcftools norm -f hg19.fasta -o ex1.step2.vcf ex1.step1.vcf

https://annovar.openbioinformatics.org/en/latest/articles/VCF/

2）ANNOVAR数据库下载

数据库可以通过annovar内置命令、wget、ftp等方式下载。

# 从annovar镜像下载软件开发者维护的数据库清单（包括数据名、最后更新日期和数据库大小）
perl annotate_variation.pl -downdb -buildver hg19 -webfrom annovar avdblist .
# 从annovar镜像下载refGene数据库
perl annotate_variation.pl -downdb -buildver hg19 -webfrom annovar refGene humandb/

3）ANNOVAR注释

ANNOVAR有以下三种注释格式方法。

注释方法	功能	使用	常见数据库
Gene-based	鉴别变异是否造成蛋白编码变化	-protocal g	refGene,ensGene,knownGene
Region-based	鉴别变异是否发生在特定基因组区域	-protocal r	cytoBand,
Filter-based	鉴别变异是否匹配特定数据库	-protocal f	gnomad*,exac03,1000g*,clinvar*

推荐阅读：ANNOVAR-注释软件用法详解

VEP

1）下载和安装

VEP可以通过版本库或压缩包两种形式下载。

# 版本库（包含历史版本）
git clone https://github.com/Ensembl/ensembl-vep.git
cd ensembl-vep
# 压缩包（仅特定版本）
curl -L -O https://github.com/Ensembl/ensembl-vep/archive/release/109.zip
unzip 109.zip
cd ensembl-vep-release-109/

VEP依赖gcc, g++, make, perl(5.10+), Archive::Zip, DBD::mysql, DBI，其中perl模块可以通过cpan安装。

#

然后，运行INSTALL.pl安装VEP。

# 自动安装API和缓存
perl INSTALL.pl --AUTO ac
##--AUTO a(API+Bio::DB::HTS/htslib); c(cache); p(plugins, require --PLUGINS)
##--CACHE_VERSION [version] 默认最新版。尽管可以通过该参数改变，但不能保证API和cacha适配。
##--CACHEDIR [dir] 默认安装cacha到~/.vep。如果更改，调用VEP时要使用--dir_cache。
##--DESTDIR [dir] 默认安装API模块到当前目录下的Bio文件夹。如果更改，将更该路径添加到PERL5LIB或者调用VEP时使用-I。
##--PLUGINS 逗号分隔的plugins名。
##--PLUGINSDIR [dir] 默认安装到--CACHEDIR下的Plugins子目录。如更改，调用VEP时要使用--dir_plugins。

# 手动安装
## 非必要组件，如输出json格式所需的JSON，读取bigWig所需的Bio::DB::BigFile
## 速度提升，如PerlIO::gzip, ensembl-xs
## 以下是一个安装Set::IntervalTree的示例
mkdir -p $HOME/cpanm
export PERL5LIB=$PERL5LIB:$HOME/cpanm/lib/perl5
cpanm -l $HOME/cpanm Set::IntervalTree

帮助文档

2）运行VEP

基本选项
通过--cache, --off-line, --database任一指定使用API方式。其中前两者性能较优，推荐使用。

# run VEP with default options
./vep --cache -i input.txt -o output.txt 
##--config [filename] 从配置文件读取选项配置。

VEP支持的输入格式通过--format指定，包括：ensembl, vcf, hgvs, id, region, spdi。
VCF支持的输出格式通过--vcf, --tab, --json任一指定。

Flag	Description
–assembly [name]	如果存在多个，选择assembly版本。
–force_overwrite	如果output文件存在，VEP默认抛错失败。
–stats_file [filename]	Summary stats file name（html格式）。
–no_stats	不产生stats文件。
–stats_text	产生纯文本stats文件，替换HTML。
–fork [num_forks]	Enable forking

cache常见选项

Flag	Description
–cache	启用cache
–dir [directory]	指定cache/plugin的上层目录
–dir_cache [directory]	指定cache目录
–dir_plugins [directory]	指定plugin目录
–offline	启用offline模式，禁用数据库连接，需要cache或GFF/GTF文件进行注释

其它注释源常见选项

Flag	Description
–plugin [plugin name]	plugin名，plugin需要预先安装在缓存目录下的Plugin子目录
–custom [filename]	向输出添加定制注释，文件需要tabix索引或者bigWig格式

帮助文档

3）注释源 & 定制注释

VEP的注释源包括cache、GFF/GTF和database（public+local）。其中使用cache的方式效率最高，推荐使用。

cache：①使用cache是最快和最有效的方式，在大多数情况下只有一个初始网络连接，还可以使用-offline消除所有网络连接。使用-offline时无法使用--hgvs, --format hgvs, --format id, --check_sv，其它功能如custom annotations和plugins则不受影响。②强烈建议使用与VEP版本一致的cache。③链接表格展示了人类cache（106）的内容。

# 通过脚本安装cache
perl INSTALL.pl --AUTO c
# 手动下载安装
curl -O https://ftp.ensembl.org/pub/release-109/variation/indexed_vep_cache/homo_sapiens_vep_109_GRCh38.tar.gz
tar xzf homo_sapiens_vep_109_GRCh38.tar.gz

# cache文件建索引，可节省时间（依赖Bio::DB::HTS, tabix）
perl convert_cache.pl --dir /path/cache --species homo_sapiens --version 106_GRCh38 --bgzip /path/bgzip --tabix /path/tabix --force_overwrite

定制注释：①VEP支持的格式包括GFF, GTF, VCF, BED, bigWig, BED-like等。②文件需要预先去除注释行、按按染色体和位置排序、使用bgzip压缩、使用tabix索引。③文件可以在命令行中配置，--custom Filename, Short_name, File_type, Annotation_type, Force_report_coordinates, VCF_fields，分别代表文件路径、输出的key名、文件格式（bed, gff, gtf, vcf, bigwig）、注释类型（exact, overlap）、是否输出坐标、指定INFO字段的注释内容。

4）插件

VEP可以通过用Perl编写的插件模块来增加软件的功能，如扩展、过滤和操作输出等。

# 通过脚本安装插件
perl INSTALL.pl -a p -g list

可用的VEP插件包括Pathogenicity predictions(CADD,dbNSFP,FATHMM,), Nearby features, Structural variant data等18类。通过--plugin使用插件。关于插件的工作原理和开发见本段开头链接。

5. 常见流程框架用法

5.1 脚本

# 略

5.2 makefile

targets : prerequisites
    command
    ...

target：一个对象文件，可执行文件或标签。
prerequisites：生成该 target 所依赖的文件或 target 。
command：该 target 要执行的命令。

5.3 Snakemake

Snakemake 通过一系列 rules 定义工作流。rules 包含 input, output, shell 等元素。

rule NAME:
    input: "path/to/inputfile", "path/to/other/inputfile"
    output: "path/to/outputfile", "path/to/another/outputfile"
    shell: "somecommand {input} {output}"
    ... ...

Snakemake 通过 input 和 output 确定工作流的依赖关系和执行顺序。

Snakefiles and Rules

5.4 Bpipe

bpipe 的基本单元是 stage, 使用 run 定义流程的执行顺序。

stage_one = {
  exec "shell command to execute stage one"
}

stage_two = {
  exec "shell command to execute stage two"
}

Bpipe.run {
 stage_one + stage_two
}

Bpipe Overview

5.5 Argo

argo 是基于 K8S CRD 实现的工作流工具。argo 通过一个 yaml 格式的配置文件来定义工作流。

apiVersion: argoproj.io/v1alpha1
kind: Workflow
metadata:
  generateName: steps-
spec:
  entrypoint: hello-hello-hello
  templates:
  - name: hello-hello-hello   # template-1
    steps:
    - - name: hello1             ##
        template: whalesay
        arguments:
          parameters:
          - name: message
            value: "hello1"
    - - name: hello2a                ###
        template: whalesay
        arguments:
          parameters:
          - name: message
            value: "hello2a"
      - name: hello2b                ### 
        template: whalesay
        arguments:
          parameters:
          - name: message
            value: "hello2b"
  - name: whalesay            # template-2
    inputs:
      parameters:
      - name: message
    container:
      image: docker/whalesay
      command: [cowsay]
      args: ["{{inputs.parameters.message}}"]

Argo Exam

5.6 Cromwell (WDL)

workflow example{
  File firstInput
  scatter (idx in range(3)) {
    call stepa { input: in=firstInput }
  }
  call stepb { input: in=stepa.out }
  call stepc { input: in=stepb.out }
  output {
    File result = stepc.out
  }
}
task stepa {
  File in
  command { cat ${in} > outputa.txt && cat outputa.txt}
  output { File out = "outputa.txt" }
  runtime {
    docker: "swr.cn-north-1.myhuaweicloud.com/op_svc_gcs_container/cromwell:1.5.8"
    cpu: "1"
    memory: "2G"
  }
}
task stepb {
  Array[File] in
  command { cat ${write_lines(in)} > outputb.txt && cat outputb.txt}
  output { File out = "outputb.txt" }
  runtime {
    docker: "swr.cn-north-1.myhuaweicloud.com/op_svc_gcs_container/cromwell:1.5.8"
    cpu: "1"
    memory: "2G"
  }
}
task stepc {
  File in
  command { cat ${in} > outputc.txt && cat outputc.txt}
  output { File out = "outputc.txt" }
  runtime {
    docker: "swr.cn-north-1.myhuaweicloud.com/op_svc_gcs_container/cromwell:1.5.8"
    cpu: "1"
    memory: "2G"
  }
}

Cromwell 示例