bp-quant 0.5.4

Description:

bpquant 0.5.4

bp-quant
Quantification of reads at defined positions to verify custom input sequences.
Given a gene fusion or splicing junction of interest, this tool can quantify
RNA-seq reads supporting the breakpoint (or splice junction) by quantifying
reads that map to the breakpoint (junction reads) and read pairs that span the
breakpoint (spanning pairs).
Workflow

Input:

Table with target sequences and breakpoints position (CSV/TSV format)
fastq files / BAM file (BAM input only works in combination with STAR as aligner)


Convert target sequences to FASTA format (bp_quant csv2fasta)
Map reads against sequences using STAR/Bowtie2/BWA

Generate Index of sequences as reference (bp_quant index)
Map reads (bp_quant align)


Count reads using bp_quant count
Output:

Table with read counts per input sequence
Table with info on each read (read_info.tsv)



Dependencies
Look at this environment.yml for the individual dependencies (Hint: File can also be used to build conda environment).
Installation
# Create conda environment from above or add the dependencies to your path
pip install bp-quant

Alternatively, bp-quant can be installed from github repository:
git clone https://github.com/TRON-Bioinformatics/easyquant.git

cd easyquant

# If you have conda installed you can simply install the environment like this
conda env create -f environment.yml --prefix conda_env/
conda activate conda_env/

python -m build
pip install dist/*.whl

Usage
usage: bp_quant pipeline [-h] [-1 FQ1] [-2 FQ2] [-b BAM] -s SEQ_TAB -o
OUTPUT_FOLDER [-d BP_DISTANCE] [--allow_mismatches]
[--interval_mode] [-m {star,bowtie2,bwa}]
[-t NUM_THREADS] [--star_cmd_params STAR_CMD_PARAMS]
[--keep_aln | --keep_all]

Runs the complete bpquant pipeline

optional arguments:
-h, --help show this help message and exit
-1 FQ1, --fq1 FQ1 Specify path to Read 1 (R1) FASTQ file
-2 FQ2, --fq2 FQ2 Specify path to Read 2 (R2) FASTQ file
-b BAM, --bam_file BAM
Specify path to input BAM file as alternative to FASTQ
input
-s SEQ_TAB, --sequence_tab SEQ_TAB
Specify the reference sequences as table with colums
name, sequence, and position
-o OUTPUT_FOLDER, --output_folder OUTPUT_FOLDER
Specify the folder to save the results into.
-d BP_DISTANCE, --bp_distance BP_DISTANCE
Threshold in base pairs for the required overlap size
of reads on both sides of the breakpoint for
junction/spanning read counting
--allow_mismatches Allow mismatches within the region around the
breakpoint determined by the bp_distance parameter
--interval_mode Specify if interval mode shall be used
-m {star,bowtie2,bwa}, --method {star,bowtie2,bwa}
Specify alignment software to generate the index
-t NUM_THREADS, --threads NUM_THREADS
Specify number of threads to use for the alignment
--star_cmd_params STAR_CMD_PARAMS
Specify STAR commandline parameters to use for the
alignment
--keep_aln Do not delete alignment files during clean up step
--keep_all Do not perform clean up step after re-quantification

Copyright (c) 2023 TRON gGmbH (See LICENSE for licensing details)

Use case with example data
We use toy example data from the folder example_data. It consists of a table
with input sequences and positions, as well as two fastq files / one BAM file.
Fastqs as input:
bp_quant pipeline \
-1 example_data/example_rna-seq_R1_001.fastq.gz \
-2 example_data/example_rna-seq_R1_001.fastq.gz \
-s example_data/CLDN18_Context_seq.csv \
-d 10 \
-o example_out \
-m star \
-t 6
[--interval_mode]

BAM as input:
bp_quant pipeline \
-b example_data/example_rna-seq.bam \
-s example_data/CLDN18_Context_seq.csv \
-d 10 \
-o example_out \
-m star \
-t 6
[--interval_mode]

Input
Table with input sequences
As input a CSV/TSV table should be given holding the target sequence
with unique names and the relative position(s) of the breakpoint(s)/intervals or fusion junction.
Example format of the input table:



name
sequence
position




seq1
AACCGCCACCG
5


seq2
GTCCGTTGGCG
5


seq3
AACCGCCCTGT
5


seq4
CGGCATCATCG
0,5,10



Fastq files / BAM file
Paired fastq files or an unsorted BAM file (no multimappers) as input is required
to successfully determine read classes as described below.
Config file
Output format
The main output consists of the file:

<OUTPUT_FOLDER>/quantification.tsv contains raw read counts for each sequence

The output of the example data <OUTPUT_FOLDER>/quantification.tsv using a mismatch ratio of 0.05 (default) should look like this:



name
pos
junc
span
anch
a
b




CLDN18_1
400
570
684
25
1109
3932


CLDN18_2
361
0
1
0
1
2968


CLDN18_total
400
593
688
25
4315
4990


CLDN18_1_fake
400
0
3
0
3
3946


CLDN18_2_fake
361
0
0
0
0
3943


HPRT1
400
107
215
25
1259
974



Using the interval mode the output will look slightly different:



name
interval
overlap_interval_end_reads
span_interval_end_pairs
within_interval
coverage_perc
coverage_mean
coverage_median




CLDN18_1
0_400
570
684
1109
0.89
182.71
137.5


CLDN18_1
400_786
0
0
3932
1.0
519.51
563.5


CLDN18_2
0_361
0
1
1
0.14
0.14
0.0


CLDN18_2
361_747
0
0
2968
1.0
392.15
425.5


CLDN18_total
0_400
593
688
4315
1.0
586.16
682.0


CLDN18_total
400_786
0
0
4990
1.0
659.15
757.5


CLDN18_1_fake
0_400
0
3
3
0.16
0.38
0.0


CLDN18_1_fake
400_786
0
0
3946
1.0
521.23
551.5


CLDN18_2_fake
0_361
0
0
0
0.0
0.0
0.0


CLDN18_2_fake
361_747
0
0
3943
1.0
520.83
551.0


HPRT1
0_400
107
215
1259
1.0
167.77
175.0


HPRT1
400_793
0
0
974
0.98
126.40
101.0



Hint: This is just an example to illustrate the design of the table. Results may differ.
Columns in output file
While not using interval mode

name name of the input sequence
pos position of interest relative to input sequence
junc reads overlapping the position of interest
span read pairs spanning the position of interest
anch maximal number of bases next to position of interest that are overlaped by a single read
a reads mapping to sequence left of the position of interest
b reads mapping to sequence right of the position of interest

While using interval mode

name name of the input sequence
interval interval of interest relative to input sequence
overlap_interval_end_reads reads overlapping the end of the interval by at least BP_DISTANCE bases
span_interval_end_pairs read pairs spanning the end of the interval
within_interval reads mapping fully onto the interval
coverage_perc percentual coverage of the interval by aligned reads
coverage_mean average coverage per base for the interval (fold coverage)
coverage_median median coverage per base for the interval

Things to consider
When choosing the aligner you have to take into account that there are several differences among them:

STAR:

end-to-end alignment with no soft-clipping
very slow for small reference sequences
several available alignment parameters to optimize your results, which can be used while starting the pipeline


bowtie2:

end-to-end alignment might lead to insertions where the context sequence starts/ends
faster than STAR for short reference sequences (index creation parameters are calculated automatically)