pyvcf

1. 读取vcf 与写出

-读取

import vcf
vcf_reader = vcf.Reader(open('vcf/test/example-4.0.vcf', 'r'))
for record in vcf_reader:
  print (record)

• 写出

>>> vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
>>> for record in vcf_reader:
... vcf_writer.write_record(record)

2. recode 属性

• record 里面有8个子属性，分别对应vcf head 行

Record.CHROM 
Record.POS
Record.ID
Record.REF
Record.ALT
Record.QUAL
Record.FILTER
Record.INFO # 返回一个字典，可以用Record.INFO['type'],Record.INFO['DP'] 键值继续提取
Record.FORMAT  #返回format列 字符串  如果你的vcf文件中没有FORMAT 返回 "GT:DP:RO:QR:AO:QA:GL"
print i.samples  # 返回的是三个样 call object 组成的列表。

# 使用方法
print(record.POS)

3. 提取基因型

for record in vcf_reader:
        for j in record.samples:
            print(j['GT']) 

4. fetch

• 如果你不想从头到尾循环某个文件。只想取某一部分，可以用fetch, 但是前体是用tabix对文件index， tabix前要用bgzip压缩
• bgzip testpyvcf.vcf #得到testpyvcf.vcf.gz文件
• tabix -p vcf testpyvcf.vcf.gz #得到testpyvcf.vcf.gz的index文件：testpyvcf.vcf.gz.tbi

 import vcf 
 myvcf = vcf.Reader(filename='testpyvcf.vcf.gz')
 for i in myvcf.fetch('Chr1', 1111, 444444):      #第一个是序列名，第二个起始，第三个end,包括。

1. https://buildmedia.readthedocs.org/media/pdf/pyvcf/latest/pyvcf.pdf

看了这个感觉还是没解决我的问题。