Scripts Reference¶

Documentation for Python scripts in workflow/scripts/.

Overview¶

transfer_tags.py¶

Transfer specified BAM tags from one file to another with optional renaming.

Usage¶

1
2
3
4
5
6

python transfer_tags.py \
    --tags ML MM \
    --rename ML=CL MM=CM \
    --source source.bam \
    --target target.bam \
    --output output.bam

Arguments¶

Behavior¶

1. Reads source BAM, stores specified tags by read ID
2. Iterates target BAM, transfers matching tags
3. Outputs only primary alignments (filters secondary/supplementary)

get_charging_table.py¶

Extract charging likelihood tag values from BAM to TSV format.

Usage¶

1

python get_charging_table.py --tag CL input.bam output.tsv.gz

Arguments¶

Output Format¶

1
2

read_id                                  tRNA              charging_likelihood
00a1b2c3-4567-89ab-cdef-0123456789ab    tRNA-Ala-AGC-1-1  245

get_trna_charging_cpm.py¶

Aggregate per-read charging data to tRNA-level CPM counts.

Usage¶

1
2
3
4

python get_trna_charging_cpm.py \
    --input charging_prob.tsv.gz \
    --output charging.cpm.tsv.gz \
    --ml-threshold 200

Arguments¶

Classification Logic¶

1
2
3
4

if charging_likelihood >= threshold:
    classification = "charged"
else:
    classification = "uncharged"

Output Format¶

get_charging_summary.py¶

Generate tRNA-gene level charging statistics.

Usage¶

1
2
3
4

python get_charging_summary.py \
    -b sample.bam \
    -t trna_table.txt \
    -o output.tsv

Arguments¶

Lookup Table Format¶

4 whitespace-delimited columns (no header):

1
2

uncharged_seq    charged_seq    isodecoder    gene_id
tRNA-Ala-AGC-u   tRNA-Ala-AGC-c Ala-AGC       tRNA-Ala-1

Output Format¶

get_bcerror_freqs.py¶

Calculate per-position basecalling error frequencies.

Usage¶

1

python get_bcerror_freqs.py sample.bam reference.fa output.tsv.gz

Arguments¶

Metrics Calculated¶

Output Format¶

get_align_stats.py¶

Summarize alignment statistics across multiple BAM files.

Usage¶

1
2
3
4
5

python get_align_stats.py \
    -o stats.tsv.gz \
    -a unmapped aligned classified \
    -i sample_id \
    -b ubam.bam aligned.bam classified.bam

Arguments¶

Classes¶

CountArray: Memory-efficient integer count tracking

• push(index): Increment count at index
• nth(n): Get nth value
• mean(): Calculate mean
• quantile(p): Calculate percentile

ReadStats: Aggregate read metrics

• Tracks read length, base quality, MAPQ distributions
• Returns summary with mean/median for each metric

Output Format¶

extract_signal_metrics.py¶

Extract signal metrics using Remora API for model interpretation.

Usage¶

1
2
3
4
5
6
7
8

python extract_signal_metrics.py \
    --pod5_dir sample.pod5 \
    --bam aligned.bam \
    --kmer 9mer_levels.txt \
    --sample_name sample1 \
    [--region chr:start-end] \
    [--bed regions.bed] \
    | gzip > output.tsv.gz

Arguments¶

Dependencies¶

Requires custom Remora fork:

1

pip install git+https://github.com/rnabioco/remora.git@metrics_missing_ok

Output Format¶

filter_reads.py¶

Filter BAM alignments by quality criteria.

Usage¶

1
2
3
4
5
6

python filter_reads.py \
    -i input.bam \
    -o output.bam \
    -5 24 \
    -3 23 \
    -s

Arguments¶

Filter Flags¶

The script adds a zf tag with bitwise filter flags:

FilterStats Class¶

Tracks filtering reasons with methods:

• add(reason): Record filter reason
• summary(): Get filter statistics
• write_stats(path): Write statistics to file

generate_squiggy_session.py¶

Generate a Squiggy session JSON file for loading pipeline outputs in Positron IDE.

Usage¶

1
2
3
4
5
6
7

python generate_squiggy_session.py \
    --samples sample1 sample2 \
    --output-dir /path/to/output \
    --fasta /path/to/reference.fa \
    --output squiggy-session.json \
    [--session-name "My session"] \
    [--no-checksums]

Arguments¶

Output Format¶

JSON with the following structure:

compute_odds_ratios.py¶

Compute per-tRNA pairwise modification odds ratios from modkit calls and charging data.

Usage¶

1
2
3
4
5
6

python compute_odds_ratios.py \
    --modkit mod_calls.tsv.gz \
    --charging charging_prob.tsv.gz \
    --output odds_ratios.tsv.gz \
    --ml-threshold 200 \
    --min-coverage 10

Arguments¶

Behavior¶

1. Loads modkit per-read modification calls and binarizes (canonical vs modified)
2. Loads charging probabilities and binarizes at threshold
3. For each tRNA, pivots to per-read matrix (rows=reads, columns=positions)
4. Adds charging as position 999
5. Tests all pairwise combinations via 2x2 contingency tables
6. Applies Haldane correction for zero cells, Fisher's exact test
7. Applies BH correction across all results

Output Format¶

compute_seq_similarity.py¶

Compute pairwise sequence similarity matrix for reference FASTA sequences.

Usage¶

1

python compute_seq_similarity.py reference.fa output.tsv

Arguments¶

Behavior¶

1. Reads all sequences from FASTA using pysam
2. Computes all-vs-all pairwise Needleman-Wunsch global alignment using parasail
3. Calculates percent identity = matches / max(len_seq1, len_seq2) * 100
4. Writes square TSV similarity matrix

Output Format¶

Square TSV matrix with sequence names as row and column headers. Values are percent identity (0-100). Diagonal entries are always 100.0.

collapse_gtrndb_fasta.py¶

Collapse redundant sequences in a GtRNAdb FASTA reference.

Usage¶

1
2
3
4
5

python collapse_gtrndb_fasta.py \
    -i ecoliK12-mature-tRNAs.fa \
    -o ecoliK12-collapsed.fa \
    -m ecoliK12-mapping.tsv \
    [--keep-unparsed]

Arguments¶

Behavior¶

1. Parses GtRNAdb FASTA headers using end-anchored regex to extract tRNA identity
2. Groups sequences by isodecoder family (amino acid + anticodon)
3. Strips trailing CCA for sequence comparison
4. Identifies redundant gene copies with identical sequences
5. Writes collapsed FASTA (first occurrence as representative) and TSV mapping

Mapping Output Format¶

Common Patterns¶

Reading BAM Files¶

1
2
3
4
5
6
7
8
9

import pysam

with pysam.AlignmentFile(bam_path, "rb") as bam:
    for read in bam.fetch():
        # Access tags
        tag_value = read.get_tag("CL")
        # Access alignment info
        ref_name = read.reference_name
        read_id = read.query_name

Writing Compressed Output¶

1
2
3
4
5

import gzip

with gzip.open(output_path, "wt") as f:
    f.write("column1\tcolumn2\n")
    f.write(f"{value1}\t{value2}\n")

Using pandas¶

1
2
3
4

import pandas as pd

df = pd.read_csv(input_path, sep="\t", compression="gzip")
df.to_csv(output_path, sep="\t", index=False, compression="gzip")