Counting and extracting invariant sites question. about raxml_ascbias HOT 5 CLOSED

RaxML considers those invariant sites because the N could be an A, T, G, or C. My script should remove those two sites, which is what it needs to do because RaxML will crash otherwise. It's not the normal definition of monomorphic sites but that's how RaxML handles them.

from raxml_ascbias.

I'm trying to think what else could be causing the inflated counts because their sum (sum of A/C/G/Ts) is basically the number of sites of the entire bacterial genome: roughly 2.2M. If all the sites are invariant that would imply the isolates are clonal which they are not.

from raxml_ascbias.

It could be gaps. I'm not sure on that but if you want I can take a look this weekend. I seem to remember that this script might not be perfect in all datasets, and i think maybe gaps throw it off. I'm normally very thorough with my coding and testing practices but this was one of my early scripts when I was learning python and I've noticed it sometimes has issues. I also wrote this as a quick and dirty script for a specific need that I had.

If it's the same dataset you sent me the other day I'll go ahead and look at it. I'm currently not super speedy about getting to this stuff lately because I have a newborn baby at home but I'll try to work on it this weekend.

from raxml_ascbias.

It could be gaps. I'm not sure on that but if you want I can take a look this weekend. I seem to remember that this script might not be perfect in all datasets, and i think maybe gaps throw it off. I'm normally very thorough with my coding and testing practices but this was one of my early scripts when I was learning python and I've noticed it sometimes has issues. I also wrote this as a quick and dirty script for a specific need that I had.

If it's the same dataset you sent me the other day I'll go ahead and look at it. I'm currently not super speedy about getting to this stuff lately because I have a newborn baby at home but I'll try to work on it this weekend.

Congrats on the baby! I totally understand. I can already tell you that the whole genome SNP alignment was generated in snippy and there are gaps and Ns. I'm awful with python so kudos to you. I might need to send you a dropbox link though to another file. I'll probably send that later tonight. Thanks!

from raxml_ascbias.

Okay I am closing this for real now. You have been a big help. After thinking about this more...I found the solution. I used snp-sites -cbp outputprefix input.phylip option to extract the monomorphic sites ACGT only and export in a phylip file. Then I used the phylip file as input into your script. The results are the exact same as the internal script. So everything is fine. Enjoy time with the new baby. You can ignore the .aln file I sent. Take care!

from raxml_ascbias.