My aim is to parse an experimental transcript set (obtained by RNAseq) to check which splice junctions are already reported in a reference annotation and which ones are new.
I tried Regtools
:
- feed the experimental alignment file (BAM: aligned transcripts to reference genome) to Regtools junctions extract
- use Regtools junctions annotate to compare with the reference annotation (GTF)
Exact commands I used for this are:
[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed
- where the sorted BAM file has index in the same folder[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out
The alignments have been generated with Minimap2.
It did not behave as expected: this gives me as output a file with only junctions annotated as 'new' ($11==N
), while other annotations should be possible according to the manual:
DA - This exon-exon junction is present in the transcriptome provided by the user(GTF) The ends of this junction are hence known donor and known acceptor sites according to the GTF file.
NDA - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are known donor and known acceptor sites according to the GTF file.
D - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The donor of this junction is a known donor but the acceptor is novel.
A - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The acceptor of this junction is a known acceptor but the donor is novel.
N - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are hence not known donor/acceptor sites according to the GTF file.
My questions are:
- am I doing something wrong?
- do you have any experience with
Regtools
and can help debugging this behavior? - do you have other suggestion to achieve the output I am seeking?