@@ -119,6 +119,122 @@ formatted file on stdout::
119119
120120Note that the file open mode needs to changed from ``r `` to ``rb ``.
121121
122+
123+ ====================================
124+ Working with VCF/BCF-formatted files
125+ ====================================
126+
127+ To iterate through a VCF/BCF formatted file use
128+ :class: `~pysam.VariantFile `::
129+
130+ from pysam import VariantFile
131+
132+ bcf_in = VariantFile("test.bcf") # auto-detect input format
133+ bcf_out = VariantFile('-', 'w', header=bcf_in.header)
134+
135+ for rec in bcf_in.fetch('chr1', 100000, 200000):
136+ bcf_out.write(rec)
137+
138+ :meth: `_pysam.VariantFile.fetch() ` iterates over
139+ :class: `~pysam.VariantRecord ` objects which provides access to
140+ simple variant attributes such as :class: `~pysam.VariantRecord.contig `,
141+ :class: `~pysam.VariantRecord.pos `, :class: `~pysam.VariantRecord.ref `::
142+
143+ for rec in bcf_in.fetch():
144+ print(rec.pos)
145+
146+ but also to complex attributes such as the contents to the
147+ :class: `~pysam.VariantRecord.info `, :class: `~pysam.VariantRecord.format `
148+ and :term: `genotype ` columns. These
149+ complex attributes are views on the underlying htslib data structures
150+ and provide dictionary-like access to the data::
151+
152+ for rec in bcf_in.fetch():
153+ print(rec.info)
154+ print(rec.info.keys())
155+ print(rec.info["DP"])
156+
157+ The :py:attr: `~pysam.VariantFile.header ` attribute
158+ (:class: `~pysam.VariantHeader `) provides access information
159+ stored in the :term: `vcf ` header. The complete header can be printed::
160+
161+ >>> print(bcf_in.header)
162+ ##fileformat=VCFv4.2
163+ ##FILTER=<ID=PASS,Description="All filters passed">
164+ ##fileDate=20090805
165+ ##source=myImputationProgramV3.1
166+ ##reference=1000GenomesPilot-NCBI36
167+ ##phasing=partial
168+ ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples
169+ With Data">
170+ ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
171+ ##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
172+ ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
173+ ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build
174+ 129">
175+ ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
176+ ##FILTER=<ID=q10,Description="Quality below 10">
177+ ##FILTER=<ID=s50,Description="Less than 50% of samples have data">
178+ ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
179+ ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
180+ ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
181+ ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
182+ ##contig=<ID=M>
183+ ##contig=<ID=17>
184+ ##contig=<ID=20>
185+ ##bcftools_viewVersion=1.3+htslib-1.3
186+ ##bcftools_viewCommand=view -O b -o example_vcf42.bcf
187+ example_vcf42.vcf.gz
188+ #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA0000
189+
190+ Individual contents such as contigs, info fields, samples, formats can
191+ be retrieved as attributes from :py:attr: `~pysam.VariantFile.header `::
192+
193+ >>> print(bcf_in.header.contigs)
194+ <pysam.cbcf.VariantHeaderContigs object at 0xf250f8>
195+
196+ To convert these views to native python types, iterate through the views::
197+
198+ >>> print(list((bcf_in.header.contigs)))
199+ ['M', '17', '20']
200+ >>> print(list((bcf_in.header.filters)))
201+ ['PASS', 'q10', 's50']
202+ >>> print(list((bcf_in.header.info)))
203+ ['NS', 'DP', 'AF', 'AA', 'DB', 'H2']
204+ >>> print(list((bcf_in.header.samples)))
205+ ['NA00001', 'NA00002', 'NA00003']
206+
207+ Alternatively, it is possible to iterate through all records in the
208+ header returning objects of type :py:class: `~pysam.VariantHeaderRecord `:: ::
209+
210+ >>> for x in bcf_in.header.records:
211+ >>> print(x)
212+ >>> print(x.type, x.key)
213+ GENERIC fileformat
214+ FILTER FILTER
215+ GENERIC fileDate
216+ GENERIC source
217+ GENERIC reference
218+ GENERIC phasing
219+ INFO INFO
220+ INFO INFO
221+ INFO INFO
222+ INFO INFO
223+ INFO INFO
224+ INFO INFO
225+ FILTER FILTER
226+ FILTER FILTER
227+ FORMAT FORMAT
228+ FORMAT FORMAT
229+ FORMAT FORMAT
230+ FORMAT FORMAT
231+ CONTIG contig
232+ CONTIG contig
233+ CONTIG contig
234+ GENERIC bcftools_viewVersion
235+ GENERIC bcftools_viewCommand
236+
237+
122238==================================================
123239Using samtools and bcftools commands within Python
124240==================================================
@@ -294,120 +410,6 @@ become available through named access, for example::
294410.. and :class:`pysam.IteratorIndelCaller`.
295411
296412
297- ====================================
298- Working with VCF/BCF formatted files
299- ====================================
300-
301- To iterate through a VCF/BCF formatted file use
302- :class: `~pysam.VariantFile `::
303-
304- from pysam import VariantFile
305-
306- bcf_in = VariantFile("test.bcf") # auto-detect input format
307- bcf_out = VariantFile('-', 'w', header=bcf_in.header)
308-
309- for rec in bcf_in.fetch('chr1', 100000, 200000):
310- bcf_out.write(rec)
311-
312- :meth: `_pysam.VariantFile.fetch() ` iterates over
313- :class: `~pysam.VariantRecord ` objects which provides access to
314- simple variant attributes such as :class: `~pysam.VariantRecord.contig `,
315- :class: `~pysam.VariantRecord.pos `, :class: `~pysam.VariantRecord.ref `::
316-
317- for rec in bcf_in.fetch():
318- print(rec.pos)
319-
320- but also to complex attributes such as the contents to the
321- :class: `~pysam.VariantRecord.info `, :class: `~pysam.VariantRecord.format `
322- and :term: `genotype ` columns. These
323- complex attributes are views on the underlying htslib data structures
324- and provide dictionary-like access to the data::
325-
326- for rec in bcf_in.fetch():
327- print(rec.info)
328- print(rec.info.keys())
329- print(rec.info["DP"])
330-
331- The :py:attr: `~pysam.VariantFile.header ` attribute
332- (:class: `~pysam.VariantHeader `) provides access information
333- stored in the :term: `vcf ` header. The complete header can be printed::
334-
335- >>> print(bcf_in.header)
336- ##fileformat=VCFv4.2
337- ##FILTER=<ID=PASS,Description="All filters passed">
338- ##fileDate=20090805
339- ##source=myImputationProgramV3.1
340- ##reference=1000GenomesPilot-NCBI36
341- ##phasing=partial
342- ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples
343- With Data">
344- ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
345- ##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
346- ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
347- ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build
348- 129">
349- ##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
350- ##FILTER=<ID=q10,Description="Quality below 10">
351- ##FILTER=<ID=s50,Description="Less than 50% of samples have data">
352- ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
353- ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
354- ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
355- ##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
356- ##contig=<ID=M>
357- ##contig=<ID=17>
358- ##contig=<ID=20>
359- ##bcftools_viewVersion=1.3+htslib-1.3
360- ##bcftools_viewCommand=view -O b -o example_vcf42.bcf
361- example_vcf42.vcf.gz
362- #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA0000
363-
364- Individual contents such as contigs, info fields, samples, formats can
365- be retrieved as attributes from :py:attr: `~pysam.VariantFile.header `::
366-
367- >>> print(bcf_in.header.contigs)
368- <pysam.cbcf.VariantHeaderContigs object at 0xf250f8>
369-
370- To convert these views to native python types, iterate through the views::
371-
372- >>> print(list((bcf_in.header.contigs)))
373- ['M', '17', '20']
374- >>> print(list((bcf_in.header.filters)))
375- ['PASS', 'q10', 's50']
376- >>> print(list((bcf_in.header.info)))
377- ['NS', 'DP', 'AF', 'AA', 'DB', 'H2']
378- >>> print(list((bcf_in.header.samples)))
379- ['NA00001', 'NA00002', 'NA00003']
380-
381- Alternatively, it is possible to iterate through all records in the
382- header returning objects of type :py:class: `~pysam.VariantHeaderRecord `:: ::
383-
384- >>> for x in bcf_in.header.records:
385- >>> print(x)
386- >>> print(x.type, x.key)
387- GENERIC fileformat
388- FILTER FILTER
389- GENERIC fileDate
390- GENERIC source
391- GENERIC reference
392- GENERIC phasing
393- INFO INFO
394- INFO INFO
395- INFO INFO
396- INFO INFO
397- INFO INFO
398- INFO INFO
399- FILTER FILTER
400- FILTER FILTER
401- FORMAT FORMAT
402- FORMAT FORMAT
403- FORMAT FORMAT
404- FORMAT FORMAT
405- CONTIG contig
406- CONTIG contig
407- CONTIG contig
408- GENERIC bcftools_viewVersion
409- GENERIC bcftools_viewCommand
410-
411413===============
412414Extending pysam
413415===============
0 commit comments