Running Isoquant in transcripts mode from reads and from bams

[apollo994]

Hi,

Thanks for the active development of Isoquant!
This is more of a request for documentation than a bug report.

I'm testing the isoquant (version 3.7.1) feature to build models from trusted transcripts.
My goal is to understand:

• Weather is better to provide fasta of the transcripts or pre-map them
• Whether I should add a synthetic polyA to my transcripts

My workflow is the following.
This is not a real use case but more of a control.

1. I build models from ont reads using isoquant in ont mode
2. I extract the fasta sequence of the transcripts from the gtf
3. I add an artificial polyA to each transcript sequence (I also test with non-polyadenilated, named simple below)
   At this point, I have 11332 fasta sequences that correspond to transcripts.

At this point, my workflow branches into mapping on my own (a) and letting IsoQuant do the job (b).
1a. I map the sequencies (simple and polyA) with minimap2 -ax splice -t 16 "$REF" "$READS" > "$SAM"
2a. I run Isoquant from these mappings isoquant.py --reference "$REF" --bam "$BAM[simple,polyA]" --data_type transcripts
or
1b. I run Isoquant from reads isoquant.py --reference "$REF" --fastq "$FASTA[simple,polyA]" --data_type transcripts

In the slide below, you can see the results (where I'm also testing ont and pacbio presets).

In brief, if I provide BAM there is no difference between polyA or simple.
If I provide sequences of transcripts I do see a difference in the models
In all cases the transcripts that I get are less than the input ones.

My questions are:

• How and at what stage isoquant handle polyA?
• Why don't I get the same number of transcripts as the input? Especially in the case of bam as input?

Thanks again for opening and maintaining Isoquant.

Fabio

[andrewprzh]

Dear @apollo994

This is an interesting experiment and I'll try to give some insignts. However, there are many steps and it's hard to say what happens at which step. In general I would not expect IsoQuant to simply convert BAM file to GTF.

1. PolyA tails are used by the transcript model construction algorithm to estimate transcript end site coordinates. However, if the percentage of polyA is low, IsoQuant doesn't use them. Hence, both experiments (simple and polyA) yield a similar number of transcripts.

2. Although in this mode IsoQuant applies the lowest possible thresholds, there is no guarantee that every "read" will be converted to GTF. For example, IsoQuant ignores reads mapped with low quality.
   3Transcript discovery algorithm has a lot steps and filtering strategies, thus some transcripts may be considered unreliable, e.g. based on canonical splice sites. By default, IsoQuant will not output transcripts for which strand cannot be detected using splice sites or polyA tail.

3. BAM files vs FASTA. IsoQuant uses slightly different minimap options, so experiments carried out with FASTA files are also affected by the alignment step.

Best
Andrey