mgkit.align module

Module dealing with BAM/SAM files

class mgkit.align.SamtoolsDepth(file_handle, num_seqs=10000, max_size=1000000, max_size_dict=None, calc_density=False, dtype='uint32')

Bases: object

Changed in version 0.4.2: several optimisations and changes to support a scanning approach, instead of lookup table. No exception is raised when a sequence is not found in the file, instead assuming that the coverage is 0

Changed in version 0.4.0: uses pandas.SparseArray now. It should use less memory, but needs pandas version > 0.24

New in version 0.3.0.

This a class that helps use the results from read_samtools_depth().

There are 2 modes of operations:

1. Request a region coverage via SamtoolsDepth.region_coverage()

2. Advance the samtools depth file until a sequence coverage is read

The method 1) was the default in MGKit < 0.4.2, when the user would request a region coverage and this class would advance the reading until the requested region is found. While scanning, all the sequences encountered are kept as pandas Series with SparseArray declared, reindexed from 0 to the max_size_dict values if provided, or max_size. The advantage is that the it’s easier to use, but with bigger datasets will keep high memory usage. This is the case for case 2) which involves calling directly SamtoolsDepth.advance_file() returning the sequence ID read. Then the region coverage can be requested the same way as before, but won’t involve scanning the file.

The use of a SparseArray in a pandas Series allows the use of samtools depth files that weren’t produced with samtools depth -aa and save memory.

Note

The amount of memory used is fairly high, especially with the increasing number of sequences in a GFF/Depth file. It is recoommended to use max_size_dict to 1) only creates arrays for sequences needed and 2) use arrays of smaller size

Warning

In MGKit 0.4.2, the internal dictionary to keep the SparseArray(s) is a weakref.WeakValueDictionary, to improve the release of memory. That on Linux seems to release the array too fast. Upgrade to MGKit version 0.4.3 if this class is needed

advance_file()

New in version 0.4.2.

Requests the scan of the file and returns the sequence found, the coverage can then be retrivied with SamtoolsDepth.region_coverage()

Returns

the return value is None if the file is fully read, otherwise the the sequnece ID found is returned

Return type

(str, None)

closed = False

data = None

density = None

drop_sequence(seq_id)

New in version 0.4.2.

Remove the sequence passed from the internal dictionary

file_handle = None

max_size = None

max_size_dict = None

not_found = None

region_coverage(seq_id, start, end)

Changed in version 0.4.2: now using SamtoolsDepth.advance_file() to scan the file

Returns the mean coverage of a region. The start and end parameters are expected to be 1-based coordinates, like the correspondent attributes in mgkit.io.gff.Annotation or mgkit.io.gff.GenomicRange.

If the sequence for which the coverage is requested is not found, the depth file is read (and cached) until it is found.

Parameters

• seq_id (str) – sequence for which to return mean coverage

• start (int) – start of the region

• end (int) – end of the region

Returns

mean coverage of the requested region

Return type

float

mgkit.align.add_coverage_info(annotations, bam_files, samples, attr_suff='_cov')

Changed in version 0.3.4: the coverage now is returned as floats instead of int

Adds coverage information to annotations, using BAM files.

The coverage information is added for each sample as a ‘sample_cov’ and the total coverage as as ‘cov’ attribute in the annotations.

Note

The bam_files and sample variables must have the same order

Parameters

• annotations (iterable) – iterable of annotations

• bam_files (iterable) – iterable of pysam.Samfile instances

• sample (iterable) – names of the samples for the BAM files

mgkit.align.covered_annotation_bp(files, annotations, min_cov=1, progress=False)

New in version 0.1.14.

Returns the number of base pairs covered of annotations over multiple samples.

Parameters

• files (iterable) – an iterable that returns the alignment file names

• annotations (iterable) – an iterable that returns annotations

• min_cov (int) – minumum coverage for a base to counted

• progress (bool) – if True, a progress bar is used

Returns

a dictionary whose keys are the uid and the values the number of bases that are covered by reads among all samples

Return type

dict

mgkit.align.get_region_coverage(bam_file, seq_id, feat_from, feat_to)

Return coverage for an annotation.

Note

feat_from and feat_to are 1-based indexes

Parameters

• bam_file (Samfile) – instance of pysam.Samfile

• seq_id (str) – sequence id

• feat_from (int) – start position of feature

• feat_to (int) – end position of feature

Return int

coverage array for the annotation

mgkit.align.read_samtools_depth(file_handle, num_seqs=10000, seq_ids=None)

Changed in version 0.4.2: the function returns lists instead of numpy arrays for speed (at least in my tests it seems ~4x increase)

Changed in version 0.4.0: now returns 3 array, instead of 2. Also added seq_ids to skip lines

Changed in version 0.3.4: num_seqs can be None to avoid a log message

New in version 0.3.0.

Reads a samtools depth file, returning a generator that yields the array of each base coverage on a per-sequence base.

Note

There’s no need anymore to use samtools depth -aa, because the function returns the position array and this can be used to create a Pandas SparseArray which can be reindexed to include missing positions (with values of 0)

Valid for version < 0.4.0:

The information on position is not used, to use numpy and save memory. samtools depth should be called with the -aa option:

`samtools depth -aa bamfile`

This options will output both base position with 0 coverage and sequneces with no aligned reads

Parameters

• file_handle (file) – file handle of the coverage file

• num_seqs (int or None) – number of sequence that fires a log message. If None, no message is triggered

• seq_ids (dict, set) – a hashed container like a dictionary or set with the sequences to return

Yields

tuple – the first element is the sequence identifier, the second one is the list with the positions, the third element is the list with the coverages