目录

  • 0. 基础环境
  • 1. 质控
    • trimmomatic
    • fastqc
  • 2. 比对 & 比对后处理
    • bwa
    • TMAP
  • 3. 变异识别（SNP/InDel）
    • samtools
    • gatk
    • bedtools
  • 4. 变异注释
    • ANNOVAR
      • 1）VCF预处理（推荐）
      • 2）ANNOVAR数据库下载
      • 3）ANNOVAR注释
    • snpEff
    • VEP
  • 5. 常见流程框架用法
    • 5.1 脚本
    • 5.2 makefile
    • 5.3 Snakemake
    • 5.4 Bpipe
    • 5.5 Argo
    • 5.6 Cromwell (WDL)
  • 附录
  • VCF文件格式

本文主要介绍软件和数据库在基因检测中的应用，不追求对软件尽可能完整的介绍。

本文第1-3节介绍基因检测流程的主要环节；第4节介绍了流程搭建的常见框架。

0. 基础环境

  • 操作系统：ubuntu 18.04 LTS
  • （推荐）包和环境管理：miniconda（安装miniconda2）
################################
# （读者请忽略）因为我的电脑显示有点问题，此处调了一下分辨率。
################################
## cvt 1600 900
## xrandr --newmode "1600x900_60.00"  118.25  1600 1696 1856 2112  900 903 908 934 -hsync +vsync
## xrandr --addmode DP-1 "1600x900_60.00"
## xrander --output DP-1 --mode 1600x900_60.00

1. 质控

trimmomatic

fastqc

2. 比对 & 比对后处理

该环节主要是将 reads 比对到参考基因组上。

bwa

  1. 简介
    建立 index。
    基于 BWT 算法，将 reads 比对到参考基因组。
    最新版本 bwa-mem2，Intel实验室对计算效率进行了优化。
  2. 编译安装
# download from https://sourceforge.net/projects/bio-bwa/
tar jxvf bwa-0.7.17.tar.bz2
cd bwa-0.7.17
make
# 添加到系统路径# 如果ubuntu没有安装make等命令，运行以下命令安装：
## sudo apt install make
## sudo apt install gcc
## sudo apt-get install zlib1g-dev
  1. 常见用法

TMAP

用于 IonTorrent 平台数据的短序列比对。

3. 变异识别（SNP/InDel）

samtools

使用 miniconda 安装

gatk

bedtools

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-注释软件用法详解

snpEff

VEP

5. 常见流程框架用法

5.1 脚本

# 略

5.2 makefile

targets : prerequisitescommand...

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-hellotemplates:- name: hello-hello-hello   # template-1steps:- - name: hello1             ##template: whalesayarguments:parameters:- name: messagevalue: "hello1"- - name: hello2a                ###template: whalesayarguments:parameters:- name: messagevalue: "hello2a"- name: hello2b                ### template: whalesayarguments:parameters:- name: messagevalue: "hello2b"- name: whalesay            # template-2inputs:parameters:- name: messagecontainer:image: docker/whalesaycommand: [cowsay]args: ["{{inputs.parameters.message}}"]

Argo Exam

5.6 Cromwell (WDL)

workflow example{File firstInputscatter (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 incommand { 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] incommand { 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 incommand { 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 示例

附录

VCF文件格式

VCF最初由千人基因组项目开发和使用，其规格和扩展目前由GA4GH负责。以v4.2为例说明格式规格。

  • 元信息：##key=values
    • ##fileformat=VCFv4.2
    • ##INFO=<ID=ID,Number=number,Type=type,Description=“description”,…>
    • ##FILTER=<ID=ID,Description=“description”>
    • ##FORMAT=<ID=ID,Number=number,Type=type,Description=“description”>
    • ##ALT=<ID=type,Description=“description”>
    • ##assembly=url
    • ##contig=<ID=ctg1,URL=ftp://somewhere.org/assembly.fa,…>
    • …
  • 标题行：包括8个必需列，以制表符分隔。
    1. #CHROM
    2. POS
    3. ID
    4. REF
    5. ALT
    6. QUAL
    7. FILTER
    8. INFO
    9. （当基因型数据存在时）FORMAT sample_IDs
  • 数据行：包括8个必需字段，以制表符分割。缺失值以点号填充。
    1. CHROM：一个来自参考基因组的ID，指向##assembly的一个contig。一个特定CHROM的所有条目应该在VCF文件中形成一个连续区块。
    2. POS：参考位置（1-based）。在每个CHROM中，位置按数字递增排序。允许有多个具有相同POS的记录。端粒通过使用位置0或N+1来表示，其中N是相应的染色体或contig的长度。
    3. ID：分号分隔的唯一标识符列表。如果在dbSNP中，鼓励使用rs编号。
    4. REF：参考碱基。[ATGCN]+。POS值是字符串中第一个碱基的位置。①REF和ALT必须包括indel之前的碱基（同时体现在POS中），除非indel发生在第一个位置，此时必须包括indel之后的碱基；②对于复杂的替换或其他事件，例如所有等位基因在其字符串中至少有一个碱基，不需要（但允许）这种碱基填充；③如果任何一个ALT是一个符号化的等位基因（一个带角括号的字符串），则填充是必需的。
    5. ALT：变异碱基。由[ACTG*]组成、逗号分隔的变异碱基列表；或一个带角括号的字符串；或一个断尾替换字符串。'*'表示该等位基因由于上游的删除而丢失。如果没有变异碱基，则使用缺失值。
    6. QUAL：ALT的phred质量分数。
    7. FILTER：过滤状态。通过所有的过滤，则为PASS；否则显示分号分隔的、未通过过滤的原因代码。
    8. INFO：附加信息。常见的有AC, AF, DP, MQ等。
    9. FORMAT
    10. sample_IDs

https://github.com/samtools/hts-specs/blob/master/VCFv4.2.pdf