Reference loci missing from extended annotation

[apollo994]

Hi,

I'm using IsoQuant to quantify the amount of novel isoforms detected by two RNA seq experiment starting from a reference annotation.

This is my experimental set up:

• Reference genome of a butterfly
• Reference annotation automatically generated with Helixer (contains only one isoform per gene)
• RNA seq with ONT of tissue 1
• RNA seq with ONT of tissue 2

Here are the question I wanna answer:

• how many novel genes
• how many novel isoforms
• overlap between novel isoforms in the two tissue.

IsoQuant was my top choice as it allows to run with/without reference annotation, it provides the SQANTI output and it was relatively easy to set up.

Here is an example of my IsoQuant command for one of the tissues:

software/IsoQuant-3.3.1/isoquant.py \
	--reference reference.fa \
    --genedb helixer_ann.gff3 \
    --fastq RNAseq1.fastq \
    --data_type nanopore \
    -o out_folder/  \
    --sqanti_output --threads 16

This resulted in two *.extended_annotation.gtf that I checked with gff compare. Here is the result:

#= Summary for dataset: helixer_ann.gff3 
#     Query mRNAs :   16102 in   16102 loci  (14527 multi-exon transcripts)
#            (0 multi-transcript loci, ~1.0 transcripts per locus)

#= Summary for dataset: RNAseq1.extended_annotation.gtf 
#     Query mRNAs :   20115 in   14893 loci  (18740 multi-exon transcripts)
#            (3155 multi-transcript loci, ~1.4 transcripts per locus)

#= Summary for dataset: RNSseq2.extended_annotation.gtf 
#     Query mRNAs :   18330 in   13136 loci  (17166 multi-exon transcripts)
#            (3068 multi-transcript loci, ~1.4 transcripts per locus)

How is possible that the total number of loci decreases?
In principle the extended annotation should contains all reference + novel.

I tought that one reason could be due to close genes that get "fused" in one locus by adding the RNA seq knowledge.
This is not the case, as by inspecting the gtf on IGV I could find many cases where the reference locus was completely missing in the extended annotation, or the annotated isoform was shorter.

This is also confirmed by extracting genes in the same interval from the three gtf:

for a in anns:
    print (a)
    print(anns[a]['ilHelHell1.1_Chr7',5452963:5485951].gene_id.unique())

helix
['_ilHelHell1.1_Chr7_000057' '_ilHelHell1.1_Chr7_000058'
 '_ilHelHell1.1_Chr7_000059' '_ilHelHell1.1_Chr7_000060'
 '_ilHelHell1.1_Chr7_000061']
RNAseq1
['_ilHelHell1.1_Chr7_000059']
RNAseq2
['_ilHelHell1.1_Chr7_000058']

Do you have any explanation for this?
Please let me know if you need any further information to troubleshoot this.

In the meantime, many thanks for sharing your project!

F

[andrewprzh]

Dear @apollo994

This is known bug in generation of extended annotation, it will be fixed in the nearest release!

Best
Andrey

[andrewprzh]

Finally released new version 3.4, which fixes this issue.