I'm trying to estimate counts of numbers of substitutions along branches in a rooted triplet. I've used goldman_q_dna_triple() to get a q matrix, and I'd like to convert this into count data.
I thought I might be able to use Rates.toProbs then Probs.toCounts, but goldman_q_dna_triple() returns an ndarray, rather than a Rates instance, so I can't use Rates.toProbs.
Is there a simple solution to this? Can I convert the ndarray to a Rates instance easily? Or can I use another method to get counts from the output of goldman_q_dna_triple() (given that I know my branchlengths from elsewhere)?
Thanks for any help,
Probably asked too soon here. Thanks to some help from another poster, and reading the source code, I think I have this figured out. In case others are wondering the same, here's the solution I landed on: using Rates.toProbs() and Counts.fromTriple().
I have three species, two ingroup and an outgroup. I wanted a Q matrix, P matrix, and counts of changes for all three. I used this code:
spp1_seq = alignment.getGappedSeq(spp1)
spp2_seq = alignment.getGappedSeq(spp2)
outg_seq = alignment.getGappedSeq(outg)
spp1_seqm = ModelDnaSequence(spp1_seq)
spp2_seqm = ModelDnaSequence(spp2_seq)
outg_seqm = ModelDnaSequence(outg_seq)
# we get 3 things: the Q matrix, the count data, and the P matrix (probabilities) the ._data mean we just get an nd array
spp1_q = goldman_q_dna_triple(spp1_seq, spp2_seq, outg_seq)
spp1_counts = Counts.fromTriple(spp1_seqm, spp2_seqm, outg_seqm, DnaPairs)._data
spp1_p = Rates.toProbs(Rates(spp1_q, DnaPairs))._data
spp2_q = goldman_q_dna_triple(spp2_seq, spp1_seq, outg_seq)
spp2_counts = Counts.fromTriple(spp2_seqm, spp1_seqm, outg_seqm, DnaPairs)._data
spp2_p = Rates.toProbs(Rates(spp2_q, DnaPairs))._data
Log in to post a comment.
Sign up for the SourceForge newsletter:
You seem to have CSS turned off.
Please don't fill out this field.