Using the permulations method on gene trees with paralogs?

[00-kelvin]

Hi there!

I'd like to make use of your permulations method to generate null distributions for an analysis I'm working on about gene loss and duplication in spiders. Because one of my aims is to investigate duplication events, I am looking at many genes in which there is more than one copy for some species. This means the gene trees often have multiple tips per species. And, even without considering paralogs, I suspect that the individual gene trees do not all necessarily match the species tree in terms of topology. (These are the Resolved Gene Trees generated by Orthofinder.)

I understand from this discussion that the readTrees function wants all genes to have the same topologies. So one solution I imagine would be to use a tool like OrthoSnap to extract single-copy trees from each of my orthogroups, then use your estimatePhangornTreeAll wrapper function to generate a readTree-compatible species tree from all the individual pruned gene trees.

However, since for my purposes all I really need is the permulated phenotype vectors to use as the null for a separate method, I'm wondering if there is a way for me to convert an existing species tree with branch lengths in Newick file format to a format that could be provided directly to the simBinPheno and simBinPhenoVec functions, rather than having the readTrees function generate the master tree for me? This would save me the effort of pruning all my gene trees to be single copy and using them to make a master tree, which after all I am not sure would produce as reliable of a species tree/branch lengths as the one I already have (also from OrthoFinder).

Alternatively, if you have suggestions of any other strategy by which to generate the permulated trees when using genes which have multiple copies in some species, please let me know!

Thank you,
Calvin

[ekopania]

Hi Calvin,

Yes, you can run simBinPheno and simBinPhenoVec with only a species tree, it just has to be formatted like a master tree in a trees object:

speciesTree<-read.tree("path_to_tree_in_Newick_format") testTrees<-list() testTrees$masterTree<-speciesTree simBinPheno(testTrees, root, phenvec) simBinPhenoVec(testTrees, root, phenvec)

It should be fine that there are no RER trees in your trees object, as long as it is formatted like a list and contains your species tree labeled as masterTree, since the simulation functions do not use any other information from the trees object.

-Emily

[00-kelvin]

Hi Emily,

Thank you so much, this is really helpful!

I was wondering if you had a sense for how long it should take to generate a single permulated tree, given a species tree of ~100 species? I have just tried running the program in the way you suggest, and it has been stalled (with no error message or output at all) on the simBinPheno function for at least 10 minutes. With the goal being to generate ~10,000 of these permulated trees, I wondered if this was due to some error rather than just the normal processing time for the function. (my script below for reference)

Thanks again!
Calvin

library(RERconverge)
library(phangorn)

rerpath <- find.package('RERconverge')

treefile <- "~/orb-selection/assets/SpeciesTree_N5_outgroupDM_brlen.nwk"

#define the root species
root <- "Drosophila_melanogaster"

#read in the species tree
speciesTree <- read.tree(treefile)

#make empty list - don't need gene trees for permulation
testTrees <- list()

#add species tree to "list" of trees
testTrees$masterTree<-speciesTree 

#vector of foreground tips (orbweavers)
orbweavers <- readLines("~/orb-selection/assets/orbweavers-list.txt")
allspecs <- testTrees$masterTree$tip.label
phenvec=ifelse(allspecs %in% orbweavers, 1, 0)
names(phenvec)=allspecs

permbintree <- simBinPheno(testTrees, root, phenvec)
permbinpheno <- simBinPhenoVec(testTrees, root, phenvec)

[00-kelvin]

Update: still hung up after several hours. Attaching my tree, foreground list and updated script in case it's of any use. Would be very grateful for any advice!

permulations.zip

[ekopania]

Hi Calvin,

The simBinPheno and simBinPhenoVec functions run older versions of permulations that try to match the exact tree structure of the foreground species, which can hang up if there are a lot of foreground taxa or if there are many sister species among the foreground taxa.

For your data, I would recommend using the new categorical permulations, which can be run as follows:
testCatPerm10<-categoricalPermulations(testTrees, phenvec, rm="ER", rp="auto", ntrees=10)

rm="ER" uses a rate matrix that assumes equal transition rates between states. 'rp="auto"' is to calculate the prior probability distribution of states at the root. Both of these are for castor's fit_mk function, which is called by categoricalPermulations. There is more information here if you would like to change the default parameters: https://www.rdocumentation.org/packages/castor/versions/1.6.1/topics/fit_mk

ntrees is the number of permulations.

The output will be a list. To get your permulated trees, use testCatPerm10$trees.

Hope this helps!

-Emily

[00-kelvin]

Thank you so much, this worked!!

One small note: the output tip values for each tree seem to default to "1" and "2" even if I provide values "0" and "1" in the phenotype vector. I'm guessing this is because it's designed to take as input any number of categories (presumably numbered 1, 2, 3...), not just 1s and 0s? Just thought I'd mention in case this wasn't the intended behavior.

Otherwise the tool worked perfectly, thank you again!

-Calvin

[ekopania]

Great, I'm glad it worked!

I believe outputting "1" and "2" is the expected behavior for categoricalPermulations, and you're correct that this is a quirk of the function being written for any categorical traits, not just binary traits. Thanks for letting us know!