seismicrna.idmut package

seismicrna.idmut.run(fasta: str | Path = Sentinel.UNSET, input_path: Iterable[str | Path] = Sentinel.UNSET, *, out_dir: str | Path = './out', branch: str = '', min_reads: int = 1000, min_mapq: int = 25, phred_enc: int = 33, min_phred: int = 25, batch_size: int = 65536, insert3: bool = True, ambindel: bool = True, overhangs: bool = True, clip_end5: int = 4, clip_end3: int = 4, sep_strands: bool = False, rev_label: str = '-rev', write_read_names: bool = False, idmut_pos_table: bool = True, idmut_read_table: bool = False, idmut_cx: bool = True, num_cpus: int = 4, brotli_level: int = 10, force: bool = False, keep_tmp: bool = False, tmp_pfx='./tmp')

Compute relationships between references and aligned reads.

Parameters:

  • out_dir (str | pathlib._local.Path) – Write all output files to this directory [keyword-only, default: ‘./out’]

  • branch (str) – Create a new branch of the workflow with this name [keyword-only, default: ‘’]

  • min_reads (int) – Discard alignment maps with fewer than this many reads [keyword-only, default: 1000]

  • min_mapq (int) – Discard reads with mapping qualities below this threshold [keyword-only, default: 25]

  • phred_enc (int) – Specify the Phred score encoding of FASTQ and SAM/BAM/CRAM files [keyword-only, default: 33]

  • min_phred (int) – Mark base calls with Phred scores lower than this threshold as ambiguous [keyword-only, default: 25]

  • batch_size (int) – Limit batches to at most this many reads [keyword-only, default: 65536]

  • insert3 (bool) – Mark each insertion on the base to its 3’ (True) or 5’ (False) side [keyword-only, default: True]

  • ambindel (bool) – Mark all ambiguous insertions and deletions (indels) [keyword-only, default: True]

  • overhangs (bool) – Retain the overhangs of paired-end mates that dovetail [keyword-only, default: True]

  • clip_end5 (int) – Clip this many bases from the 5’ end of each read [keyword-only, default: 4]

  • clip_end3 (int) – Clip this many bases from the 3’ end of each read [keyword-only, default: 4]

  • sep_strands (bool) – Separate each alignment map into forward- and reverse-strand reads [keyword-only, default: False]

  • rev_label (str) – With –sep-strands, add this label to each reverse-strand reference [keyword-only, default: ‘-rev’]

  • write_read_names (bool) – Write the name of each read in a second set of batches (necessary for the options –drop-read or –drop-read-file) [keyword-only, default: False]

  • idmut_pos_table (bool) – Tabulate relationships per position for idmut data [keyword-only, default: True]

  • idmut_read_table (bool) – Tabulate relationships per read for idmut data [keyword-only, default: False]

  • idmut_cx (bool) – Use a fast (C extension module) version of the idmut algorithm; the slow (Python) version is still avilable as a fallback if the C extension cannot be loaded, and for debugging/benchmarking [keyword-only, default: True]

  • num_cpus (int) – Use up to this many CPUs simultaneously [keyword-only, default: 4]

  • brotli_level (int) – Compress pickle files with this level of Brotli (0 - 11) [keyword-only, default: 10]

  • force (bool) – Force all tasks to run, overwriting any existing output files [keyword-only, default: False]

  • keep_tmp (bool) – Keep temporary files after finishing [keyword-only, default: False]

  • tmp_pfx – Write all temporary files to a directory with this prefix [keyword-only, default: ‘./tmp’]

Subpackages

Submodules

class seismicrna.idmut.batch.FullReadBatch(*, batch: int, **kwargs)

Bases: ReadBatch, ABC

property max_read

Maximum possible value for a read index.

property read_indexes

Map each read number to its index in self.read_nums.

property read_nums

Read numbers.

abstractmethod classmethod simulate(*args, **kwargs) → Self

Simulate a batch.

class seismicrna.idmut.batch.IDmutMutsBatch(*, region: Region, sanitize: bool = True, muts: dict[int, dict[int, list[int] | np.ndarray]], masked_read_nums: np.ndarray | list[int] | None = None, **kwargs)

Bases: FullReadBatch, MutsBatch, ABC

property read_weights

Weights for each read when computing counts.

class seismicrna.idmut.batch.IDmutRegionMutsBatch(*, region: Region, **kwargs)

Bases: IDmutMutsBatch, RegionMutsBatch

classmethod simulate(ref: str, pmut: pd.DataFrame, uniq_end5s: np.ndarray, uniq_end3s: np.ndarray, pends: np.ndarray, paired: bool, read_length: int, p_rev: float, min_mut_gap: int, injected_mut_probs: dict[int, float], mut_collisions: str, num_reads: int, seed: int | None, **kwargs)

Simulate a batch.

Parameters:

  • ref (str) – Name of the reference.

  • pmut (pd.DataFrame) – Rate of each type of mutation at each position.

  • uniq_end5s (np.ndarray) – Unique read 5’ end coordinates.

  • uniq_end3s (np.ndarray) – Unique read 3’ end coordinates.

  • pends (np.ndarray) – Probability of each set of unique end coordinates.

  • paired (bool) – Whether to simulate paired-end or single-end reads.

  • read_length (int) – Length of each read segment (paired-end reads only).

  • p_rev (float) – Probability that mate 1 is reversed (paired-end reads only).

  • min_mut_gap (int) – Minimum number of positions between two mutations.

  • injected_mut_probs (dict[int, float]) – Mapping of offset (positions 5’ of an existing mutation) to injection probability; passed directly to inject_close_muts. An empty dict disables injection.

  • mut_collisions (str) – How to handle reads with mutations closer than min_mut_gap: “drop” to remove such reads, or “merge” to merge them.

  • num_reads (int) – Number of reads in the batch.

  • seed (int | None) – Random seed for reproducibility; None for no fixed seed.

class seismicrna.idmut.batch.ReadNamesBatch(*, names: list[str] | np.ndarray, **kwargs)

Bases: FullReadBatch

property num_reads

Number of reads.

classmethod simulate(branches: dict[str, str], batch: int, num_reads: int, formatter: ~typing.Callable[[int, int], str] = <function format_read_name>, **kwargs)

Simulate a batch.

Parameters:

  • branches (dict[str, str]) – Branches of the workflow.

  • batch (int) – Batch number.

  • num_reads (int) – Number of reads in the batch.

  • formatter (Callable[[int, int], str]) – Function to generate the name of each read: must accept the batch number and the read number and return a string.

seismicrna.idmut.batch.format_read_name(batch_num: int, read_num: int)

Format a read name.

class seismicrna.idmut.dataset.AverageDataset(report_file: str | Path, verify_times: bool = True)

Bases: Dataset, ABC

Dataset of population average data.

property best_k

Best number of clusters.

property ks

Numbers of clusters.

class seismicrna.idmut.dataset.IDmutDataset(report_file: str | Path, verify_times: bool = True)

Bases: AverageDataset, ABC

Dataset of relationships.

class seismicrna.idmut.dataset.IDmutMutsDataset(report_file: str | Path, verify_times: bool = True)

Bases: IDmutDataset, LoadedDataset, MutsDataset

Dataset of mutations from the IDmut step.

get_batch(batch: int)

Get a specific batch of data.

classmethod get_batch_type()

Type of batch.

classmethod get_report_type()

Type of report.

property paired

Whether the reads are paired-end.

property pattern

Pattern of mutations to count.

property refseq

Sequence of the reference.

property region

Region of the dataset.

class seismicrna.idmut.dataset.NamesDataset(report_file: str | Path, verify_times: bool = True)

Bases: AverageDataset, ABC

classmethod kind()
class seismicrna.idmut.dataset.PoolDataset(*args, **kwargs)

Bases: TallDataset, ABC

Pooled dataset of relationships.

class seismicrna.idmut.dataset.PoolMutsDataset(*args, **kwargs)

Bases: IDmutDataset, PoolDataset, MutsDataset, MergedRegionDataset

Load pooled batches of relationships.

classmethod get_dataset_load_func()

Function to load one constituent dataset.

classmethod get_report_type()

Type of report.

property region

Region of the dataset.

class seismicrna.idmut.dataset.PoolReadNamesDataset(*args, **kwargs)

Bases: NamesDataset, PoolDataset

Pooled Dataset of read names.

classmethod get_dataset_load_func()

Function to load one constituent dataset.

classmethod get_report_type()

Type of report.

class seismicrna.idmut.dataset.ReadNamesDataset(report_file: str | Path, verify_times: bool = True)

Bases: NamesDataset, LoadedDataset

Dataset of read names from the IDmut step.

classmethod get_batch_type()

Type of batch.

classmethod get_report_type()

Type of report.

property pattern

Pattern of mutations to count.

class seismicrna.idmut.io.IDmutBatchIO(*args, region: Region, **kwargs)

Bases: IDmutMutsBatch, MutsBatchIO, RefBrickleIO, IDmutIO

classmethod from_region_batch(batch: IDmutRegionMutsBatch, *, sample: str, branches: dict[str, str])

Create an instance from an IDmutRegionMutsBatch.

classmethod get_file_seg_type()

Type of the last segment in the path.

classmethod simulate(*args, sample: str, branches: dict[str, str], **kwargs)

Simulate a batch.

to_region_batch(region: Region)

Create an IDmutRegionMutsBatch from this instance.

class seismicrna.idmut.io.IDmutFile

Bases: HasRefFilePath, ABC

classmethod get_step()

Step of the workflow.

class seismicrna.idmut.io.IDmutIO

Bases: IDmutFile, RefFileIO, ABC

class seismicrna.idmut.io.ReadNamesBatchIO(*, names: list[str] | np.ndarray, **kwargs)

Bases: ReadNamesBatch, ReadBatchIO, RefBrickleIO, IDmutIO

classmethod get_file_seg_type()

Type of the last segment in the path.

class seismicrna.idmut.io.RefseqIO(*args, refseq: DNA, **kwargs)

Bases: RefBrickleIO, IDmutIO

classmethod get_file_seg_type()

Type of the last segment in the path.

property refseq
seismicrna.idmut.io.from_reads(reads: Iterable[tuple[str, tuple[tuple[list[int], list[int]], dict[int, int]]]], *, sample: str, branches: dict[str, str], ref: str, refseq: DNA, batch: int, write_read_names: bool, drop_empty_reads: bool = True)

Gather reads into a batch of relationships.

class seismicrna.idmut.lists.IDmutList(*, sample: str, branches: Iterable[str], ref: str, data: pd.DataFrame, **kwargs)

Bases: List, IDmutFile, ABC

class seismicrna.idmut.lists.IDmutPositionList(*, sample: str, branches: Iterable[str], ref: str, data: pd.DataFrame, **kwargs)

Bases: PositionList, IDmutList

classmethod get_table_type()

Type of table that this type of list can process.

seismicrna.idmut.main.check_duplicates(xam_files: list[Path])

Check if any combination of sample, reference, and branches occurs more than once.

seismicrna.idmut.main.run(fasta: str | Path = Sentinel.UNSET, input_path: Iterable[str | Path] = Sentinel.UNSET, *, out_dir: str | Path = './out', branch: str = '', min_reads: int = 1000, min_mapq: int = 25, phred_enc: int = 33, min_phred: int = 25, batch_size: int = 65536, insert3: bool = True, ambindel: bool = True, overhangs: bool = True, clip_end5: int = 4, clip_end3: int = 4, sep_strands: bool = False, rev_label: str = '-rev', write_read_names: bool = False, idmut_pos_table: bool = True, idmut_read_table: bool = False, idmut_cx: bool = True, num_cpus: int = 4, brotli_level: int = 10, force: bool = False, keep_tmp: bool = False, tmp_pfx='./tmp')

Compute relationships between references and aligned reads.

Parameters:

  • out_dir (str | pathlib._local.Path) – Write all output files to this directory [keyword-only, default: ‘./out’]

  • branch (str) – Create a new branch of the workflow with this name [keyword-only, default: ‘’]

  • min_reads (int) – Discard alignment maps with fewer than this many reads [keyword-only, default: 1000]

  • min_mapq (int) – Discard reads with mapping qualities below this threshold [keyword-only, default: 25]

  • phred_enc (int) – Specify the Phred score encoding of FASTQ and SAM/BAM/CRAM files [keyword-only, default: 33]

  • min_phred (int) – Mark base calls with Phred scores lower than this threshold as ambiguous [keyword-only, default: 25]

  • batch_size (int) – Limit batches to at most this many reads [keyword-only, default: 65536]

  • insert3 (bool) – Mark each insertion on the base to its 3’ (True) or 5’ (False) side [keyword-only, default: True]

  • ambindel (bool) – Mark all ambiguous insertions and deletions (indels) [keyword-only, default: True]

  • overhangs (bool) – Retain the overhangs of paired-end mates that dovetail [keyword-only, default: True]

  • clip_end5 (int) – Clip this many bases from the 5’ end of each read [keyword-only, default: 4]

  • clip_end3 (int) – Clip this many bases from the 3’ end of each read [keyword-only, default: 4]

  • sep_strands (bool) – Separate each alignment map into forward- and reverse-strand reads [keyword-only, default: False]

  • rev_label (str) – With –sep-strands, add this label to each reverse-strand reference [keyword-only, default: ‘-rev’]

  • write_read_names (bool) – Write the name of each read in a second set of batches (necessary for the options –drop-read or –drop-read-file) [keyword-only, default: False]

  • idmut_pos_table (bool) – Tabulate relationships per position for idmut data [keyword-only, default: True]

  • idmut_read_table (bool) – Tabulate relationships per read for idmut data [keyword-only, default: False]

  • idmut_cx (bool) – Use a fast (C extension module) version of the idmut algorithm; the slow (Python) version is still avilable as a fallback if the C extension cannot be loaded, and for debugging/benchmarking [keyword-only, default: True]

  • num_cpus (int) – Use up to this many CPUs simultaneously [keyword-only, default: 4]

  • brotli_level (int) – Compress pickle files with this level of Brotli (0 - 11) [keyword-only, default: 10]

  • force (bool) – Force all tasks to run, overwriting any existing output files [keyword-only, default: False]

  • keep_tmp (bool) – Keep temporary files after finishing [keyword-only, default: False]

  • tmp_pfx – Write all temporary files to a directory with this prefix [keyword-only, default: ‘./tmp’]

class seismicrna.idmut.report.BaseIDmutReport(**kwargs: Any | Callable[[Report], Any])

Bases: RefReport, IDmutIO, ABC

classmethod get_file_seg_type()

Type of the last segment in the path.

class seismicrna.idmut.report.IDmutReport(**kwargs: Any | Callable[[Report], Any])

Bases: BatchedReport, BaseIDmutReport

classmethod get_checksum_report_fields()

Checksum fields of the report.

classmethod get_param_report_fields()

Parameter fields of the report.

get_refseq_file(top: Path)
classmethod get_result_report_fields()

Result fields of the report.

class seismicrna.idmut.report.PoolReport(**kwargs: Any | Callable[[Report], Any])

Bases: BaseIDmutReport

classmethod get_param_report_fields()

Parameter fields of the report.

class seismicrna.idmut.sam.SamFileViewer(xam_input: Path, tmp_dir: Path, branch: str, batch_size: int, num_cpus: int = 1)

Bases: object

property ancestors
property branches
create_tmp_sam()

Create the temporary SAM file.

delete_tmp_sam()

Delete the temporary SAM file.

property flagstats
property indexes
iter_records(batch: int)

Iterate through the records of the batch.

property n_reads

Total number of reads.

property paired

Whether the reads are paired.

property ref
property sample
property tmp_sam_path

Get the path to the temporary SAM file.

seismicrna.idmut.sam.get_line_attrs(line: str) → tuple[str, bool, bool]

Read attributes from a line in a SAM file.

seismicrna.idmut.sam.tmp_xam_cmd(xam_in: Path, xam_out: Path, paired: bool, num_cpus: int = 1)

Collate and create a temporary XAM file.

seismicrna.idmut.sim.calc_pmut_pattern(pmut: pd.DataFrame, pattern: RelPattern, normalize: bool = True)

Calculate the rate of a given type of mutation.

Parameters:

normalize (bool) – If True (default), divide by (MATCH + fmut) to get the conditional mutation rate among informative positions. If False, return the unconditional probability sum fmut directly — use this when computing p_noclose for simulation, where non-informative (AMB) positions are treated as non-mutations by reads_noclose_muts.

seismicrna.idmut.sim.make_p_ends_2d(pends: np.ndarray, uniq_end5s: np.ndarray, uniq_end3s: np.ndarray, pmut_index) → np.ndarray

Convert a 1-D end-pair probability array to a 2-D (positions × positions) matrix required by calc_p_noclose_given_clust.

seismicrna.idmut.sim.simulate_batch(sample: str, branches: dict[str, str], ref: str, batch: int, write_read_names: bool, formatter: ~typing.Callable[[int, int], str] = <function format_read_name>, **kwargs)

Simulate a pair of IDmutBatchIO and ReadNamesBatchIO.

seismicrna.idmut.sim.simulate_batches(batch_size: int, pmut: pd.DataFrame, pclust: pd.Series, pends: np.ndarray, uniq_end5s: np.ndarray, uniq_end3s: np.ndarray, num_reads: int, min_mut_gap_weights: dict[int, float], injected_mut_probs: dict[int, float], mut_collisions: str, seed: int | None, **kwargs)

Simulate batches of reads for all clusters.

Parameters:

  • batch_size (int) – Number of reads per batch.

  • pmut (pd.DataFrame) – Mutation rate DataFrame with columns indexed by (rel, k, clust).

  • pclust (pd.Series) – Proportion of reads belonging to each cluster; index is (k, clust).

  • num_reads (int) – Total number of reads to simulate across all clusters.

  • seed (int | None) – Random seed for reproducibility; None for no fixed seed.

  • min_mut_gap_weights (dict[int, float]) – Mapping of min_mut_gap value to weight. When non-empty, reads for each cluster are split across the gap values proportionally. When empty, the gap is derived from the maximum offset in injected_mut_probs (or 0 if that is also empty).

  • injected_mut_probs (dict[int, float]) – Mapping of 5’ offset to injection probability. Passed to each simulate_cluster call and used to derive the gap when min_mut_gap_weights is empty.

  • mut_collisions (str) – How to handle reads with mutations closer than the gap.

  • **kwargs – Additional keyword arguments forwarded to simulate_cluster.

Yields:

tuple[IDmutBatchIO, ReadNamesBatchIO | None] – Pairs of idmut and (optionally) read-name batch objects.

seismicrna.idmut.sim.simulate_cluster(first_batch: int, batch_size: int, num_reads: int, seed: int | None, **kwargs)

Simulate all batches for one cluster.

seismicrna.idmut.sim.simulate_idmut(*, out_dir: Path, tmp_dir: Path, branch: str, sample: str, ref: str, refseq: DNA, write_read_names: bool, brotli_level: int, force: bool, **kwargs)

Simulate an entire IDmut step.

seismicrna.idmut.strands.generate_both_strands(ref: str, seq: DNA, rev_label: str)

Yield both the forward and reverse strand for each sequence.

seismicrna.idmut.strands.write_both_strands(fasta_in: Path, fasta_out: Path, rev_label: str)

Write a FASTA file of both forward and reverse strands.

class seismicrna.idmut.table.AverageTable

Bases: RelTypeTable, ABC

Average over an ensemble of RNA structures.

classmethod get_header_type()

Type of the header for the table.

class seismicrna.idmut.table.AverageTabulator(*, top: Path, branches: dict[str, str], sample: str, region: Region, count_ends: bool, count_pos: bool, count_read: bool, validate: bool = True)

Bases: Tabulator, ABC

property data_per_clust

Series of per-cluster data (or None if no clusters).

class seismicrna.idmut.table.FullTabulator(*, ref: str, refseq: DNA, count_ends: bool = False, **kwargs)

Bases: Tabulator, ABC

classmethod get_null_value()

The null value for a count: either 0 or NaN.

class seismicrna.idmut.table.IDmutBatchTabulator(*, get_batch_count_all: Callable, num_batches: int, num_cpus: int = 1, **kwargs)

Bases: BatchTabulator, IDmutTabulator

class seismicrna.idmut.table.IDmutCountTabulator(*, batch_counts: Iterable[tuple[Any, Any, Any, Any]], **kwargs)

Bases: CountTabulator, IDmutTabulator

class seismicrna.idmut.table.IDmutDatasetTabulator(*, dataset: MutsDataset, validate: bool = False, **kwargs)

Bases: DatasetTabulator, IDmutTabulator

classmethod init_kws()

Attributes of the dataset to use as keyword arguments in super().__init__().

class seismicrna.idmut.table.IDmutPositionTable

Bases: IDmutTable, PositionTable, ABC

class seismicrna.idmut.table.IDmutPositionTableLoader(table_file: str | Path, **kwargs)

Bases: PositionTableLoader, IDmutPositionTable

Load IDmut data indexed by position.

class seismicrna.idmut.table.IDmutPositionTableWriter(tabulator: Tabulator)

Bases: PositionTableWriter, IDmutPositionTable

class seismicrna.idmut.table.IDmutReadTable

Bases: IDmutTable, ReadTable, ABC

class seismicrna.idmut.table.IDmutReadTableLoader(table_file: str | Path, **kwargs)

Bases: ReadTableLoader, IDmutReadTable

Load IDmut data indexed by read.

class seismicrna.idmut.table.IDmutReadTableWriter(tabulator: Tabulator)

Bases: ReadTableWriter, IDmutReadTable

class seismicrna.idmut.table.IDmutTable

Bases: AverageTable, IDmutFile, ABC

classmethod get_load_function()

LoadFunction for all Dataset types for this Table.

class seismicrna.idmut.table.IDmutTabulator(*, ref: str, refseq: DNA, count_ends: bool = False, **kwargs)

Bases: FullTabulator, AverageTabulator, ABC

classmethod table_types()

Types of tables that this tabulator can write.

class seismicrna.idmut.write.RelationWriter(sam_view: SamFileViewer, fasta_file: str | Path)

Bases: object

Compute and write relationships for all reads from one sample aligned to one reference sequence.

property branches
property num_reads
property ref
property refseq
property sample
write(*, out_dir: Path, release_dir: Path, min_mapq: int, min_reads: int, min_phred: int, phred_enc: int, insert3: bool, ambindel: bool, overhangs: bool, clip_end5: int, clip_end3: int, idmut_pos_table: bool, idmut_read_table: bool, brotli_level: int, force: bool, num_cpus: int, **kwargs)

Compute relationships for every record in a XAM file.

seismicrna.idmut.write.generate_batch(batch: int, *, sam_view: SamFileViewer, top: Path, refseq: DNA, brotli_level: int, count_pos: bool, count_read: bool, write_read_names: bool, **kwargs)

Compute relationships for every SAM record in one batch.

seismicrna.idmut.write.idmut_records(records: Iterable[tuple[str, str, str]], ref: str, refseq: str, min_mapq: int, min_qual: int, insert3: bool, ambindel: bool, overhangs: bool, clip_end5: int, clip_end3: int, idmut_cx: bool)

Yield relationships for each SAM record in an iterable.

Parameters:

  • records (Iterable[tuple[str, str, str]]) – Iterable of (name, line1, line2) tuples from a SAM file, where line2 is empty for single-end reads.

  • ref (str) – Reference name expected in each SAM record.

  • refseq (str) – Full reference sequence string.

  • min_mapq (int) – Minimum acceptable mapping quality score.

  • min_qual (int) – Minimum Phred quality score (as an integer) for base calls.

  • insert3 (bool) – Whether to mark insertions on the 3’ flanking reference position (True) or the 5’ position (False).

  • ambindel (bool) – Whether to find and label ambiguous indel positions.

  • overhangs (bool) – Whether to allow paired-end mates to overhang one another.

  • clip_end5 (int) – Number of bases to clip from the 5’ end of each read.

  • clip_end3 (int) – Number of bases to clip from the 3’ end of each read.

  • idmut_cx (bool) – Whether to use the C extension for the idmut algorithm; falls back to the Python implementation if import fails.

Yields:

tuple[str, tuple] – Read name and the result of id_muts_lines for that read.

seismicrna.idmut.write.idmut_xam(xam_file: Path, *, fasta: Path, tmp_dir: Path, branch: str, batch_size: int, num_cpus: int, **kwargs)

Write the batches of relationships for one XAM file.