Count data from goldman_q_dna_triple

Help
2012-03-27
2013-04-29
  • Rob Lanfear

    Rob Lanfear - 2012-03-27

    Hi,

    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,

    Rob

     
  • Rob Lanfear

    Rob Lanfear - 2012-03-29

    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

    Cheers

    Rob

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks