基因测序VCF文件的解析与处理

yoyo 12/6/2021 基因测序VCFGWAShtslibbcftools

# 1. 前言

# 1.1 基因测序技术

基因测序技术也称作DNA测序技术，即获得目的DNA片段碱基排列顺序的技术，获得目的DNA片段的序列是进一步进行分子生物学研究和基因改造的基础。

1977年，Walter Gilbert和Frederick Sanger发明了第一台测序仪，并应用其测定了第一个基因组序列，噬菌体X174，全长5375个碱基。由此开始，人类获得了探索生命遗传本质的能力，生命科学的研究进入了基因组学的时代，到至今为止的四十年时间内，测序技术已取得了相当大的发展，从第一代发展到了第三代测序技术。Sanger所发明的测序方法被称为第一代测序技术，该技术直到现在依然被广泛使用，但是其一次只能获得一条长度在700～1000个碱基的序列，无法满足现代科学发展对生物基因序列获取的迫切需求。

高通量测序 (High-Throughput Sequencing, HTS) 是对传统Sanger测序的革命性变革，其解决了一代测序一次只能测定一条序列的限制，一次运行即可同时得到几十万到几百万条核酸分子的序列，因此也被称为第二代测序技术。

第二代测序技术虽然测序的通量大大增加，但是其获得单条序列长度很短，想要得到准确的基因序列信息依赖于较高的测序覆盖度和准确的序列拼接技术，因此最终得到的结果中会存在一定的错误信息。因此，科研人员又发明了第三代测序技术也称为单分子测序技术，该技术在保证测序通量的基础上，对单条长序列进行从头测序，能够直接得到长度在数万个碱基的核酸序列信息。

基因测序技术的发展历史

# 1.2 全基因组关联分析

GWAS(Genome-wide association study)，即全基因组关联分析，是指在全基因组范围内找出存在的序列变异，即单核苷酸多态性(SNP)，从中筛选出与疾病相关的SNPs。全基因组关联研究是一种检测特定物种中不同个体间的全部或大部分基因，从而了解不同个体间的基因变化有多大的一种方法。不同的变化带来不同的性状，如各种疾病的不同。

# 2. 解析基因变异记录VCF文件

# 2.1 文件格式说明

VCF是用于描述SNP（单个碱基上的变异），INDEL（插入缺失标记）和SV（结构变异位点）结果的生物基因数据文件，通过基因测序技术获得。

文件的开头是一堆以“##”开始的注释行，包含了文件的基本信息。然后是以“#”开头的一行，共9+n个部分，前九部分标注的是后面行每部分代表的信息，相当于表头。后面部分是样本名称，可以有多个。注释行结束后是具体的突变信息，每一行分为9+n个部分，每部分之间用制表符（‘\t’）分隔。

VCF基因测序文件示例

通常会对样本压缩成.vcf.gz格式，这个我们下面使用PyVCF库可以直接解析，无需解压。如果需要解压取样的话，可以使用如下命令：

$ gunzip test.vcf.gz                             // 解压.vcf.gz得到.vcf文件，并删除原始.vcf.gz文件
$ cat test.vcf | head -n 10000 > test_1w.vcf     // 对.vcf文件进行取样

1
2

注：使用gunzip解压成功后，会自动把原始的.vcf.gz文件删掉，如果需要保留，可以提前拷贝一份。

# 2.2 使用PyVCF库解析文件

PyVCF：通过正则表达式把vcf文件信息转换成结构化的信息，简化了vcf文件的处理过程，方便后续提取相关参数及处理。

项目地址：https://github.com/jamescasbon/PyVCF (opens new window)
安装依赖：直接安装PyVCF可能会遇到error in setup command: use_2to3 is invalid.报错，因为在setuptools 58之后的版本已经废弃了use_2to3。
```
$ pip3 install setuptools==57.5.0 
$ pip3 install PyVCF==0.6.8
```
1
2

PyVCF库使用介绍

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    print(record)
    
>>> 输出内容   
Record(CHROM=NC_006583.3, POS=548048, REF=G, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548057, REF=T, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548071, REF=C, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548075, REF=G, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548076, REF=T, ALT=[<NON_REF>])

1
2
3
4
5
6
7
8
9
10
11
12

调用vcf.Reader类处理vcf文件，vcf文件信息就被保存到vcf_reader中了。它是一个可迭代对象，它的迭代元素都是一个_Record对象的实例，保存着非注释行的一行信息，即变异位点的具体信息。通过它，我们可以很轻易地得到位点的详细信息。

vcf.gz文件可以指定读哪个染色体的哪个位置，但需要目录处有对应的.vcf.gz.tbi文件并且安装pysam库。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

vcf_reader.fetch('chr_1', 2, 5)
for record in vcf_reader:
    print(record)

1
2
3
4
5
6
7

详见：https://pyvcf.readthedocs.io/en/latest/API.html#vcf-reader (opens new window)

[1] _Record对象-位点信息的储存形式

class vcf.model._Record(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, sample_indexes, samples=None)

下面分别解释一下 CHROM、POS、ID、REF、ALT、QUAL、FILTRE、INFO、FORMAT 等9部分的含义。

CHROM：染色体名称，类型为str。
POS：位点在染色体上的位置，类型为int。
ID：一般是突变的rs号，类型为str。如果是‘.’，则为None。
REF：参考基因组在该位点上的碱基，类型为str。
ALT：在该位点的测序结果，类型为list，代表了突变的几种类型。
QUAL：该位点的测序质量，类型为int或float。
FILTER：过滤信息。将FILTER列按分号分隔形成的字符串列表，类型为list。如果未给出参数则为None。
INFO：该位点的一些测试指标。将‘=’前的参数作为键，后面的参数作为值，构建成的字典，类型为dict。
FORMAT：基因型信息。保存vcf的FORMAT列的原始形式，类型为str。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    print(record)
    print("=====")
    print(type(record.CHROM), record.CHROM)
    print(type(record.POS), record.POS)
    print(type(record.ID), record.ID)
    print(type(record.REF), record.REF)
    print(type(record.ALT), record.ALT)
    print(type(record.QUAL), record.QUAL)
    print(type(record.FILTER), record.FILTER)
    print(type(record.INFO), record.INFO)
    print(type(record.FORMAT), record.FORMAT)
    break
    
>>> 输出内容
Record(CHROM=NC_006583.3, POS=1, REF=N, ALT=[<NON_REF>])
=====
<class 'str'> NC_006583.3
<class 'int'> 1467
<class 'NoneType'> None
<class 'str'> G
<class 'list'> [<NON_REF>]
<class 'NoneType'> None
<class 'NoneType'> None
<class 'dict'> {'END': 1471}
<class 'str'> GT:DP:GQ:MIN_DP:PL

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

除了这些基本属性之外，_Record对象还有一些其他属性：

samples：把FORMAT信息作为键，后面对应的信息做为值，构建成的字典（CallData对象），以及sample名称，这两个值组成一个Call对象，共同构成samples的一个元素。这样就把sample和基因型信息给关联起来，按下标访问每一个Call对象。samples类型为list。
start：突变开始的位置。
end：突变结束的位置。
alleles：该位点所有的可能情况，由REF和ALT参数组成的列表（REF类型是str，ALT参数是_AltRecord对象的子类实例），类型是list。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    # 按下标访问Call，按.sample访问sample，按键访问FORMAT对应信息
    print(record.samples, '\n', record.samples[0].sample, '\n', record.samples[0]['GT'])  
    print(record.start, record.POS, record.end)
    print(record.REF, record.ALT, record.alleles)
    break
    
>>> 输出内容
[Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))] 
A1 
0/0
0 1 1
N [<NON_REF>] ['N', <NON_REF>]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

_Record对象方法：

对象之间比较大小方法：根据染色体名称和位置信息比较。
迭代方法：对samples里的元素进行迭代。
字符串方法：只返回CHROM，POS，REF，ALT四列信息。
genotype(name)方法，和samples按下标访问不同，这个方法提供按sample名称进行访问的功能。
add_format(fmt)， add_filter(flt)， add_info(info, value=True)：给相应的属性添加元素。
get_hom_refs()：拿到samples中该位点未突变的所有sample，返回列表。
get_hom_alts()：拿到samples中该位点100%突变的所有sample，返回列表。
get_hets()：拿到samples中该位点基因型为杂合的所有sample，返回列表。
get_unknown()：拿到samples中该位点基因型未知的所有sample，返回列表。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

# 按染色体名称和位置进行比较
record = next(vcf_reader)
record2 = next(vcf_reader)
print(record < record2)

# 按samples列表进行迭代
for i in record:
    print(i)
    
# 按sample名字进行访问
print(record.genotype('A1'))

>>> 输出内容
True
Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))
Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

_Record对象还有很多有用的方法属性：

num_called：该位点已识别的sample数目。
call_rate：已识别的sample数目占sample总数的比例。
num_hom_ref，num_hom_alt，num_het，num_unknown：四种基因型的数量
aaf：所有sample等位基因的频率（即除开REF），返回列表。
heterozygosity：该位点的杂合度，0.5为杂合突变，0为纯合突变。
var_type：突变类型，包括‘snp’，‘indel’，‘sv’（structural variant），‘unknown’。
var_subtype：更加细化的突变类型，如‘indel’包括‘del’，‘ins’，‘unknown’。
is_snp，is_indel，is_sv，is_transition，is_deletion：判断突变是不是snp，indel，sv，transition，deletion等。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

record = next(vcf_reader)
print(record)
print(record.samples)
print(record.num_called)
print(record.call_rate)
print(record.num_hom_ref)
print(record.aaf)
print(record.num_het)
print(record.heterozygosity)
print(record.var_type)
print(record.var_subtype)
print(record.is_snp)
print(record.is_indel)

>>> 输出内容
[Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))]
1
1.0
1
[0.0]
0
0.0
unknown
unknown
False
False

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

[2] Reader对象-处理vcf文件，构建结构化信息

class Reader(fsock=None, filename=None, compressed=None, prepend_chr=False, strict_whitespace=False, encoding='ascii')

在读取vcf文件时，总共有六个参数可供选择。

fsock：目标文件的文件对象，可以用open(文件名)得到这个文件对象。
filename：文件名，当fsock和filename同时存在时，优先考虑fsock。
compressed：是否要解压，不提供参数时由程序自行判断（以文件名是否以.gz结尾判断是否需要解压）。
prepend_chr：在保存染色体名称时，是否加前缀‘chr’，默认不加，如果vcf文件的染色体名称本来没有前缀‘chr’，可设置为True，自动加上。
strict_whitespace：是否严格以制表符‘\t’分隔数据。True则表示严格按制表符分，False表示可以夹杂空格。
encoding：文件编码。

头文件信息主要保存在Reader对象的属性中，包括alts，contigs，filters，formats，infos，metadata。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
print(vcf_reader.alts)
print(vcf_reader.alts['NON_REF'].id)
print(vcf_reader.alts['NON_REF'].desc)

>>> 输出内容
OrderedDict([('NON_REF', Alt(id='NON_REF', desc='Represents any possible alternative allele at this location'))])
NON_REF
Represents any possible alternative allele at this location

1
2
3
4
5
6
7
8
9
10
11

# 3. 处理VCF文件并生成索引

# 3.1 工具概述

# 3.1.1 htslib工具

HTSlib 是一个统一的 C 库的实现，用于访问SAM、CRAM 和 VCF等常见文件格式，用于高通量测序数据，是samtools和bcftools使用的核心库。

项目地址：https://github.com/samtools/htslib (opens new window)
安装说明：https://github.com/samtools/htslib/blob/develop/INSTALL (opens new window)

# 3.1.2 bcftools工具

bcftools是一组用于变异获取和操作VCF和BCF的实用工具集合。参数众多，功能也十分强大。

项目地址：https://github.com/samtools/bcftools (opens new window)
官方文档：http://samtools.github.io/bcftools/bcftools.html (opens new window)

# 3.2 安装工具

htslib是bcftools使用的核心库，单独安装bcftools会缺失很多环境，建议先安装htslib，然后再安装bcftools。

以下在Debian 11 x86_64系统安装htslib是bcftools工具（不同系统的安装方式不同），基础编译环境如下：

$ sudo apt-get update  
$ sudo apt-get install autoconf
$ sudo apt-get install -y make

1
2
3

# 3.2.1 安装htslib

$ git clone http://github.com/samtools/htslib.git
$ cd htslib
$ sudo apt-get install autoconf automake make gcc perl zlib1g-dev libbz2-dev liblzma-dev libcurl4-gnutls-dev libssl-dev
$ autoreconf -i
$ git submodule update --init --recursive
$ ./configure    
$ make
$ make install

1
2
3
4
5
6
7
8

参考：https://github.com/samtools/htslib (opens new window)、https://github.com/samtools/htslib/blob/develop/INSTALL (opens new window)

# 3.2.2 安装bcftools

$ git clone http://github.com/samtools/bcftools.git
$ cd bcftools
$ autoheader && autoconf && ./configure --enable-libgsl --enable-perl-filters
$ make
$ make install

1
2
3
4
5

# 3.3 生成索引随机读取VCF位点

建立索引可以加快数据读取速度，而且建立索引后可以实现随机读取VCF位点。bcftools工具index命令用于对VCF文件建立索引，要求输入的VCF文件必须是使用bgzip压缩之后的文件，支持.csi和.tbi两种索引，默认情况下建立的索引是.csi格式。

[1] 将vcf文件压缩转换为vcf.gz

$ bcftools view your_file.vcf -Oz -o your_file.vcf.gz

[2] 为vcf.gz构建索引

$ bcftools index your_file.vcf.gz

注：必须要先压缩成gz格式，直接用vcf否则会报错index: the file is not BGZF compressed, cannot index: your_file.vcf。在我的实验环境下，123s对一个具有23亿（由于测序VCF里只有突变位置，所以实际远小于23亿）基因位点的生物样本建立了索引。

[3] 读取任意范围的位点信息

$ tabix your_file.vcf.gz chr_1:1-100000000

注：chr_1是染色体名，1-100000000是位点范围。随机读取的性能挺强的，在我的实验环境下，13秒取出了3600w基因位点信息。

# 4. 基因分析工具

# 4.1 Plink工具

Plink是一个免费的开源全基因组关联分析工具集，旨在以高效率计算的方式执行一系列基本的、大规模的分析。Plink的重点分析对象是基因型或者表型数据，可为后续的可视化、注释和结果存储提供一些支持。

Plink官网：http://zzz.bwh.harvard.edu/plink/ (opens new window)
Plink功能：数据管理、质控、人群分层检测、关联分析-SNP、多marker预测、haplotype分析、拷贝数变异分析、异位显性分析、关联分析-基因、基因环境互作分析、meta分析等。

Plink工具

# 4.2 Admixture工具

Admixture是一款快速分析群体遗传结构的工具，一般进行群体分析所使用的软件为Structure，但Structure的运行速度较慢，如今Admixture凭借其高速的运算速度逐渐成为群体遗传结构分析的主流软件。

Admixture官网：https://bioinformaticshome.com/tools/descriptions/ADMIXTURE.html (opens new window)
群体遗传结构（Structure）指遗传变异在物种或群体中的一种非随机分布。按照地理分布或其他标准可将一个群体分为若干亚群，处于同一亚群内的不同个体亲缘关系较高，而亚群与亚群之间则亲缘关系稍远。群体结构分析有助于理解进化过程，并且可以通过基因型和表型的关联研究确定个体所属的亚群。

Admixture工具