Python逐行读取VCF文件

VCF（Variant Call Format）是一种常用的生物数据文件格式，用于存储单个或多个样本的基因变异信息。本文介绍如何使用Python逐行读取VCF文件。

准备工作

在开始之前，请确保已经安装了Python，并且安装了相应的依赖库。本文例子中使用的依赖库有pandas和numpy。

首先，我们需要下载一个VCF文件来进行测试。可以在1000 Genomes Project网站上找到一些免费的VCF文件。这里选择ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz文件。

然后，我们需要使用gunzip命令解压文件：

gunzip ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

解压后，我们用文本编辑器打开这个文件可以看到如下内容：

##fileformat=VCFv4.1
##fileDate=20130502
##source=1000GenomesPhase3Pipeline
##reference=human_g1k_v37.fasta.gz
##contig=<ID=22,length=51304566,assembly=b37>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NA00001 NA00002 NA00003 NA00004
22      16053794        rs138196087     C       T       .       PASS    AC=1;AF=0.000199681;AN=5008;DP=58;NS=2504;EAS_AF=0;AMR_AF=0.001;AFR_AF=0;EUR_AF=0;SAS_AF=0;AA=T|||;VT=SNP   GT:AD:DP:GQ:PL  0|0:26,0:26:78:0,78,936 0|0:22,0:22:66:0,66,792 0|0:25,0:25:75:0,75,900 0|0:30,0:30:90:0,90,1080
22      16053987        rs181691356     A       G       .       PASS    AC=1;AF=0.000199681;AN=5008;DP=57;NS=2504;EAS_AF=0;AMR_AF=0;AFR_AF=0;EUR_AF=0;SAS_AF=0;AA=G|||;VT=SNP   GT:AD:DP:GQ:PL  0|0:24,0:24:72:0,72,864 0|0:22,0:22:63:0,63,756 0|0:29,0:29:87:0,87,1044 0|0:29,0:29:87:0,87,1044
22      16054053        rs201431404     A       G       .       PASS    AC=1;AF=0.000199681;AN=5008;DP=63;NS=2504;EAS_AF=0;AMR_AF=0;AFR_AF=0;EUR_AF=0;SAS_AF=0;AA=G|||;VT=SNP   GT:AD:DP:GQ:PL  0|0:25,0:25:75:0,75,900 0|0:25,0:25:75:0,75,900 0|0:24,0:24:71:0,71,852 0|0:29,0:29:87:0,87,1044
...

可以看到，这个文件包括了注释行、格式行和数据行等部分。每一行通过\n进行换行。

逐行读取VCF文件

Python中提供了open()函数来打开文件并且读取文件内容。我们可以使用with语句来自动关闭文件。

with open('ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf', 'r') as f:
    for line in f:
        print(line)

这个程序会一行一行地读取文件内容，然后将每一行打印出来。

过滤注释行和格式行

在VCF文件中，注释行以##开头，格式行以#开头。我们可以根据行的开头来过滤掉这些行。

with open('ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf', 'r') as f:
    for line in f:
        if line.startswith('##'):
            continue
        if line.startswith('#'):
            header = line.strip().split('\t') # 获取列名
            continue
        print(line)

这个程序会先过滤掉注释行，然后获取格式行的列名，最后打印出数据行。

解析数据行

在VCF文件中，数据行有许多列，包括染色体位置、ID、引用序列、变异序列、质量等级、过滤结果、信息等级以及基因型等等。我们可以使用pandas库来解析这些数据。

import pandas as pd

with open('ALL.chr22.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf', 'r') as f:
    rows = []
    for line in f:
        if line.startswith('##'):
            continue
        if line.startswith('#'):
            header = line.strip().split('\t')
            continue
        fields = line.strip().split('\t')
        row = dict(zip(header, fields))
        rows.append(row)

df = pd.DataFrame(rows)

print(df.head())

这个程序会解析出每一行的数据，并且将数据存储到一个pandas的DataFrame中。最后打印出DataFrame的前5行。

总结

本文介绍了如何使用Python逐行读取VCF文件，并且解析其中的注释行、格式行和数据行。此外，还使用pandas库解析了数据行，并且将其存储到一个DataFrame中。