Improving the class2tree resolution

[gedankenstuecke]

The function class2tree converts a list of taxonomic classifications into a tree. So far it only takes ranked levels into account. This leads to extremely multifurcations, as there are only so many ranked taxonomic levels, which limits the use of the resulting trees.

In comparison, this is how the NCBI Taxonomy tree for the same taxa looks like:

Based on an idea by @trvinh (also pinging as he's the one who had the issue in the first place) I started to explore the idea of using the total vectors of classifications for calculating the pairwise distances. The quick & dirty idea:

Calculating the distance between two taxa A and B as (A ∪ B) - (A ∩ B), reasoning that the closer the intersection of the two classifications approaches the union of it, the smaller the distance. An example implementation of that is here.

This yields a better resolved tree, alas the branchings are not fitting the taxonomy any longer. The reason for that in the example: Drosophila and Anopheles have a very detailed, shared sub-classification (starting at the phylum level), which joins them to the exclusion of basically all other taxa.

I guess my tl;dr is: The current class2tree-solution is too low resolution to resolve meaningful trees. My proposed alternative has it's own issues. What can we do to make it better? 😄

[sckott]

Thanks for this @gedankenstuecke Will think about this and try out your solution

pinging @Edild wrt class2tree

[gedankenstuecke]

Thanks, and I haven't heard from @trvinh, but he might have some new insights as well! 👍

[sckott]

any thoughts on this @Edild ?

[gedankenstuecke]

Let me push to @trvinh as well, I vaguely remember him telling me something about how he fixed this.

[sckott]

pinging @jarioksa in case he has any opinions

[trvinh]

I did try to create a matrix that contains equal number of taxonomy ranks for all taxa, including noranks. The missing rank IDs in the matrix were complemented by the previous available IDs. Then using that matrix to sort the taxa, I got what I expected to have.
However, I have to use a perl script to create that matrix. It would be great if it can be done directly in R :)

[gedankenstuecke]

Ah, let me see whether I got this: So if you have the taxonomy-IDs like this

Species A: 1,2,3,4,5,6,7,8,9
Species B: 1,2,3,10,11

with 4,5,7,8 being unranked you will fill them up to

Species A: 1,2,3,4,5,6,7,8,9
Species B: 1,2,3,3,3,10,10,10,11

?

[trvinh]

yes, almost correct :)
However, I also fix the relative positions of all possible main ranks from strain to superkingdom (regardless to the taxa of interest), and also the relative positions of each individual taxon rank IDs.

In your example, as 9,10,11 are main ranks, I already knew their relative positions: 9->10->11, therefore

Species A: 1,2,3, 4 , 5 , 6 , 7 , 8 , 9 ,(9),(9)
Species B: 1,2,3,(3),(3),(3),(3),(3),(3), 10, 11

[gedankenstuecke]

Ah, okay! Is your perl script somewhere on GitHub so that people can have a look and potentially adapt it to R?

[trvinh]

yup, here it is
sorry for the dirty code :P

[sckott]

@gedankenstuecke any thoughts on this?

[gedankenstuecke]

@sckott Due to thesis writing I didn't have any time to look into this in depth. I remember @trvinh and I discussed moving his solution from Perl to R, but I have no clue what the status is for that.

[sckott]

okay, thanks for update. will wait and see if something happens - will leave open for now

[gedankenstuecke]

Thanks! I will be back in the office on Monday and then @trvinh and I can discuss how to best convert his code. I think for the tool he's writing it would be necessary/ideal to get a native R-solution as well.

[sckott]

Great, 👌

[jarioksa]

I quite can't follow this. The current solution may not be satisfactory, but I think it is based on the data: if you have taxa that are separated at the same level, they will have "multifurcations". That is, if you have only information about families, all genera within the family are separated at the same level. You can produce fancier plots, but they are artifactual. If I get something reproducible, I may think differently. Do you have some information in addition to classification which should be taken into account?

[trvinh]

@jarioksa: I hope that I understood your point correctly. Yes, you are right, if we have only information about families, we will have multifurcation for taxa within a family.
With our current idea, if the information about the splits are missing, we still have multifurcations. We are not trying to remove multifurcations by creating artifactual trees (plots). Instead, we just want to solve them as most as we can by using more information (here, the IDs of "norank" levels between "main" levels strain, species, genus, ....).

[jarioksa]

So it's about not using existing information. Got to see how the information is extracted.

[trvinh]

So, now we can reconstruct the NCBI taxonomy tree successfully :)

Our solution is that, we first index ALL available levels (also including the unranked) from input data in ascending order. Then create a data frame with columns are all levels and rows are the taxa. After all, we arrange the IDs for each taxon to that data frame. For a certain level where we found a missing ID, we assign its value (ID) with the ID of previous level. By that all taxa will have the same number of level IDs (all rows have the same length). That data frame is used as input for taxa2dist() .

Since I am still beginner in R, my code looks very ugly :-P I'm really appreciate if you could improve it. Many thanks in advance!!

[sckott]

thanks @trvinh !

@jarioksa can you chime in with any thoughts you have on the above solution/code?

[jarioksa]

Got to have a look at the details of the code. However, the solution is the correct one.

[gedankenstuecke]

I suspect there might be some ways to make the code much more R-like. Given that @trvinh and I are rather new to R (coming from Python and Perl) we learned quite a bit about R data types in the process. 😂

[sckott]

@gedankenstuecke and @trvinh - are you meaning to modify code more? you can send a PR, or would you rather i just add it?

[trvinh]

@sckott : we will send a PR next week.

[sckott]

sounds good

[gedankenstuecke]

Closed through #634 😊