stitch-reference-transcripts

Stitch reference transcript sequence onto the degraded ends of aligned RNA sequences. The stitched sequences no longer match the read coordinates of any call-rna-transcript-vars output computed on the same BAM file: realign the stitched FASTA and re-run call-rna-transcript-vars before translate-transcripts.

Usage

Template:

exacto stitch-reference-transcripts \
    --bam-file <bam_file> \
    --reference-genome-fasta-file <reference_genome_fasta_file> \
    --reference-gene-annotation-file <reference_gene_annotation_file> \
    --reference-gene-annotation-source <reference_gene_annotation_source> \
    --reference-gene-annotation-assembly <reference_gene_annotation_assembly> \
    --reference-gene-annotation-version <reference_gene_annotation_version> \
    --output-dir <output_dir> \
    --output-prefix <output_prefix> \
    [--min-mapping-quality MIN_MAPPING_QUALITY] \
    [--min-num-splice-junction-matches MIN_NUM_SPLICE_JUNCTION_MATCHES] \
    [--pas-search-size PAS_SEARCH_SIZE] \
    [--pas-hexamers PAS_HEXAMERS [PAS_HEXAMERS ...]] \
    [--pas-start-offset-range PAS_START_OFFSET_RANGE] \
    [--polya-tail-min-adenosine-fraction POLYA_TAIL_MIN_ADENOSINE_FRACTION] \
    [--polya-window-size POLYA_WINDOW_SIZE] \
    [--polya-min-consecutive-adenosine POLYA_MIN_CONSECUTIVE_ADENOSINE] \
    [--num-threads NUM_THREADS]

Example:

exacto stitch-reference-transcripts \
    --bam-file pass2/tumor_pass2_exacto_rna_consensus.sorted.bam \
    --reference-genome-fasta-file reference_genome.fasta \
    --reference-gene-annotation-file gencode.gtf.gz \
    --reference-gene-annotation-source gencode \
    --reference-gene-annotation-assembly hg38 \
    --reference-gene-annotation-version v45 \
    --output-dir stitched/ \
    --output-prefix tumor

Description

Stitch reference transcript sequence onto the degraded ends of aligned RNA sequences. The stitched sequences no longer match the read coordinates of any call-rna-transcript-vars output computed on the same BAM file: realign the stitched FASTA and re-run call-rna-transcript-vars before translate-transcripts.

NoteAt a glance

Inputs: *.bam (aligned RNA sequences with cs tags — the realigned determine-rna-consensus FASTA), *.fasta (reference genome), *.gtf.gz (gene annotation)

Outputs: {prefix}_exacto_stitched_transcripts.tsv, {prefix}_exacto_stitched_transcripts.fasta (one record per transcript, named as in the BAM)

Typical next step: Realign the stitched FASTA with minimap2, run call-rna-transcript-vars on it, then pass the stitched TSV to translate-transcripts — the stitched sequences no longer match the pre-stitch coordinates

Required arguments

Flag	Type	Description
--bam-file	str	Input BAM file of the aligned RNA sequences to stitch (e.g. the realigned consensus sequences of exacto determine-rna-consensus). Every mapped read is one transcript, identified by its read name.
--reference-genome-fasta-file	str	Reference genome FASTA file.
--reference-gene-annotation-file	str	Reference gene annotation file.
--reference-gene-annotation-source	str	Reference gene annotation source.
--reference-gene-annotation-assembly	str	Reference gene annotation assembly (e.g. ‘hg38’).
--reference-gene-annotation-version	str	Reference gene annotation version (e.g. ‘v41’).
--output-dir	str	Output directory.
--output-prefix	str	Output prefix.

Optional arguments

Flag	Type	Default	Description
--min-mapping-quality	int	0	Minimum mapping quality. A read whose best-mapped record falls below this is left out of the output.
--min-num-splice-junction-matches	int	1	A reference transcript match with fewer splice junction matches than this is treated as no match. A single-exon reference transcript is exempt.
--pas-search-size	int	40	Number of read bases, ending at an intronic 3’ terminal aligned base, searched for a polyadenylation signal hexamer.
--pas-hexamers	str	AATAAA, ATTAAA	Polyadenylation signal hexamers, in transcript orientation. Pass multiple as: –pas-hexamers AATAAA ATTAAA.
--pas-start-offset-range	int	(10, 40)	Permitted distance in bases, inclusive, from a polyadenylation signal hexamer’s first base to the end of the searched read sequence. Pass as: –pas-start-offset-range 10 40.
--polya-tail-min-adenosine-fraction	float	0.8	Read bases past the 3’ terminal aligned base are accepted as a polyA tail when at least this fraction of them is A; otherwise they may be the remnant of a fusion or aberrant splicing, and the 3’ end is left alone.
--polya-window-size	int	20	Number of genomic bases on each side of an intronic 3’ terminal aligned base scanned for an A-tract.
--polya-min-consecutive-adenosine	int	6	An A-run of at least this length in that genomic window calls internal priming. An intronic 3’ end with a polyadenylation signal and no such A-run is kept as a genuine intronic polyadenylation site and is not stitched.
--num-threads	int	4	Number of threads.