Multiple Trees in File

[mmckibben]

Good afternoon! Sorry to leave another response, I have been trying for several days to trouble shoot this on my own with no luck.

When running the sample data, the output from the align subcommand gives me CBF3_MSA.fasta and CBF3.tree instead of CBF3_tree.tree. I checked the .log files and I don't see an error happening. I then tried executing the predict subcommand using this CBF3.tree and it gives me the error that there are multiple trees in this file. I have tried looking at the Hyphy/Biopython side and the PASTA side of the problem, and noticed there was an extra semicolon in the tree. I deleted that manually and tried running again with the same error occurring. At this point I assume some step in the PASTA align/tree component is being skipped and just not giving me an error. I have tried reinstalling things and messing with the versions for each dependency without any luck. Bellow are the end portions of the logs:

Align:

PASTA INFO: TreeShrink option has been turned off!
PASTA INFO: Step 4. Realigning with decomposition strategy set to mincluster
PASTA INFO: Step 4. Alignment obtained. Tree inference beginning...
PASTA INFO: realignment accepted despite the score not improving.
PASTA INFO: current score: -17739.33, best score: -17698.421
PASTA INFO: TreeShrink option has been turned off!
PASTA INFO: Writing resulting alignment to /tmp/pastajob_1572975010.129370.marker001.BAD_Mutations_PastaInput_KxK2_1.aln
PASTA INFO: Writing resulting tree to /tmp/pastajob_1572975010.129370.tre
PASTA INFO: Writing resulting likelihood score to /tmp/pastajob_1572975010.129370.score.txt
PASTA INFO: The resulting alignment (with the names in a "safe" form) was first written as the file "/tmp/pastajob_1572975010.129370_temp_iteration_1_seq_alignment.txt"
PASTA INFO: The resulting tree (with the names in a "safe" form) was first written as the file "/tmp/pastajob_1572975010.129370_temp_iteration_1_tree.tre"
PASTA INFO: Total time spent: 70.6753809452s
===2019-11-05 10:31:22,141 - LRT_Predict===
DEBUG stderr:
===2019-11-05 10:31:22,169 - Pasta_Align===
DEBUG Sanitizing alignment in /tmp/BAD_Mutations_BackTranslated_Aw55AW.fasta
===2019-11-05 10:31:22,169 - Pasta_Align===
DEBUG Sanitizing tree in /tmp/pastajob_1572975010.129370.tre
===2019-11-05 10:31:22,243 - LRT_Predict===
INFO Nucleotide alignment in /tmp/BAD_Mutations_BackTranslated_Aw55AW.fasta
===2019-11-05 10:31:22,243 - LRT_Predict===
INFO Tree in /tmp/pastajob_1572975010.129370.tre
===2019-11-05 10:31:22,261 - LRT_Predict===
INFO MSA copied to Output_Dir/CBF3_MSA.fasta
===2019-11-05 10:31:22,261 - LRT_Predict===
INFO Tree copied to Output_Dir/CBF3.tree<

Predict:

===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable BASE to /scratch/BAD_Mutations_Data
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable TARGET_SPECIES to hordeum_vulgare
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable EVAL_THRESHOLD to 0.05
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable BASH to /bin/bash
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable GZIP to /bin/gzip
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable SUM to /usr/bin/sum
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable TBLASTX to /usr/bin/tblastx
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable PASTA to /usr/local/bin/run_pasta.py
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Setting variable HYPHY to /usr/bin/hyphymp
===2019-11-07 08:44:27,234 - Configuration_Handler===
DEBUG Command line and config options merged. Values: {'sum_path': '/usr/bin/sum', 'target': 'hordeum_vulgare', 'gzip_path': '/bin/gzip', 'loglevel': 'DEBUG', 'config': 'BAD_Mutations_Config.txt', 'evalue': '0.05', 'tree': 'Output_Dir/CBF3.tree', 'substitutions': 'Test_Data/CBF3.subs', 'base': '/scratch/BAD_Mutations_Data', 'hyphy_path': '/usr/bin/hyphymp', 'action': 'predict', 'output': 'Predictions_Dir', 'fasta': 'Test_Data/CBF3.fasta', 'tblastx_path': '/usr/bin/tblastx', 'bash_path': '/bin/bash', 'alignment': 'Output_Dir/CBF3_MSA.fasta', 'pasta_path': '/usr/local/bin/run_pasta.py'}
===2019-11-07 08:44:27,234 - LRT_Predict===
DEBUG Checking if BAD_Mutations_Config.txt exists.
===2019-11-07 08:44:27,235 - LRT_Predict===
DEBUG Checking if Output_Dir/CBF3.tree exists.
Traceback (most recent call last):
File "./BAD_Mutations.py", line 372, in
main()
File "./BAD_Mutations.py", line 330, in main
arguments_valid, msg = parse_args.validate_args(config_opts, loglevel)
File "/home/arizonica/BAD_Mutations/lrt_predict/General/parse_args.py", line 347, in validate_args
if not parse_input.valid_tree(args['tree'], log):
File "/home/arizonica/BAD_Mutations/lrt_predict/General/parse_input.py", line 29, in valid_tree
p = Phylo.read(f, 'newick')
File "/usr/local/lib/python2.7/dist-packages/Bio/Phylo/_io.py", line 73, in read
"There are multiple trees in this file; use parse() instead.")
ValueError: There are multiple trees in this file; use parse() instead.<

Thank you again for your assistance, I very much appreciate your time!

[mmckibben]

update: I changed the biopython portion from read() to parse() and was able to run the tree, however I am getting p-values of 0 for every SNP in the test files.

[TomJKono]

Hello,

The first issue (multiple trees in file) looks like a real bug - I will start working on the tree sanitization script to address that. The second issue regarding P-values of 0 seems like more of a complicated case. I will look into that once the tree sanitizing is in place and I can reliably reproduce the problematic files. Thank you, and apologies for the inconvenience!

[mmckibben]

Thank you for your help, please let me know if you would like more information or need any help! As an update I tried running the tree again with 42 species instead of all the 90 (what currently downloads) and now get proper p-values, however all of them are significant, even for those in the test data which should not be.

[TomJKono]

Hello,

Thank you for your patience! I have pushed some edits to the dev branch that should address the tree parsing issue. With regard to the P-value issue, can you share the output of the predict and compile commands? I can make sure that the columns of the HYPHY output are being parsed and handled properly. Thank you!

[mmckibben]

Thanks again for the help! I have attached all the files and logs associated. Please note if I run the predict command on the provided tree and MSA fasta file I get similar p-values to those in the sample output files. I also tried running both commands on some of my own data and it gives highly significant p-values (e^-10) for all possible snps, which seems somewhat suspicious to me.

CBF3.tree.txt
CBF3_MSA.fasta.txt
CBF3_Alignment.log.txt
CBF3_Predictions.log.txt
messages.log.txt

[TomJKono]

Hello,

Thank you for sharing the files! It looks like the alignments and tree files are generating properly and the HYPHY log does not show any problematic output. The issue might be that the logistic regression model for generating P-values does not really translate to analyses that use a different set of query genomes than the ones for which it was developed. My hunch with that is that the total substitution rate observed depends on the number of sequences in the alignment.

I think the easiest way to get a usable prediction from the output of BAD_Mutations when using a different set of sequences, you may be able to use the heuristics that were originally used in the Chun and Fay 2009 (https://genome.cshlp.org/content/19/9/1553.abstract) paper:

Deleterious mutations were predicted by nonsynonymous SNPs that disrupt significantly constrained codons defined by the LRT (P < 0.001) and a number of subsequent filters (Supplemental Table S1). First, positions with low power, <10 eutherian mammals, were eliminated. Second, a small number of sites with dN significantly greater than dS were discarded. Finally, positions where the derived deleterious allele occurred in another eutherian species were eliminated.

I'm sorry that the best solution for now is to go back to heuristic approaches for prediction. We applied similar heuristics in our 2016 paper, which you can see implemented in this old script: https://github.com/MorrellLAB/Deleterious_Mutations/blob/master/Analysis_Scripts/Count_Deleterious_By_Sample.py (function defined on lines 25-63). I'll look into getting a new logistic regression model with a more modern set of genomes for the next release of BAD_Mutations. My apologies!