NGSTk API Documentation

Overview

NGSTk (Next-Generation Sequencing Toolkit) is a toolkit class that provides helper functions for building command strings used in NGS pipelines. It can be configured with a YAML configuration file to specify custom tool paths, or it will use tools from the system PATH.

Key Features

• Command Building: Generate command strings for common NGS tools
• Configuration Management: Use custom tool paths via YAML config
• Tool Integration: Built-in support for common tools like samtools, bedtools, etc.
• Pipeline Integration: Works seamlessly with PipelineManager

Installation

NGSTk is included with pypiper:

pip install pypiper

Quick Example

from pypiper.ngstk import NGSTk

# Initialize NGSTk
tk = NGSTk()

# Generate a command
cmd = tk.samtools_index("sample.bam")
# Returns: "samtools index sample.bam"

API Reference

NGSTk Class

NGSTk

class NGSTk(config_file: str | None = None, pm: Any = None)

Build shell command strings for common NGS processing operations.

Example

tk = NGSTk() tk.samtools_index("sample.bam") # => "samtools index sample.bam"

tk = NGSTk("pipeline_config.yaml") tk.samtools_index("sample.bam") # uses configured samtools path

Tool paths come from the config's "tools" section; unconfigured tools default to their name (assuming they're on $PATH).

Parameters:

• config_file (str | None): Path to pipeline YAML config file. (default: None)
• pm (Any): PipelineManager to associate with this toolkit. (default: None)

Attributes

• parameters
• pm
• tools
• ziptool_cmd

Properties

ziptool

@property
def ziptool -> str

Compression command: 'pigz' if available with multiple cores, else 'gzip'.

Methods

add_track_to_hub

def add_track_to_hub(sample_name: str, track_url: str, track_hub: str, colour: str, five_prime: str = '') -> list[str]

bam2fastq

def bam2fastq(input_bam: str, output_fastq: str, output_fastq2: str | None = None, unpaired_fastq: str | None = None) -> str

Build command to convert BAM to FASTQ via Picard SamToFastq.

bam_conversions

def bam_conversions(bam_file: str, depth: bool = True) -> str

Build command to sort and index a BAM file (optionally with depth).

Parameters:

• bam_file (str): Path to BAM file.
• depth (bool): Also calculate per-position coverage. (default: True)

bam_to_bed

def bam_to_bed(input_bam: str, output_bed: str) -> str

bam_to_bigwig

def bam_to_bigwig(input_bam: str, output_bigwig: str, genome_sizes: str, genome: str, tagmented: bool = False, normalize: bool = False, norm_factor: int = 1000) -> list[str]

Convert a BAM file to a bigWig file.

Parameters:

• input_bam (str): path to BAM file to convert
• output_bigwig (str): path to which to write file in bigwig format
• genome_sizes (str): path to file with chromosome size information
• genome (str): name of genomic assembly
• tagmented (bool): flag related to read-generating protocol (default: False)
• normalize (bool): whether to normalize coverage (default: False)
• norm_factor (int): number of bases to use for normalization (default: 1000)

Returns:

• list[str]: list[str]: sequence of commands to execute

bam_to_fastq

def bam_to_fastq(bam_file: str, out_fastq_pre: str, paired_end: bool) -> str

Build Picard SamToFastq command for BAM to FASTQ conversion.

bam_to_fastq_awk

def bam_to_fastq_awk(bam_file: str, out_fastq_pre: str, paired_end: bool, zipmode: bool = False) -> tuple[str, str, str | None]

Build fast awk-based BAM to FASTQ conversion command.

Faster than Picard/bedtools but assumes paired reads are properly ordered with no singletons.

bam_to_fastq_bedtools

def bam_to_fastq_bedtools(bam_file: str, out_fastq_pre: str, paired_end: bool) -> tuple[str, str, str | None]

Build bedtools bamtofastq command for BAM to FASTQ conversion.

bowtie2_map

def bowtie2_map(input_fastq1: str, output_bam: str, log: str, metrics: str, genome_index: str, max_insert: int | str, cpus: int | str, input_fastq2: str | None = None) -> str

calc_frip

def calc_frip(input_bam: str, input_bed: str, threads: int = 4) -> str

Calculate fraction of reads in peaks.

A file of with a pool of sequencing reads and a file with peak call regions define the operation that will be performed. Thread count for samtools can be specified as well.

Parameters:

• input_bam (str): sequencing reads file
• input_bed (str): file with called peak regions
• threads (int): number of threads samtools may use (default: 4)

Returns:

• float (str): fraction of reads in peaks defined in given peaks file

calculate_frip

def calculate_frip(input_bam: str, input_bed: str, output: str, cpus: int = 4) -> str

center_peaks_on_motifs

def center_peaks_on_motifs(peak_file: str, genome: str, window_width: int, motif_file: str, output_bed: str) -> str

check_command

def check_command(command: str) -> bool

Check if a command is callable on the system.

check_fastq

def check_fastq(input_files: str | list[str], output_files: str | list[str], paired_end: bool) -> Callable

Return a follow function that validates FASTQ conversion read counts.

Example

cmd, prefix, out = tk.input_to_fastq(bam, name, pe, folder) pm.run(cmd, out, follow=tk.check_fastq(bam, out, pe))

Parameters:

• input_files (str | list[str]): Original input file(s) before conversion.
• output_files (str | list[str]): FASTQ output file(s) from conversion.
• paired_end (bool): Whether data is paired-end.

Returns:

• Callable: Callable that compares read counts and reports stats.

check_trim

def check_trim(trimmed_fastq: str, paired_end: bool, trimmed_fastq_R2: str | None = None, fastqc_folder: str | None = None) -> Callable

Return a follow function that counts trimmed reads and optionally runs FastQC.

Parameters:

• trimmed_fastq (str): Path to trimmed reads file.
• paired_end (bool): Whether data is paired-end.
• trimmed_fastq_R2 (str | None): Path to R2 trimmed file for paired-end. (default: None)
• fastqc_folder (str | None): If set, run FastQC and place output here. (default: None)

Returns:

• Callable: Callable for use as pm.run() follow function.

count_concordant

def count_concordant(aligned_bam: str) -> str

Count reads aligned concordantly exactly once (YT:Z:CP flag).

count_fail_reads

def count_fail_reads(file_name: str, paired_end: bool) -> int

Count reads that failed platform/vendor quality checks (SAM flag 512).

count_flag_reads

def count_flag_reads(file_name: str, flag: int | str, paired_end: bool) -> str

Count reads with a specific SAM flag value.

count_lines

def count_lines(file_name: str) -> str

Count lines in a file using wc -l.

count_lines_zip

def count_lines_zip(file_name: str) -> str

Count lines in a gzipped file using zcat | wc -l.

count_mapped_reads

def count_mapped_reads(file_name: str, paired_end: bool) -> str | int

Count mapped reads in a BAM/SAM file (excludes unmapped, flag -F4).

count_multimapping_reads

def count_multimapping_reads(file_name: str, paired_end: bool) -> int

Count reads flagged as multimapping (SAM flag 256).

count_reads

def count_reads(file_name: str, paired_end: bool) -> str | int | float

Count reads in a BAM/SAM/FASTQ file.

Paired-end reads count as 2. Assumes paired-end FASTQs are split into separate R1/R2 files (divides line count by 2 instead of 4).

Parameters:

• file_name (str): Path to BAM/SAM/FASTQ file.
• paired_end (bool): Whether the file contains paired-end reads.

count_unique_mapped_reads

def count_unique_mapped_reads(file_name: str, paired_end: bool) -> int

Count mapped reads (by name, deduplicated) in a BAM/SAM file.

count_unique_reads

def count_unique_reads(file_name: str, paired_end: bool) -> int

Count unique reads (by name) in a BAM/SAM file. Paired-end counts as 2.

count_uniquelymapping_reads

def count_uniquelymapping_reads(file_name: str, paired_end: bool) -> str

Count reads that mapped to a unique position (exclude flag 256).

fastqc

def fastqc(file: str, output_dir: str) -> str

Build FastQC command for a reads file.

fastqc_rename

def fastqc_rename(input_bam: str, output_dir: str, sample_name: str) -> list[str]

Build commands to run FastQC and rename output by sample name.

filter_peaks_mappability

def filter_peaks_mappability(peaks: str, alignability: str, filtered_peaks: str) -> str

filter_reads

def filter_reads(input_bam: str, output_bam: str, metrics_file: str, paired: bool = False, cpus: int = 16, Q: int = 30) -> list[str]

Build commands to dedup, quality-filter, and remove multimappers.

genome_wide_coverage

def genome_wide_coverage(input_bam: str, genome_windows: str, output: str) -> str

get_chrs_from_bam

def get_chrs_from_bam(file_name: str) -> list[str]

Extract chromosome names from a BAM file header via samtools.

get_file_size

def get_file_size(filenames: str | list[str]) -> float

Get total size of file(s) in megabytes.

get_fragment_sizes

def get_fragment_sizes(bam_file: str) -> Any

get_frip

def get_frip(sample: Any) -> dict

Calculates the fraction of reads in peaks for a given sample.

Parameters:

• sample (pipelines.Sample): Sample object with "peaks" attribute.

get_input_ext

def get_input_ext(input_file: str) -> str

Detect input file type: ".bam", ".fastq.gz", or ".fastq".

get_mitochondrial_reads

def get_mitochondrial_reads(bam_file: str, output: str, cpus: int = 4) -> list[str]

get_peak_number

def get_peak_number(sample: Any) -> Any

Counts number of peaks from a sample's peak file.

Parameters:

• sample (pipelines.Sample): Sample object with "peaks" attribute.

get_read_type

def get_read_type(bam_file: str, n: int = 10) -> tuple[str, int]

Gets the read type (single, paired) and length of bam file.

Parameters:

• bam_file (str): Bam file to determine read attributes.
• n (int): Number of lines to read from bam file. (default: 10)

Returns:

• tuple[str, int]: str, int: tuple of read type and read length

homer_annotate_pPeaks

def homer_annotate_pPeaks(peak_file: str, genome: str, motif_file: str, output_bed: str) -> str

homer_find_motifs

def homer_find_motifs(peak_file: str, genome: str, output_dir: str, size: int = 150, length: str = '8,10,12,14,16', n_motifs: int = 12) -> str

htseq_count

def htseq_count(input_bam: str, gtf: str, output: str) -> list[str]

index_bam

def index_bam(input_bam: str) -> str

input_to_fastq

def input_to_fastq(input_file: str | list[str], sample_name: str, paired_end: bool, fastq_folder: str, output_file: str | None = None, multiclass: bool = False, zipmode: bool = False) -> list

Build command to convert any input (.bam/.fastq.gz/.fastq) to FASTQ.

Parameters:

• input_file (str | list[str]): Path(s) to input file(s).
• sample_name (str): Sample name for output file naming.
• paired_end (bool): Whether data is paired-end.
• fastq_folder (str): Directory for output FASTQ files.
• output_file (str | None): Explicit output path (auto-derived if None). (default: None)
• multiclass (bool): Internal flag for recursive R1/R2 handling. (default: False)
• zipmode (bool): Output as .fastq.gz instead of .fastq. (default: False)

Returns:

• list: List of [command, fastq_prefix, output_file].

kallisto

def kallisto(input_fastq: str, output_dir: str, output_bam: str, transcriptome_index: str, cpus: int, input_fastq2: str | None = None, size: int = 180, b: int = 200) -> list[str]

link_to_track_hub

def link_to_track_hub(track_hub_url: str, file_name: str, genome: str) -> None

macs2_call_peaks

def macs2_call_peaks(treatment_bams: str | Iterable[str], output_dir: str, sample_name: str, genome: str, control_bams: str | Iterable[str] | None = None, broad: bool = False, paired: bool = False, pvalue: float | None = None, qvalue: float | None = None, include_significance: bool | None = None) -> str

Use MACS2 to call peaks.

Parameters:

• treatment_bams (str | Iterable[str]): Paths to files with data to regard as treatment.
• output_dir (str): Path to output folder.
• sample_name (str): Name for the sample involved.
• genome (str): Name of the genome assembly to use.
• control_bams (str | Iterable[str]): Paths to files with data to regard as control (default: None)
• broad (bool): Whether to do broad peak calling. (default: False)
• paired (bool): Whether reads are paired-end (default: False)
• pvalue (float | NoneType): Statistical significance measure to pass as --pvalue to peak calling with MACS (default: None)
• qvalue (float | NoneType): Statistical significance measure to pass as --qvalue to peak calling with MACS (default: None)
• include_significance (bool | NoneType): Whether to pass a statistical significance argument to peak calling with MACS; if omitted, this will be True if the peak calling is broad or if either p-value or q-value is specified; default significance specification is a p-value of 0.001 if a significance is to be specified but no value is provided for p-value or q-value. (default: None)

Returns:

• str: Command to run.

macs2_call_peaks_atacseq

def macs2_call_peaks_atacseq(treatment_bam: str, output_dir: str, sample_name: str, genome: str) -> str

macs2_plot_model

def macs2_plot_model(r_peak_model_file: str, sample_name: str, output_dir: str) -> list[str]

make_dir

def make_dir(path: str) -> None

Create directory and all intermediates, no error if exists.

make_sure_path_exists

def make_sure_path_exists(path: str) -> None

Alias for make_dir

mark_duplicates

def mark_duplicates(aligned_file: str, out_file: str, metrics_file: str, remove_duplicates: str = 'True') -> str

merge_bams

def merge_bams(input_bams: list[str], merged_bam: str, in_sorted: bool | str = 'TRUE', tmp_dir: str | None = None) -> str | int

Build Picard MergeSamFiles command to combine BAM files.

merge_bams_samtools

def merge_bams_samtools(input_bams: list[str], merged_bam: str) -> str

merge_fastq

def merge_fastq(inputs: list[str], output: str, run: bool = False, remove_inputs: bool = False) -> str | None

Merge multiple FASTQ files into one via cat.

merge_or_link

def merge_or_link(input_args: list, raw_folder: str, local_base: str = 'sample') -> str | list[str]

Standardize inputs by linking or merging .bam/.fastq/.fastq.gz files.

Example

local = tk.merge_or_link([["s1_R1.fq.gz"], ["s1_R2.fq.gz"]], "raw/", "sample1")

For single files, creates a symlink. For multiple files of the same type, merges them (cat for fastq, samtools merge for bam).

Parameters:

• input_args (list): List of input file paths or list of lists (R1/R2).
• raw_folder (str): Directory for the merge/link output.
• local_base (str): Base name for output file (usually sample name). (default: 'sample')

move_file

def move_file(old: str, new: str) -> str

parse_bowtie_stats

def parse_bowtie_stats(stats_file: str) -> dict

Parses Bowtie2 stats file, returns dict with values.

Parameters:

• stats_file (str): Bowtie2 output file with alignment statistics.

parse_duplicate_stats

def parse_duplicate_stats(stats_file: str) -> dict

Parses sambamba markdup output, returns dict with values.

Parameters:

• stats_file (str): sambamba output file with duplicate statistics.

parse_qc

def parse_qc(qc_file: str) -> dict

Parse phantompeakqualtools (spp) QC table and return quality metrics.

Parameters:

• qc_file (str): Path to phantompeakqualtools output file, which contains sample quality measurements.

picard_mark_duplicates

def picard_mark_duplicates(input_bam: str, output_bam: str, metrics_file: str, temp_dir: str = '.') -> list[str]

plot_atacseq_insert_sizes

def plot_atacseq_insert_sizes(bam: str, plot: str, output_csv: str, max_insert: int = 1500, smallest_insert: int = 30) -> None

Heavy inspiration from here: https://github.com/dbrg77/ATAC/blob/master/ATAC_seq_read_length_curve_fitting.ipynb

preseq_coverage

def preseq_coverage(bam_file: str, output_prefix: str) -> str

preseq_curve

def preseq_curve(bam_file: str, output_prefix: str) -> str

preseq_extrapolate

def preseq_extrapolate(bam_file: str, output_prefix: str) -> str

remove_file

def remove_file(file_name: str) -> str

run_spp

def run_spp(input_bam: str, output: str, plot: str, cpus: int) -> str

Run the SPP read peak analysis tool.

Parameters:

• input_bam (str): Path to reads file
• output (str): Path to output file
• plot (str): Path to plot file
• cpus (int): Number of processors to use

Returns:

• str: Command with which to run SPP

sam_conversions

def sam_conversions(sam_file: str, depth: bool = True) -> str

Build command to convert SAM to sorted/indexed BAM (optionally with depth).

Parameters:

• sam_file (str): Path to SAM file.
• depth (bool): Also calculate per-position coverage. (default: True)

sambamba_remove_duplicates

def sambamba_remove_duplicates(input_bam: str, output_bam: str, cpus: int = 16) -> str

samtools_index

def samtools_index(bam_file: str) -> str

Index a bam file.

samtools_view

def samtools_view(file_name: str, param: str, postpend: str = '') -> str

Run samtools view with given parameters and optional post-processing pipe.

shift_reads

def shift_reads(input_bam: str, genome: str, output_bam: str) -> str

simple_frip

def simple_frip(input_bam: str, input_bed: str, threads: int = 4) -> str

skewer

def skewer(input_fastq1: str, output_prefix: str, output_fastq1: str, log: str, cpus: int | str, adapters: str, input_fastq2: str | None = None, output_fastq2: str | None = None) -> list[str]

Build skewer adapter-trimming commands with file renaming.

slurm_footer

def slurm_footer() -> str

slurm_header

def slurm_header(job_name: str, output: str, queue: str = 'shortq', n_tasks: int = 1, time: str = '10:00:00', cpus_per_task: int = 8, mem_per_cpu: int = 2000, nodes: int = 1, user_mail: str = '', mail_type: str = 'end') -> str

slurm_submit_job

def slurm_submit_job(job_file: str) -> int

sort_index_bam

def sort_index_bam(input_bam: str, output_bam: str) -> list[str]

spp_call_peaks

def spp_call_peaks(treatment_bam: str, control_bam: str, treatment_name: str, control_name: str, output_dir: str, broad: str | bool, cpus: int, qvalue: float | None = None) -> str

Build command for R script to call peaks with SPP.

Parameters:

• treatment_bam (str): Path to file with data for treatment sample.
• control_bam (str): Path to file with data for control sample.
• treatment_name (str): Name for the treatment sample.
• control_name (str): Name for the control sample.
• output_dir (str): Path to folder for output.
• broad (str | bool): Whether to specify broad peak calling mode.
• cpus (int): Number of cores the script may use.
• qvalue (float): FDR, as decimal value (default: None)

Returns:

• str: Command to run.

topHat_map

def topHat_map(input_fastq: str, output_dir: str, genome: str, transcriptome: str, cpus: int | str) -> str

trimmomatic

def trimmomatic(input_fastq1: str, output_fastq1: str, cpus: int | str, adapters: str, log: str, input_fastq2: str | None = None, output_fastq1_unpaired: str | None = None, output_fastq2: str | None = None, output_fastq2_unpaired: str | None = None) -> str

validate_bam

def validate_bam(input_bam: str) -> str

Build Picard ValidateSamFile command.

zinba_call_peaks

def zinba_call_peaks(treatment_bed: str, control_bed: str, cpus: int, tagmented: bool = False) -> str