s = s.lower()

if s not in {'false', 'true', '0', '1'}:

raise ValueError('Not a valid boolean string')

parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter)

ref_args_group = parser.add_argument_group('Reference data')

input_args_group = parser.add_argument_group('Input data')

output_args_group = parser.add_argument_group('Output naming')

pipeline_args_group = parser.add_argument_group('Pipeline options')

algo_args_group = parser.add_argument_group('Algorithm settings')

output_setup_args_group = parser.add_argument_group('Output configuration')

align_args_group = parser.add_argument_group('Aligner settings')

filer_args_group = parser.add_argument_group('Read filtering options')

sc_args_group = parser.add_argument_group('Single-cell/spatial-related options:')

show_full_help = '--full_help' in cmd_args

def add_option_to_group(opt_group, *args, **kwargs): # show command only with --full-help

opt_group.add_argument(*args, **kwargs)

def add_additional_option_to_group(opt_group, *args, **kwargs): # show command only with --full-help

if not show_full_help:

kwargs['help'] = argparse.SUPPRESS

opt_group.add_argument(*args, **kwargs)

def add_hidden_option(*args, **kwargs): # show command only with --full-help

kwargs['help'] = argparse.SUPPRESS

parser.add_argument(*args, **kwargs)

parser.add_argument("--full_help", action='help', help="show full list of options")

parser.add_argument("--test", action=TestMode, nargs=0, help="run IsoQuant on toy dataset")

add_hidden_option('--debug', action='store_true', default=False,

help='debug log output.')

# REFERENCE

add_option_to_group(ref_args_group,"--reference", "-r",

help="reference genome in FASTA format (can be gzipped)",

type=str)

add_option_to_group(ref_args_group,"--genedb", "-g", help="gene database in gffutils DB format or GTF/GFF "

"format (optional)", type=str)

add_option_to_group(ref_args_group,'--complete_genedb', action='store_true', default=False,

help="use this flag if gene annotation contains transcript and gene features (e.g. for official annotations)")

add_additional_option_to_group(ref_args_group, "--discard_chr", nargs="+", help="chromosome IDs to ignore",

type=str, default=[])

add_additional_option_to_group(ref_args_group, "--process_only_chr", nargs="+", help="chromosome IDs to process",

type=str, default=None)

add_additional_option_to_group(ref_args_group, "--index", help="genome index for specified aligner (optional)",

type=str)

# INPUT READS

add_option_to_group(input_args_group, "--data_type", "-d", type=str, choices=DATA_TYPE_ALIASES.keys(),

help="type of data to process")

input_args = input_args_group.add_mutually_exclusive_group()

input_args.add_argument('--bam', nargs='+', type=str,

help='sorted and indexed BAM file(s), each file will be treated as a separate sample')

input_args.add_argument('--fastq', nargs='+', type=str,

help='input FASTQ/FASTA file(s) with reads, each file will be treated as a separate sample')

input_args.add_argument('--unmapped_bam', nargs='+', type=str,

help='unmapped BAM file(s), each file will be treated as a separate sample')

input_args.add_argument('--yaml', type=str, help='yaml file containing all input files, one entry per sample')

add_option_to_group(input_args_group, '--illumina_bam', nargs='+', type=str,

help='sorted and indexed file(s) with Illumina reads from the same sample')

add_option_to_group(input_args_group, "--read_group", nargs='+', type=str,

help="read grouping rules (space-separated); supported values: "

"tag:TAG (BAM tag), "

"file:FILE:READ_COLUMNS:GROUP_COLUMNS:DELIM (TSV file), "

"read_id:DELIM (read ID suffix), "

"file_name (filename), "

"barcode_spot (barcodes to spots/cell via --barcode2spot), "

"barcode_barcode (barcodes to spots via --barcode2barcode), "

"barcode (barcode--barcoded_reads), "

"no_auto (only explicitly stated rules), "

"none (disable grouping)")

add_additional_option_to_group(input_args_group, "--read_assignments", nargs='+', type=str,

help="reuse read assignments (binary format)", default=None)

# INPUT PROPERTIES

add_option_to_group(input_args_group, '--stranded', type=str, choices=SUPPORTED_STRANDEDNESS,

help="reads strandness type [none]", default="none")

add_option_to_group(input_args_group, '--polya_trimmed', default=PolyATrimmed.none.name, type=str,

help="define reads which had polyA tail trimmed or take polyA status from an external "

"source: none, all, stranded, tag:<TAG> (e.g. Dorado's tag:pt), "

"list:<FILE> (read ids, optional strand column), "

"flnc:<FILE> (isoseq refine flnc.report.csv) [%s]" % PolyATrimmed.none.name)

add_option_to_group(input_args_group, '--fl_data', action='store_true', default=False,

help="reads represent FL transcripts; both ends of the read are considered to be reliable")

# OUTPUT

add_option_to_group(output_args_group, "--output", "-o",

help="output folder, will be created automatically [isoquant_output]",

type=str, default="isoquant_output")

add_option_to_group(output_args_group,'--prefix', '-p', type=str,

help='experiment name; to be used for folder and file naming [OUT]',

default="OUT")

add_option_to_group(output_args_group,'--labels', '-l', nargs='+', type=str,

help='sample/replica labels to be used as column names; input file names are used '

'if not set; must be equal to the number of input files given via --fastq/--bam')

# PIPELINE STEPS

add_option_to_group(pipeline_args_group, "--analysis", nargs='+', type=str,

choices=ANALYSIS_CHOICES, default=None,

metavar="ANALYSIS",

help="analyses to run (space-separated); supported values: "

"quantification/quant, transcript_discovery/td, "

"exon_quantification/ex_quant, fusion [auto]")

add_option_to_group(pipeline_args_group, "--threads", "-t", help="number of threads [16]", type=int,

default="16")

# deprecated

add_additional_option_to_group(pipeline_args_group, "--count_exons",

help="deprecated: use --analysis exon_quantification",

action='store_true', default=False)

add_additional_option_to_group(pipeline_args_group, "--count_intron_retentions",

help="deprecated: use --analysis exon_quantification",

action='store_true', default=False)

add_additional_option_to_group(pipeline_args_group, "--no_model_construction", action="store_true",

default=False, help="deprecated: omit transcript_discovery from --analysis")

resume_args = pipeline_args_group.add_mutually_exclusive_group()

resume_args.add_argument("--resume", action="store_true", default=False,

help="resume failed run, specify output folder, input options are not allowed")

resume_args.add_argument("--force", action="store_true", default=False,

help="force to overwrite the previous run")

add_additional_option_to_group(pipeline_args_group, '--clean_start', action='store_true', default=False,

help='Do not use previously generated index, feature db or alignments.')

add_additional_option_to_group(pipeline_args_group, "--run_aligner_only", action="store_true", default=False,

help="align reads to reference without running further analysis")

add_additional_option_to_group(pipeline_args_group, "--no_gtf_check", help="do not perform GTF checks",

dest="gtf_check",

action='store_false', default=True)

add_additional_option_to_group(pipeline_args_group, "--high_memory",

help="increase RAM consumption (store alignment and the genome in RAM)",

action='store_true', default=False)

add_additional_option_to_group(pipeline_args_group, "--keep_tmp", help="do not remove temporary files "

"in the end", action='store_true',

default=False)

# OUTPUT SETUP

add_option_to_group(output_setup_args_group, '--check_canonical', action='store_true', default=False,

help="report whether splice junctions are canonical")

add_option_to_group(output_setup_args_group, "--sqanti_output", help="produce SQANTI-like TSV output",

action='store_true', default=False)

add_option_to_group(output_setup_args_group, "--bam_tags",

help="comma separated list of BAM tags to be imported to read_assignments.tsv",

type=str)

add_option_to_group(output_setup_args_group, "--large_output", nargs='+', type=str,

choices=LARGE_OUTPUT_TYPES,

default=["read_info", "read2transcripts"],

help="large output files to generate [read_info read2transcripts]")

add_additional_option_to_group(output_setup_args_group, "--no_report",

help="do not write the run summary (SAMPLE.summary.html and "

"SAMPLE.summary_stats.json)",

action='store_true', default=False)

add_additional_option_to_group(output_setup_args_group, "--genedb_output", help="output folder for converted gene "

"database, will be created automatically "

" (same as output by default)", type=str)

add_additional_option_to_group(output_setup_args_group, "--no_gzip", help="do not gzip large output files",

dest="gzipped", action='store_false', default=True)

add_additional_option_to_group(output_setup_args_group, "--normalization_method",

type=str, choices=[e.name for e in NormalizationMethod],

help="TPM normalization method",

default=NormalizationMethod.simple.name)

add_additional_option_to_group(output_setup_args_group, "--counts_format", type=str, nargs='+',

choices=[e.name for e in GroupedOutputFormat],

help="output format for grouped counts",

default=[GroupedOutputFormat.default.name])

add_additional_option_to_group(output_setup_args_group, "--emit_read_ids", action='store_true', default=False,

help="add read id columns to exon splice-site counts output")

add_additional_option_to_group(output_setup_args_group, "--old_exon_count_format", action='store_true', default=False,

help="output old exon inclusion/exclusion counts (deprecated)")

# ALIGNER

add_additional_option_to_group(align_args_group, "--aligner", choices=SUPPORTED_ALIGNERS,

help="use this alignmer method, can be [minimap2]", type=str)

add_additional_option_to_group(align_args_group, "--no_junc_bed", action="store_true", default=False,

help="do NOT use annotation for read mapping")

add_additional_option_to_group(align_args_group, "--junc_bed_file", type=str,

help="annotation in BED format produced by minimap's paftools.js gff2bed "

"(will be created automatically if not given)")

add_additional_option_to_group(align_args_group, "--indexing_options", type=str,

help="additional options that will be passed to the aligner indexer")

add_additional_option_to_group(align_args_group, "--mapping_options", type=str,

help="additional options that will be passed to the aligner")

# READ FILTERING

add_additional_option_to_group(filer_args_group, "--use_secondary",

help="use secondary alignments (slower processing)",

action='store_true', default=False)

add_additional_option_to_group(filer_args_group, "--min_mapq",

help="ignore alignments with MAPQ < this (including secondary alignments) [None]", type=int)

add_additional_option_to_group(filer_args_group, "--inconsistent_mapq_cutoff",

help="ignore inconsistent alignments with MAPQ < cutoff (works only with the reference annotation) [5]",

type=int, default=5)

add_additional_option_to_group(filer_args_group, "--simple_alignments_mapq_cutoff",

help="ignore alignments with 1 or 2 exons and MAPQ < cutoff "

"(works only in annotation-free mode) [1]", type=int, default=1)

add_additional_option_to_group(filer_args_group, "--max_coverage_small_chr",

help="process only a fraction of reads for high-coverage loci on small chromosomes (e.g. MT), "

"improves running time and RAM [1000000]",

type=int, default=1000000)

add_additional_option_to_group(filer_args_group, "--max_coverage_normal_chr",

help="process only a fraction of reads for high-coverage loci on usual chromosomes, "

"improves running time and RAM [-1]",

type=int, default=-1)

# SC ARGUMENTS

add_option_to_group(sc_args_group, "--mode", "-m", type=str,

choices=ISOQUANT_MODES + list(DEPRECATED_MODE_ALIASES),

help="IsoQuant mode [%s]" % IsoQuantMode.bulk.name, default=IsoQuantMode.bulk.name,

metavar="{%s}" % ",".join(ISOQUANT_MODES))

add_option_to_group(sc_args_group, "--split_molecules", type=str, choices=SPLIT_MOLECULES_CHOICES,

default=None,

help="split reads containing several cDNA molecules ligated end to end: "

"%s splits wherever the protocol allows it, %s also fails for protocols "

"that cannot [%s]" % (SPLIT_MOLECULES_AUTO, SPLIT_MOLECULES_TRUE,

SPLIT_MOLECULES_AUTO))

add_option_to_group(sc_args_group, '--barcode_whitelist', type=str, nargs='+',

help='file(s) with barcode whitelist(s) for barcode calling, '

'or "%s" to detect cell barcodes from the data' % AUTO_BARCODES)

add_option_to_group(sc_args_group, "--barcoded_reads", type=str, nargs='+',

help='TSV file(s) with barcoded reads')

add_option_to_group(sc_args_group, "--barcoded_bam", action='store_true', default=False,

help='extract barcodes and UMIs from BAM tags (CB/UB by default)')

add_option_to_group(sc_args_group, "--barcode2spot", type=str,

help='TSV file mapping barcode to cell type / spot id.')

add_option_to_group(sc_args_group, "--molecule", type=str,

help='molecule definition file (MDF) for custom_sc mode: '

'defines molecule structure for universal barcode extraction')

add_option_to_group(sc_args_group, "--n_cells", type=str,

help='expected number of cell-associated barcodes, or "%s" to estimate it '

'from the barcode count distribution. When set, the whitelist is '

'treated as a pool to select cell barcodes from; when omitted, the '

'whitelist is taken to be the cell barcodes themselves' % AUTO_BARCODES)

add_additional_option_to_group(sc_args_group, "--n_cells_interval", type=int, default=25,

help='percentage by which the number of selected barcodes may differ '

'from --n_cells [25]')

add_additional_option_to_group(sc_args_group, "--barcode_correction", type=str,

choices=[e.name for e in BarcodeCorrectionMethod],

default=BarcodeCorrectionMethod.auto.name,

help="override how the cell barcode list is obtained: %s matches reads "

"against the whitelist as given, %s selects cell barcodes from read "

"counts first [%s: decided by --n_cells]"

% (BarcodeCorrectionMethod.whitelist.name,

BarcodeCorrectionMethod.detect.name,

BarcodeCorrectionMethod.auto.name))

add_additional_option_to_group(sc_args_group, "--barcode_tag", type=str, default="CB",

help='BAM tag for cell barcode [CB]')

add_additional_option_to_group(sc_args_group, "--umi_tag", type=str, default="UB",

help='BAM tag for UMI [UB]')

add_additional_option_to_group(sc_args_group, "--strip_barcode_suffix", action='store_true', default=False,

help='remove suffix after dash from barcodes extracted from BAM tag')

add_additional_option_to_group(sc_args_group, "--barcode2barcode", type=str,

help='TSV file mapping barcode to spot IDs for UMI deduplication')

# ALGORITHM

add_additional_option_to_group(algo_args_group, "--report_novel_unspliced", "-u", type=bool_str,

help="report novel monoexonic transcripts (true/false), "

"default: false for ONT, true for other data types")

add_additional_option_to_group(algo_args_group, "--report_canonical", type=str,

choices=[e.name for e in StrandnessReportingLevel],

help="reporting level for novel transcripts based on canonical splice sites "

"[%s]" % StrandnessReportingLevel.auto.name,

default=StrandnessReportingLevel.auto.name)

add_additional_option_to_group(algo_args_group, "--polya_requirement", type=str,

choices=[e.name for e in PolyAUsageStrategies],

help="require polyA tails to be present when reporting transcripts; "

"default: auto (requires polyA only when polyA percentage is >= 70%%)",

default=PolyAUsageStrategies.auto.name)

# Alternative-polyA/TSS isoform discovery is default ON when an annotation is

# given (--genedb). --novel_apa (default off) extends it to novel

# (non-reference) transcripts, not only known ones.

add_additional_option_to_group(algo_args_group, "--novel_apa", action="store_true", default=False,

help="generate novel isoforms with identical intron chain but distinct polyA sites")

# Splice-site correction shifts novel transcript-model junctions onto canonical

# motifs using clustered read deletions near junctions; default ON.

add_additional_option_to_group(algo_args_group, "--no_splice_site_correction", action="store_true",

default=False,

help="disable canonical splice-site correction for novel transcript models")

add_additional_option_to_group(algo_args_group, "--transcript_quantification", choices=COUNTING_STRATEGIES,

help="transcript quantification strategy [%s]" % CountingStrategy.unique_only.name,

type=str, default=CountingStrategy.unique_only.name)

add_additional_option_to_group(algo_args_group, "--gene_quantification", choices=COUNTING_STRATEGIES,

help="gene quantification strategy [%s]" % CountingStrategy.unique_splicing_consistent.name,

type=str, default=CountingStrategy.unique_splicing_consistent.name)

add_additional_option_to_group(algo_args_group, "--matching_strategy",

choices=["exact", "precise", "default", "loose"],

help="read-to-isoform matching strategy from the most strict to least",

type=str, default=None)

add_additional_option_to_group(algo_args_group, "--splice_correction_strategy",

choices=["none", "default_pacbio", "default_ont",

"conservative_ont", "all", "assembly"],

help="read alignment correction strategy to use", type=str, default=None)

add_additional_option_to_group(algo_args_group, "--model_construction_strategy",

choices=["reliable", "default_pacbio", "sensitive_pacbio", "pacbio_all",

"fl_pacbio", "default_ont", "sensitive_ont", "all", "assembly"],

help="transcript model construction strategy to use", type=str, default=None)

add_additional_option_to_group(algo_args_group, "--delta", type=int, default=None,

help="delta for inexact splice junction comparison [auto]")

add_additional_option_to_group(algo_args_group, "--use_replicas", type=bool_str, default=True,

help="require novel transcripts to be confirmed by multiple files "

"when file_name grouping is used [true]")

# REST

add_hidden_option("--graph_clustering_distance", type=int, default=None,

help="intron graph clustering distance, "

"splice junctions less that this number of bp apart will not be differentiated")

add_hidden_option("--cage", help="bed file with CAGE peaks", type=str, default=None)

add_hidden_option("--cage-shift", type=int, default=50, help="interval before read start to look for CAGE peak")

# PolyA / TSS training-data collection (developer-only).

# When either flag is given to a normal IsoQuant run, the matching

# terminal-position counter switches from inference to dumping per-peak

# features (the eight FEATURE_COLUMNS plus `chromosome` and a `true_peak`

# 0/1 label computed against the annotated transcript end) to the CSV

# path supplied as the flag value. The XGBoost model is not consulted;

# rows are accumulated, not filtered. Feed the CSV to

# `misc/train_polya_tss_model.py` to fit a fresh model.

# python isoquant.py … --genedb GENEDB --collect_polya_training peaks.csv

# python misc/train_polya_tss_model.py --features peaks.csv \

# --output isoquant_lib/data/model_polya.json

# For TSS, add --fl_data and --collect_tss_training tss_peaks.csv.

# These flags are intentionally hidden from --help/--full_help; the run

# emits a clearly marked warning when they are used.

add_hidden_option("--collect_polya_training", type=str, default=None,

help="developer: dump per-peak features + true_peak label for polyA training to this CSV path.")

add_hidden_option("--collect_tss_training", type=str, default=None,

help="developer: dump per-peak features + true_peak label for TSS training to this CSV path.")

isoquant_version = "4.1.0"

try:

with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), "VERSION")) as version_f:

isoquant_version = version_f.readline().strip()

except FileNotFoundError:

try:

from importlib.metadata import version as _get_version

isoquant_version = _get_version("isoquant")

except Exception:

pass

parser.add_argument('--version', '-v', action='version', version='IsoQuant ' + isoquant_version,

help="show IsoQuant version and exit")

args = parser.parse_args(cmd_args, namespace)

if args.resume:

resume_parser = argparse.ArgumentParser(add_help=False)

resume_parser.add_argument("--resume", action="store_true", default=False,

help="resume failed run, specify only output folder, "

"input options are not allowed")

resume_parser.add_argument("--output", "-o",

help="output folder, will be created automatically [isoquant_output]",

type=str, required=True)

resume_parser.add_argument('--debug', action='store_true', default=argparse.SUPPRESS,

help='debug log output.')

resume_parser.add_argument("--threads", "-t", help="number of threads to use",

type=int, default=argparse.SUPPRESS)

resume_parser.add_argument("--high_memory",

help="increase RAM consumption (store alignment and the genome in RAM)",

action='store_true', default=False)

resume_parser.add_argument("--keep_tmp", help="do not remove temporary files in the end",

action='store_true', default=argparse.SUPPRESS)

args, unknown_args = resume_parser.parse_known_args(cmd_args)

if unknown_args:

logger.error("You cannot specify options other than --output/--threads/--debug/--high_memory/--keep_tmp "

"with --resume option")

parser.print_usage()

sys.exit(IsoQuantExitCode.INCOMPATIBLE_OPTIONS)

args._cmd_line = " ".join(sys.argv)

args._version = isoquant_version

args.output_exists = os.path.exists(args.output)

if not args.output_exists:

os.makedirs(args.output)

"""Extract all BAM file paths from input_data.samples.

bam_files: list[str] = []

for sample in input_data.samples:

for lib in sample.file_list:

for in_file in lib:

bam_files.append(in_file)

if getattr(sample, "illumina_bam", None):

bam_files.extend(sample.illumina_bam)

"""Run fusion detection per sample using a shared FusionDetector.

summary = {"total": 0, "successful": 0, "failed": 0, "skipped": []}

for sample in samples:

sample_bams = fusion_input_bams(sample)

if not sample_bams:

continue

fusion_checkpoint = checkpoint_name("fusion", sample.prefix)

if checkpoints is not None and checkpoints.is_done(fusion_checkpoint):

logger.info("Fusion detection for sample %s was done during the previous run, skipping" % sample.prefix)

summary["total"] += 1

summary["successful"] += 1

continue

summary["total"] += 1

out_fname = os.path.join(sample.out_dir, sample.prefix + ".fusions.tsv")

try:

fd.clear_state() # Accumulate all of this sample's BAMs into one report

for bam_path in sample_bams:

logger.info("Running fusion detection on %s" % bam_path)

fd.bam_path = bam_path # Switch BAM

fd.detect_fusions()

fd.report(output_path=out_fname)

logger.info("Fusion candidates for sample %s written to %s" % (sample.prefix, out_fname))

except Exception as e:

logger.error("Fusion detection failed for sample %s: %s" % (sample.prefix, str(e)))

logger.debug("Traceback:", exc_info=True)

summary["failed"] += 1

summary["skipped"].append(sample.prefix)

continue

summary["successful"] += 1

# only a success is marked, outside the try: a failed sample stays retryable on resume,

# and a failure to write the checkpoint is not reported as a failed detection

if checkpoints is not None:

mark_done_with_debug_hook(checkpoints, fusion_checkpoint)

"""True if some sample still needs fusion detection (building the detector is expensive).

return any(fusion_input_bams(sample) and not checkpoints.is_done(checkpoint_name("fusion", sample.prefix))

args.param_file = os.path.join(args.output, ".params")

if args.resume:

if not os.path.exists(args.output) or not os.path.exists(args.param_file):

# logger is not defined yet

logger.error("Previous run config was not detected, cannot resume. "

"Check that output folder is correctly specified.")

sys.exit(IsoQuantExitCode.RESUME_CONFIG_NOT_FOUND)

args = load_previous_run(args)

# every run of this version writes a checkpoint record and saves its id in .params

reason = get_resume_incompatibility_reason(run_checkpoint_dir(args.output), getattr(args, "checkpoint_run_id", None),

args._version)

if reason:

logger.error("The run in %s cannot be resumed: %s. Restart it without --resume." % (args.output, reason))

sys.exit(IsoQuantExitCode.RESUME_INCOMPATIBLE)

elif args.output_exists:

if os.path.exists(args.param_file):

if args.force:

logger.warning("Output folder already contains a previous run, will be overwritten.")

else:

logger.warning("Output folder already contains a previous run, some files may be overwritten. "

"Use --resume to resume a failed run. Use --force to avoid this message.")

logger.warning("Press Ctrl+C to interrupt the run now.")

delay = 9

for i in range(delay):

countdown = delay - i

sys.stdout.write("Resuming the run in %d second%s\r" % (countdown, "s" if countdown > 1 else ""))

time.sleep(1)

logger.info("Overwriting the previous run")

time.sleep(1)

else:

logger.warning("Output folder already exists, some files may be overwritten.")

if args.genedb_output is None:

args.genedb_output = args.output

elif not os.path.exists(args.genedb_output):

os.makedirs(args.genedb_output)

if not args.genedb:

args.genedb_filename = None

elif args.genedb.lower().endswith("db"):

args.genedb_filename = args.genedb

else:

args.genedb_filename = os.path.join(args.genedb_output, os.path.splitext(os.path.basename(args.genedb))[0] + ".db")

if not check_input_params(args):

parser.print_usage()

sys.exit(IsoQuantExitCode.INVALID_PARAMETER)

# Validate --read_group none

if args.read_group:

if "none" in args.read_group and len(args.read_group) > 1:

logger.error("--read_group 'none' cannot be combined with other values")

sys.exit(IsoQuantExitCode.INVALID_PARAMETER)

# Validate --large_output values

if args.large_output:

if "none" in args.large_output and len(args.large_output) > 1:

logger.error("--large_output 'none' cannot be combined with other values")

sys.exit(IsoQuantExitCode.INVALID_PARAMETER)

for val in args.large_output:

if val not in LARGE_OUTPUT_TYPES:

logger.error("Invalid --large_output value: %s. Valid values: %s" % (val, ", ".join(LARGE_OUTPUT_TYPES)))

sys.exit(IsoQuantExitCode.INVALID_PARAMETER)

_warn_about_unusable_bam_outputs(args)

if not args.resume:

reset_checkpoints(args)

save_params(args)

"""Start a fresh run: drop every checkpoint a previous run left in this output folder.

for sample in args.input_data.samples:

shutil.rmtree(sample.checkpoint_dir, ignore_errors=True)

run_dir = run_checkpoint_dir(args.output)

CheckpointStore(run_dir, resume=False).reset()

# saved in .params below; --resume checks it against the record

"""Say up front when a requested BAM output cannot be produced, rather than mid-run."""

if not args.mode.needs_pcr_deduplication():

for output_type in ("tagged_bam", "deduplicated_bam"):

if output_type in args.large_output:

logger.warning("--large_output %s has no effect in %s mode, which has no barcodes "

"or UMIs; it will be skipped" % (output_type, args.mode.name))

elif "tagged_bam" in args.large_output and getattr(args, 'barcoded_bam', False):

logger.warning("--large_output tagged_bam is redundant with --barcoded_bam: those "

logger.info("Loading parameters from the previous run")

logger.info("Only --output/--threads/--debug/--high_memory/--keep_tmp are compatible with --resume option")

unpickler = pickle.Unpickler(open(args.param_file, "rb"), fix_imports=False)

loaded_args = unpickler.load()

for option in args.__dict__:

loaded_args.__dict__[option] = args.__dict__[option]

if loaded_args.debug:

logger.setLevel(logging.DEBUG)

logger.handlers[0].setLevel(logging.DEBUG)

for file_opt in ["genedb", "reference", "index", "bam", "fastq", "junc_bed_file",

"cage", "genedb_output", "read_assignments"]:

if file_opt in args.__dict__ and args.__dict__[file_opt]:

if isinstance(args.__dict__[file_opt], list):

args.__dict__[file_opt] = list(map(os.path.abspath, args.__dict__[file_opt]))

else:

args.__dict__[file_opt] = os.path.abspath(args.__dict__[file_opt])

if "read_group" in args.__dict__ and args.__dict__["read_group"]:

updated_specs = []

for spec in args.read_group:

vals = spec.split(":")

if len(vals) > 1 and vals[0] == 'file':

vals[1] = os.path.abspath(vals[1])

updated_specs.append(":".join(vals))

else:

updated_specs.append(spec)

args.read_group = updated_specs

polya_name, sep, polya_file = args.polya_trimmed.partition(":")

if sep and polya_name in (PolyATrimmed.list.name, PolyATrimmed.flnc.name) and polya_file:

args.polya_trimmed = polya_name + ":" + os.path.abspath(polya_file)

pickler = pickle.Pickler(open(args.param_file, "wb"), -1)

pickler.dump(args)

is_single_cell = args.mode.needs_barcode_calling()

# capture legacy stage flags before they are overwritten below

legacy_count_exons = args.count_exons

legacy_intron_retentions = args.count_intron_retentions

legacy_no_model = args.no_model_construction

if args.analysis is not None:

analyses = {ANALYSIS_ALIASES[a] for a in args.analysis}

elif is_single_cell:

analyses = {AnalysisType.quantification}

elif args.genedb:

analyses = {AnalysisType.quantification, AnalysisType.transcript_discovery}

else:

analyses = {AnalysisType.transcript_discovery}

# apply deprecated stage flags additively, on top of the chosen analyses

if legacy_count_exons:

logger.warning("--count_exons is deprecated, use --analysis exon_quantification")

analyses.add(AnalysisType.exon_quantification)

if legacy_intron_retentions:

logger.warning("--count_intron_retentions is deprecated, use --analysis exon_quantification")

if legacy_no_model:

logger.warning("--no_model_construction is deprecated, omit transcript_discovery from --analysis")

analyses.discard(AnalysisType.transcript_discovery)

# feasibility: annotation-dependent analyses cannot run without a gene database

if not args.genedb:

for analysis in (AnalysisType.quantification, AnalysisType.exon_quantification, AnalysisType.fusion):

if analysis in analyses:

logger.warning("%s requires a reference annotation (--genedb) and will not be run" % analysis.name)

analyses.discard(analysis)

# without an annotation transcript discovery is the only available analysis

if not analyses and not legacy_no_model:

logger.info("Without a reference annotation only transcript discovery is available; "

"running transcript_discovery")

analyses.add(AnalysisType.transcript_discovery)

# feasibility: model construction in single-cell/spatial modes is limited

if is_single_cell and AnalysisType.transcript_discovery in analyses:

logger.warning("transcript_discovery in single-cell/spatial mode runs after UMI deduplication: "

"it will not yield novel genes and may be incomplete; "

"consider a pseudo-bulk run for transcript discovery")

# canonical internal flags (names reused across the pipeline)

args.run_quantification = AnalysisType.quantification in analyses

args.count_exons = AnalysisType.exon_quantification in analyses or legacy_count_exons

args.count_intron_retentions = AnalysisType.exon_quantification in analyses or legacy_intron_retentions

args.fusion = AnalysisType.fusion in analyses

args.no_model_construction = AnalysisType.transcript_discovery not in analyses

# polyA/TSS site prediction is part of quantification and is also needed for

# model construction; gated by the annotation (TSS additionally needs --fl_data)

args.predict_terminal_sites = bool(args.genedb) and (args.run_quantification or not args.no_model_construction)

if not args.reference:

logger.error("Reference genome was not provided")

return False

if not args.data_type:

logger.error("Data type is not provided, choose one of " + " ".join(DATA_TYPE_ALIASES.keys()))

return False

elif args.data_type not in DATA_TYPE_ALIASES.keys():

logger.error("Unsupported data type " + args.data_type + ", choose one of: " + " ".join(DATA_TYPE_ALIASES.keys()))

return False

args.data_type = DATA_TYPE_ALIASES[args.data_type]

if not any([args.fastq, args.bam, args.unmapped_bam, args.read_assignments, args.yaml]):

logger.error("No input data was provided")

return False

if args.yaml and args.illumina_bam:

logger.error("When providing a yaml file it should include all input files, including the illumina bam file.")

return False

args.input_data = InputDataStorage(args)

# the run-level checkpoint store is <output>/checkpoints and is wiped on a fresh run

run_store_dir = os.path.normpath(run_checkpoint_dir(args.output))

for sample in args.input_data.samples:

if os.path.normpath(sample.out_dir) == run_store_dir:

logger.error("Experiment name '%s' is reserved, please choose a different prefix or name"

% os.path.basename(run_store_dir))

return False

if args.aligner is not None and args.aligner not in SUPPORTED_ALIGNERS:

logger.error(" Unsupported aligner " + args.aligner + ", choose one of: " + " ".join(SUPPORTED_ALIGNERS))

return False

if args.run_aligner_only and not args.input_data.input_type.needs_mapping():

logger.error("Data type %s cannot be mapped and thus incompatible with --run_aligner_only option." % args.input_data.input_type.name)

return False

if args.stranded not in SUPPORTED_STRANDEDNESS:

logger.error("Unsupported strandness " + args.stranded + ", choose one of: " + " ".join(SUPPORTED_STRANDEDNESS))

return False

if not args.genedb and args.sqanti_output:

args.sqanti_output = False

logger.warning("--sqanti_output option has no effect without gene annotation")

if args.no_model_construction and args.sqanti_output:

args.sqanti_output = False

logger.warning("--sqanti_output option has no effect without model construction")

if args.process_only_chr and args.discard_chr:

args.discard_chr = []

if "read_group" in args.__dict__ and args.__dict__["read_group"]:

updated_specs = []

spec_set = set()

for spec in args.read_group:

if spec in spec_set:

logger.warning("Read group %s is set twice, which has no effect, duplicate will be ignored" % spec)

continue

updated_specs.append(spec)

spec_set.add(spec)

"""Refuse to split molecules when the reads are already aligned.

if not args.split_molecules or args.input_data.input_type.needs_mapping():

return

logger.critical("Reads cannot be split into separate molecules when they are already "

"aligned (%s input): the split molecules would have to be mapped again. "

"Provide the raw reads instead, or use --split_molecules %s to analyse "

"the alignments as they are."

% (args.input_data.input_type.name, SPLIT_MOLECULES_FALSE))

if not _validate_data_type_and_input(args):

return False

if not _reject_reserved_sample_names(args):

return False

if not _validate_alignment_options(args):

return False

if not isinstance(args.mode, IsoQuantMode):

resolve_deprecated_mode(args)

args.mode = IsoQuantMode[args.mode]

resolve_split_molecules(args)

_reject_splitting_aligned_input(args)

# translate --analysis (and the deprecated stage flags) into internal booleans

resolve_analyses(args)

_apply_stage_warnings(args)

_dedup_read_group_specs(args)

validate_barcode_calling(args)

check_input_files(args)

"""Check BAM file exists and optionally has index."""

if not os.path.isfile(bam_path):

logger.critical("BAM file " + bam_path + " does not exist")

sys.exit(IsoQuantExitCode.INPUT_FILE_NOT_FOUND)

if check_index:

bamfile_in = pysam.AlignmentFile(bam_path, "rb")

if not bamfile_in.has_index():

logger.critical("BAM file " + bam_path + " is not indexed, run samtools sort and samtools index")

sys.exit(IsoQuantExitCode.BAM_NOT_INDEXED)

"""

parts = spec.split(":")

if len(parts) >= 2 and parts[0] == 'file':

return parts[1]

"""Validate a single input read file (save prefix / BAM / FASTQ)."""

if args.input_data.input_type == InputDataType.save:

saves = glob.glob(in_file + "*")

if not saves:

logger.critical("Input files " + in_file + "* do not exist")

return

if not os.path.isfile(in_file):

logger.critical("Input file " + in_file + " does not exist")

sys.exit(IsoQuantExitCode.INPUT_FILE_NOT_FOUND)

if args.input_data.input_type == InputDataType.bam:

# Check reference genome

if args.reference and not os.path.isfile(args.reference):

logger.critical("Reference genome " + args.reference + " does not exist")

sys.exit(IsoQuantExitCode.INPUT_FILE_NOT_FOUND)

# Check input reads (BAM/FASTQ/save files)

for sample in args.input_data.samples:

for lib in sample.file_list:

for in_file in lib:

_check_input_read_file(args, in_file)

# Check Illumina BAM files

if sample.illumina_bam is not None:

for illumina in sample.illumina_bam:

# Check read_group file specs

if hasattr(args, 'read_group') and args.read_group:

for spec in args.read_group:

file_path = extract_read_group_file_path(spec)

if file_path:

# Check junction BED file

if hasattr(args, 'junc_bed_file') and args.junc_bed_file:

check_file_exists(args.junc_bed_file, "Junction BED file")

# Check CAGE file (currently not supported)

if args.cage is not None:

logger.critical("CAGE data is not supported yet")

sys.exit(IsoQuantExitCode.INVALID_PARAMETER)

if not os.path.isfile(args.cage):

logger.critical("Bed file with CAGE peaks " + args.cage + " does not exist")

sys.exit(IsoQuantExitCode.INPUT_FILE_NOT_FOUND)

# Check gene database

if args.genedb is not None:

if not os.path.isfile(args.genedb):

logger.critical("Gene database " + args.genedb + " does not exist")

sys.exit(IsoQuantExitCode.GENE_DB_NOT_FOUND)

else:

args.no_junc_bed = True

# Check read assignments

if args.read_assignments is not None:

for r in args.read_assignments:

if not glob.glob(r + "*"):

logger.critical("No files found with prefix " + str(r))

_check_reference_and_reads(args)

check_barcode_input_files(args)

check_barcode_mapping_files(args)

_check_read_group_files(args)

for sample in args.input_data.samples:

sample_dir = sample.out_dir

if os.path.exists(sample_dir):

if not args.resume:

logger.warning(sample_dir + " folder already exists, some files may be overwritten")

else: