FAQ

How should I cite pysam

Pysam has not been published in print. When referring to pysam, please use the github URL: https://github.com/pysam-developers/pysam. As pysam is a wrapper around htslib and the samtools package, I suggest citing [Li.2009], [Bonfield.2021], and/or [Danecek.2021], as appropriate.

Is pysam thread-safe?

Pysam is a mix of python and C code. Instructions within python are generally made thread-safe through python’s global interpreter lock (GIL). This ensures that python data structures will always be in a consistent state.

If an external function outside python is called, the programmer has a choice to keep the GIL in place or to release it. Keeping the GIL in place will make sure that all python threads wait until the external function has completed. This is a safe option and ensures thread-safety.

Alternatively, the GIL can be released while the external function is called. This will allow other threads to run concurrently. This can be beneficial if the external function is expected to halt, for example when waiting for data to read or write. However, to achieve thread-safety, the external function needs to be implemented with thread-safety in mind. This means that there can be no shared state between threads, or if there is shared, it needs to be controlled to prevent any access conflicts.

Pysam generally uses the latter option and aims to release the GIL for I/O intensive tasks. This is generally fine, but thread-safety of all parts have not been fully tested.

A related issue is when different threads read from the same file object - or the same thread uses two iterators over a file. There is only a single file-position for each opened file. To prevent this from happening, use the option multiple_iterators=True when calling a fetch() method. This will return an iterator on a newly opened file.

pysam coordinates are wrong

pysam uses 0-based coordinates and the half-open notation for ranges as does python. Coordinates and intervals reported from pysam always follow that convention.

Confusion might arise as different file formats might have different conventions. For example, the SAM format is 1-based while the BAM format is 0-based. It is important to remember that pysam will always conform to the python convention and translate to/from the file format automatically.

The only exception is the region string in the fetch() and pileup() methods. This string follows the convention of the samtools command line utilities. The same is true for any coordinates passed to the samtools command utilities directly, such as pysam.mpileup().

Calling pysam.fetch() confuses existing iterators

The following code will cause unexpected behaviour:

samfile = pysam.AlignmentFile("pysam_ex1.bam", "rb")

iter1 = samfile.fetch("chr1")
print(next(iter1).reference_id)
iter2 = samfile.fetch("chr2")
print(next(iter2).reference_id)
print(next(iter1).reference_id)

This will give the following output:

0
1
Traceback (most recent call last):
  File "xx.py", line 8, in <module>
    print(next(iter1).reference_id)
  File "libcalignmentfile.pyx", line 2103,
  in pysam.libcalignmentfile.IteratorRowRegion.__next__
StopIteration

Note how the second iterator stops as the file pointer has moved to chr2. The correct way to work with multiple iterators is:

samfile = pysam.AlignmentFile("pysam_ex1.bam", "rb")

iter1 = samfile.fetch("chr1", multiple_iterators=True)
print(next(iter1).reference_id)
iter2 = samfile.fetch("chr2")
print(next(iter2).reference_id)
print(next(iter1).reference_id)

Here, the output is:

0
1
0

The reason for this behaviour is that every iterator needs to keep track of its current position in the file. Within pysam, each opened file can only keep track of one file position and hence there can only be one iterator per file. Using the option multiple_iterators=True will return an iterator within a newly opened file. This iterator will not interfere with existing iterators as it has its own file handle associated with it.

Note that re-opening files incurs a performance penalty which can become severe when calling fetch() often. Thus, multiple_iterators is set to False by default.

AlignmentFile.fetch does not show unmapped reads

By default, fetch() will only iterate over placed alignments in the SAM/BAM/CRAM file. Thus the following always works:

f = pysam.AlignmentFile(fname, "r")
for r in f.fetch():
    assert r.reference_name is not None

If the file contains unaligned reads, they can be included in the iteration by adding the until_eof=True flag:

f = pysam.AlignmentFile(fname, "r")
for r in f.fetch(until_eof=True):
    if r.is_unmapped:
        print("read is unmapped")

See also fetch("*") which iterates only over the unplaced unmapped reads at the end of the file.

I can’t call AlignmentFile.fetch on a file without an index

fetch() requires an index when iterating over a SAM/BAM file. To iterate over a file without an index, use until_eof=True:

bf = pysam.AlignmentFile(fname, "rb")
for r in bf.fetch(until_eof=True):
    print(r)

BAM files with a large number of reference sequences are slow

If you have many reference sequences in a BAM file, the following might be slow:

track = pysam.AlignmentFile(fname, "rb")
for aln in track.fetch():
    pass

The reason is that track.fetch() will iterate through the BAM file for each reference sequence in the order as it is defined in the header. This might require a lot of jumping around in the file. To avoid this, use:

track = pysam.AlignmentFile(fname, "rb")
for aln in track.fetch(until_eof=True):
    pass

This will iterate through reads as they appear in the file.

Weirdness with spliced reads in samfile.pileup(chr,start,end) given spliced alignments from an RNA-seq bam file

Spliced reads are reported within samfile.pileup. To ignore these in your analysis, test the flags is_del == True and indel == 0 in the PileupRead object.

I can’t edit quality scores in place

Editing reads in-place generally works, though there is one quirk to be aware of. Assigning to AlignedSegment.query_sequence will invalidate any quality scores in AlignedSegment.query_qualities. The reason is that samtools manages the memory of the sequence and quality scores together and thus requires them to always be of the same length or 0.

Thus, to in-place edit the sequence and quality scores, copies of the quality scores need to be taken. Consider trimming for example:

quals = read.query_qualities
read.query_sequence = read.query_sequence[5:10]
read.query_qualities = quals[5:10]

Why is there no SNPCaller class anymore?

SNP calling is highly complex and heavily parameterized. There was a danger that the pysam implementations might show different behaviour from the samtools implementation, which would have caused a lot of confusion.

The best way to use samtools SNP calling from python is to use the pysam.mpileup() command and parse the output directly.

I get an error ‘PileupProxy accessed after iterator finished’

Pysam works by providing proxy objects to objects defined within the C-samtools package. Thus, some attention must be paid to the lifetime of objects. The following to code snippets will cause an error:

s = AlignmentFile('ex1.bam')
for p in s.pileup('chr1', 1000,1010):
    pass

for pp in p.pileups:
    print(pp)

The iteration has finished, thus the contents of p are invalid. Another variation of this:

p = next(AlignmentFile('ex1.bam').pileup('chr1', 1000, 1010))
for pp in p.pileups:
    print(pp)

Again, the iteration finishes as the temporary iterator created by pileup goes out of scope. The solution is to keep a handle to the iterator that remains alive:

i = AlignmentFile('ex1.bam').pileup('chr1', 1000, 1010)
p = next(i)
for pp in p.pileups:
    print(pp)

Pysam won’t compile

Compiling pysam can be tricky as there are numerous variables that differ between build environments such as OS, version, python version, and compiler. It is difficult to build software that builds cleanly on all systems and the process might fail. Please see the pysam user group for common issues.

If there is a build issue, read the generated output carefully - generally the cause of the problem is among the first errors to be reported. For example, you will need to have the development version of python installed that includes the header files such as Python.h. If that file is missing, the compiler will report this at the very top of its error messages but will follow it with any unknown function or variable definition it encounters later on.

General advice is to always use the latest version on python and cython when building pysam. There are some known incompatibilities:

Python 3.4 requires cython 0.20.2 or later (see here)

ImportError: cannot import name csamtools

In version 0.10.0 and onwards, all pysam extension modules contain a lib-prefix. This facilates linking against pysam extension modules with compilers that require to start with lib. As a consequence, all code using pysam extension modules directly will need to be adapted. For example,:

cimport pysam.csamtools

will become:

cimport pysam.libcsamtools