Implementation of additional processing steps in postoga following current TOGA dir structure

[alejandrogzi]

Hi @MichaelHiller and @kirilenkobm,

I am working in the next steps of postoga. Some of the previous reported problems ocurred because I was taking the old output dir structure. Would it be possible for you to share a random output dir (human-to-some species; preferred mammal)?. I do not have access to an HPC anymore, so I can't run TOGA at the time.

Hope everything is ok,
Alejandro

[MichaelHiller]

no problem. I have tar.gz a full TOGA run with the temp dir and everything else, and this will be on https://genome.senckenberg.de/download/forAlejandroTOGAdir/ in a minute.

Thx

[kirilenkobm]

I believe worth documenting too https://github.com/hillerlab/TOGA/wiki/Output-directory-structure

[alejandrogzi]

@MichaelHiller, @kirilenkobm,

Just finished all the new changes to postoga.

I wanted to quickly note that the bug in #136 about bed2gtf:

we have tested bed2gtf on a TOGA run and get warnings/error messages like
Isoform ENST00000347757.SSX5.4667 not found. Check your isoforms file

Is now corrected. It did not have anything to do with bed2gtf directly, but with how postoga merged some dataframes in order to include all gene -> projection pairs in the isoforms file. Now if you run it through postoga, all the lost projections will be included in the resultant .gtf file and the program will run ok. Sorry for any inconvenience there.

postoga now follows the most recent dir structure pointed out by Bogdan. Also, new features have been added, including: 1) ortholog length distributions, 2) orthology score distributions, and 3) pseudo-BUSCO completeness analysis.

1. Ortholog length distributions: only adds up ~10s at runtime and works through NOEL, a program I wrote a few months back to calculate this type of stuff. All the configurations are also automated.

2. Orthology score distributions: Instead of using the deprecated transcript_quality.tsv, postoga now uses the prediction scores directly. Plotting this information could be useful to quickly see how confident is TOGA about your newly predicted transcripts.

3. pseudo-BUSCO completeness: Instead of running BUSCO from scratch, postoga uses a set of curated BUSCO datasets (under eukaryota_odb10, for example), to compare your assembly and extract its completeness (in the case you have declared some filters, the filtered bed will be used to calculate the completeness). At this moment there are only 4 DBs curated (eukaryota, vertebrata, mammalia and eutheria) and 3 about to be curated (aves, carnivora and passeriformes). The default DBs are eukaryota and vertebrata. In the case a user wants to evaluate a bird, eukaryota, vertebrata, aves and passeriformes are used.

I also have tried to fix the nomenclature problem (some users that prefer entrez IDs or gene names instead of ensembl IDs). A --source flag has now been introduced to ask the users if they want to use an ensembl, entrez or gene_names background. This option maps to your Ancestral placental DB (now with gene_names and entrez IDs) and also with the BUSCO DBs (also with gene_names and entrez IDs).

I think this wraps up most of the changes. However, is probable that some error/bugs might have slipped through my tests. All the possible feedback will be very appreciated!

Best,
Alejandro

[MichaelHiller]

Fantastic. We will give it a try.
Thx a lot
Michael

[MichaelHiller]

One question. How do you determine which BUSCO ortholog corresponds to a TOGA annotated gene?

[alejandrogzi]

@MichaelHiller,

One question. How do you determine which BUSCO ortholog corresponds to a TOGA annotated gene?

Thanks for this question. There are some caveats here, that is why I called "pseudo-BUSCO":

1. Each curated DB has 3 sources: ensembl IDS (based on GRCh38), gene names, and entrez IDs. This indirectly assumes that the human should be the reference species and the isoforms provided are ensembl IDs (if using -src ensembl; which is default). With the data TOGA provides (e.g. ENST00000420798.TRO) or normal IDs (e.g. ENST00000420798), postoga builds this for your annotation (if you are using gene_names or entrez IDs, this df would look very similar):

t_gene	helper	transcripts	orthology_class	class	pred	q_gene
ENSG00000239697	ENST00000293825.TNFSF12	ENST00000293825.TNFSF12.22	one2one	I	0.996369	reg_10433
ENSG00000185298	ENST00000329214.CCDC137	ENST00000329214.CCDC137.22	one2one	I	0.996369	reg_10200
ENSG00000111911	ENST00000229633.HINT3	ENST00000229633.HINT3.3	one2one	I	0.996369	reg_20375
ENSG00000119707	ENST00000261973.RBM25	ENST00000261973.RBM25.13	one2many	I	0.996369	reg_17756

2. At the pseudo-BUSCO step, postoga just overlaps all the genes in your annotation ["t_gene" column] (if you specified some filters, postoga will use those filtered genes) with each one of the BUSCO DBs (if mammal -> eukaryota, vertebrata, mammalia, eutheria; if birds -> eukaryota, vertebrata, aves, passeriformes).

3. The BUSCO DBs where made by first running BUSCO over the human assembly (GRCh38) (excluding aves and passeriformes, which were performed with the chicken assembly), converting protein IDs to their respective gene IDs (ensembl, gene name, entrez) and filling out the missing ones through OrthoDB API. Is important to note that this approach yielded a reasonable % of BUSCO genes annotated.; however, there are some differences between nomenclatures (ensembl vs gene names vs entrez). For some DBs, the complete set of mapped ensembl IDs is indeed present but the amount of gene names is reduced by 10-15. For other DBs, the opposite case happen: 98% of the whole dataset is correctly annotated with gene_names but a lower % is annotated with ensembl IDs. The difference is not that high, though (~2-3%).

4. From this explanation it is now understandable why the plot does not communicate % missing, % fragmented, % duplicated and % complete (single-copy).

5. This will become obsolete if a user runs TOGA with a species different than human as a reference and ensembl IDs as their isoforms source. Conversely, if a user runs TOGA with a species different than human BUT uses gene_names or entrez IDs as their isoforms source, this would probably work. This is one the bigger problems of this approach but I think it could be improved (maybe with some ideas that I am not considering).

Hope this clarify your question. Any feedback will be amazing to hear!

Alejandro

[kirilenkobm]

It seems more appropriate to move this issue to the discussions section and give it a more suitable name. :)

[MichaelHiller]

Got it. Alejandro intersects BUSCO gene IDs (Ensembl) with the Ensembl IDs that TOGA annotated. Makes sense. Thx for explaining it.

I agree, this should be documented and quickly explained, likely on the wiki.
We may recommend that this should not be used as publication-ready stats (actually our ancestral gene classification, Fig 6 in the Science pub, is much more informative), but it provides a quick overview.

[alejandrogzi]

Sweet! With this update postoga outputs both: ancestral placental (bartplot with % of ancestral classes and scatterplot with inactivated vs missing ancestral sequences, just like in your paper) and BUSCO, a quick overview of both gives you a solid initial interpretation of TOGA results.

I can make myself some time to write that and specify all features, recommendations and stuff. Please let me know if you agree with that (this would also be in the README?). Also, I wanted to note that the last postoga release had some problems with dependencies (at least when I tried it in some VMs here). Solved that and build a conda env .yml to make it work easier.

@MichaelHiller I send you an email yesterday I think, please let me know if you read it.