<p>For future googlers, Ding et al., 2017 published a <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5476949/" rel="nofollow noreferrer">review</a>: "Comparison of Alternative Splicing Junction Detection Tools Using RNA-Seq Data". </p>

<p><a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5476949/table/T1/" rel="nofollow noreferrer">Table 1</a> has the list of software in this review. </p>

CC BY-SA 3.0
2018-02-12T13:39:25.533Z
2018-02-12T14:28:27.553Z
2018-02-12T14:28:27.553Z
aechchiki
<p>You can use <a href="https://github.com/lh3/minimap2/blob/master/misc/paftools.js" rel="nofollow noreferrer">the script</a> that comes with minimap2:</p>

<pre class="lang-sh prettyprint-override"><code># install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam &gt; output.txt
</code></pre>

<p>Some example output:</p>

<pre><code>P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672
</code></pre>

<p>Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.</p>

CC BY-SA 3.0
2018-02-27T10:54:49.447Z
Hi @user172818, I have two questions on this: (1) could you please explain what the value in `$3` represent exactly - and (2) I had cases where the value in `$7` is a list of more than one tuple: do you have an explanation ?
aechchiki
CC BY-SA 3.0
2018-02-27T12:49:50.627Z
@aechchiki 1) the exon number of the read. 2) an intron on a read may overlap multiple distinct introns in annotations due to alt splicing. The list gives all.
user172818
CC BY-SA 3.0
2018-02-10T00:10:54.673Z
2018-02-12T14:52:21.027Z
2018-02-12T14:52:21.027Z
user172818
<p>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.</p>

<p>I tried <code>Regtools</code>:</p>

<ol>
<li>feed the experimental alignment file (BAM: aligned transcripts to reference genome) to <a href="https://regtools.readthedocs.io/en/latest/commands/junctions-extract/" rel="nofollow noreferrer">Regtools junctions extract</a></li>
<li>use <a href="https://regtools.readthedocs.io/en/latest/commands/junctions-annotate/" rel="nofollow noreferrer">Regtools junctions annotate</a> to compare with the reference annotation (GTF)</li>
</ol>

<p>Exact commands I used for this are: </p>

<ol>
<li><code>[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed</code> - where the sorted BAM file has index in the same folder</li>
<li><code>[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out</code></li>
</ol>

<p>The alignments have been generated with <a href="https://github.com/lh3/minimap2" rel="nofollow noreferrer">Minimap2</a>.</p>

<p>It did not behave as expected: this gives me as output a file with only junctions annotated as 'new' (<code>$11==N</code>), while other annotations should be possible according to the <a href="https://regtools.readthedocs.io/en/latest/commands/junctions-annotate/" rel="nofollow noreferrer">manual</a>: </p>

<blockquote>
  <p>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.</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>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.</p>
</blockquote>

<p>My questions are: </p>

<ul>
<li>am I doing something wrong?</li>
<li>do you have any experience with <code>Regtools</code> and can help debugging this behavior?</li>
<li>do you have other suggestion to achieve the output I am seeking?</li>
</ul>

<rna-seq><software-recommendation><rna-splicing>
Exon-exon junctions: compare experimental transcripts to reference annotation
# Exon-exon junctions: compare experimental transcripts to reference annotation

(Asked by: username_0 on 2018-02-09)

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`:

1. feed the experimental alignment file (BAM: aligned transcripts to reference genome) to [Regtools junctions extract](https://regtools.readthedocs.io/en/latest/commands/junctions-extract/)
2. use [Regtools junctions annotate](https://regtools.readthedocs.io/en/latest/commands/junctions-annotate/) to compare with the reference annotation (GTF)

Exact commands I used for this are:

1. `[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed` - where the sorted BAM file has index in the same folder
2. `[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out`

The alignments have been generated with [Minimap2](https://github.com/lh3/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](https://regtools.readthedocs.io/en/latest/commands/junctions-annotate/):

> 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?

**Tags:** rna-seq, software-recommendation, rna-splicing

---

## 2 answers

### Answer by username_1 on 2018-02-10. Score: 4 | Accepted answer

You can use the script that comes with minimap2:
# install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam > output.txt

Some example output:
P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672

Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.

> **Comments:**
> **username_0**: Hi @user172818, I have two questions on this: (1) could you please explain what the value in `$3` represent exactly - and (2) I had cases where the value in `$7` is a list of more than one tuple: do you have an explanation ?  
> **username_1**: @aechchiki 1) the exon number of the read. 2) an intron on a read may overlap multiple distinct introns in annotations due to alt splicing. The list gives all.  

---

### Answer by username_0 on 2018-02-12. Score: 1

For future googlers, Ding et al., 2017 published a [review](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5476949/): "Comparison of Alternative Splicing Junction Detection Tools Using RNA-Seq Data".

[Table 1](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5476949/table/T1/) has the list of software in this review.


---

