Repository navigation
Novel/Unannotated Isoforms #300
Description
Hi there,
So I'm trying to look at novel/unannotated isoforms detected and quantified by Bambu and I'd just like some clarification and/or advice. So basically I get an extended annotation GTF file, and upon looking at it, it contains all the Ensemble transcript IDs that were in my reference GTF, plus several hundred Ensemble transcript IDs that were not in my reference GTF. Where is Bambu getting these Ensemble transcript IDs from? I was hoping unannotated transcripts would be given some kind of separate ID instead of being associated with existing annotations.
Activity
andredsim commentedon Aug 16, 2022
Hi there,
Sorry for the late response. Bambu shouldn't be able to create new ensembl transcript IDs. All novel transcripts will be assigned a transcript id of tx.XXX where XXX is an arbitrary id. They may be assigned an ensembl gene ID if they overlap with an annotated gene, but these ID should already be present in the annotation GTF file.
If this doesn't answer your question, you could send me the files, or a screen shot of an example, so that I may look a bit deeper.
ajritter8 commentedon Aug 16, 2022
andredsim commentedon Aug 17, 2022
Hi Alex,
There are many things that could influence this. Which parameters are you using for your bambu run and is this a standard long-read run or a custom library preparation (i.e used primers to target specific genes)?
I would recommend just using the NDR (novel discovery rate) parameter and if you want a more sensitive output increasing the NDR. In a human sample we generally use between .1-.2, however if you are using a different species with less comprehensive annotations you might want to increase this value. The maximum this value can be set is 1, which will mean every novel transcript that has at least 2 reads will be annotated.
In the coming patch (see the BambuRevisionManuscript branch) we are also introducing a recommended NDR feature which will recommend an NDR value for your sample, in the event many annotations are missing.
ajritter8 commentedon Aug 18, 2022
andredsim commentedon Aug 19, 2022
Hi Alex,
This is very strange, if the NDR is 1, you should be getting loads of novel transcripts as that is maximum sensitivity.
Could you let me know how you changed minGeneFrac, and minReadCount that allowed you to get some unannotated transcripts? This might give me a clue as to what might be occuring.
Could you also please try a run where you replace makeTxDbFromGFF(gtf_file) with prepareAnnotations(gtf_file) as a sanity check for me?
In the likely event this doesn't make a difference, could you add rcOutDir = "/some/output/directory/here" to your arguments. This should output an .rds file which contains all the read classes for your sample (essentially all possible candidates before any filtering). If you could send that to me I could have a look what might be going on.
ajritter8 commentedon Aug 19, 2022
andredsim commentedon Aug 22, 2022
Hi Alex,
The default threshold for min.readCount is 2, meaning only transcripts with 1 full length read are filtered out via this threshold. Your threshold of 5 is therefore less sensitive than the default. (Note that it is >= not <=)
The default min.readFractionByGene is also 0.05.
Let me know once you can upload the rcOutDir, it gets produced before the quantification step, so it make be in the directory even before it finishes running.
I also wanted to ask about what happened with your custom annotation file when using prepareAnnotations, as it is not expected behavior that it should not import anything other than the .gtf that is provided to it. How did you know it loaded in a standard annotation and not your custom annotation?
Could you please check the annotations object you have created with for me with TxDbFromGFF
mcols(annotations)and just double check that the TXNAME and GENEID are not the same, and that there are multiple rows that have different TXNAMES but share a GENEID. We have had reports that makeTxDbFromGFF sometimes causes all transcript and gene names to be the same resulting in no multi-isoform genes. This is a big problem for bambu and might be the crux of the issue.Thanks, I hope we will be able to get this working for you soon!
Andre
ajritter8 commentedon Aug 22, 2022
andredsim commentedon Aug 23, 2022
Hi Alex,
Thanks for sending the rcFile. I havn't finished investigating this, but I will update you with what I have seen so far. I had a look at the rcFile and Bambu reports 523982 read classes (with 2 or more reads), with 14283 being already annotated and the remaining 509699 being novel candidates.
When I look at the transcript probability score of these read classes, nearly all the novel candidates are being scored very lowly, with only 2725 read classes being > 0.5.
After doing extend annotations (where additional filtering is performed including the NDR and min.readFractionByGene) I am left with 874 novel genes with default threshold and NDR = 1 in the final output.
This was with my own human genome and annotation file (which might explain different results), however while it sounds higher than what you are reporting it is still lower than I would expect so I am going to investigate further and get back to you on this.
Kind Regards,
Andre Sim
ajritter8 commentedon Aug 23, 2022
andredsim commentedon Sep 5, 2022
Hi Alex,
Sorry again for the belated response. I was a bit busy the last week and didn't get a chance to investigate this further until today.
I found that the majority of candidate novel transcripts in this sample appear to be subsets of known transcripts. That is transcripts which share a continuous run of exons with a known transcript and could also be produced by degradation of the transcripts. By default bambu does not output subset transcripts due to the high chance they may be false positive. However, if I set
opt.discovery = list(remove.subsetTx = FALSE)and output them at an NDR of 1 I get 14606 novel transcripts.As a side experiment I also set
opt.discovery = list(min.readFractionByGene = 0). This resulted in a boost of ~5000 novel transcripts at an NDR of 1. These would be transcripts that have at least 2 reads, but stem from a gene with very high expression that dwarves the expression of these novel isoform candidates. You can combine these two options like soopt.discovery = list(remove.subsetTx = FALSE, min.readFractionByGene = 0)which at an NDR of 1 results in 335863 novel transcripts accounting for nearly every possible candidate.Therefore if you are looking to boost the number of novel transcripts you could consider allowing the subset transcripts to be output as well. (I would recommend an NDR of 0.1-0.2 and not 1).
To answer your question. It is technically possible to filter the rcFile yourself but it is not very user friendly and the quantification statistics are not yet included at that step.
You can view the rcFile data with rowData(rcFile).
I believe it is unlikely that writeBambuOutput() is the cause. To check this however you can view the annotations before they are output by looking at the rowRanges of the se output from bambu. rowRanges(se). This should have the same number of annotations in it as the output GTF file.
Let me know if this doesn't solve your issue and I can reopen this.
Kind Regards,
Andre Sim
ajritter8 commentedon Oct 11, 2022