mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Add HYPHY branch lengths tool.
This commit is contained in:
Daniel Blankenberg
@@ -343,4 +343,7 @@
|
||||
<tool file="mdea/AlphaSubst.xml" />
|
||||
<tool file="stats/kaks.xml" />
|
||||
</section>
|
||||
<section name="HYPHY" id="hyphy">
|
||||
<tool file="hyphy/hyphy_branch_lengths_wrapper.xml" />
|
||||
</section>
|
||||
</toolbox>
|
||||
|
||||
@@ -0,0 +1,53 @@
|
||||
#Dan Blankenberg
|
||||
#takes commandline tree def and input multiple fasta alignment file and runs the branch length ananlysis
|
||||
import os, sys
|
||||
import hyphy_util
|
||||
|
||||
#Retrieve hard coded hyphy path, this will need to be the same across the cluster
|
||||
HYPHY_PATH = hyphy_util.HYPHY_PATH
|
||||
HYPHY_EXECUTABLE = hyphy_util.HYPHY_EXECUTABLE
|
||||
|
||||
#Read command line arguments
|
||||
input_filename = os.path.abspath(sys.argv[1].strip())
|
||||
output_filename = os.path.abspath(sys.argv[2].strip())
|
||||
tree_contents = sys.argv[3].strip()
|
||||
nuc_model = sys.argv[4].strip()
|
||||
base_freq = sys.argv[5].strip()
|
||||
model_options = sys.argv[6].strip()
|
||||
|
||||
#Set up Temporary files for hyphy run
|
||||
#set up tree file
|
||||
tree_filename = hyphy_util.get_filled_temp_filename(tree_contents)
|
||||
|
||||
#Guess if this is a single or multiple FASTA input file
|
||||
found_blank = False
|
||||
is_multiple = False
|
||||
for line in open(input_filename):
|
||||
line = line.strip()
|
||||
if line == "": found_blank = True
|
||||
elif line.startswith(">") and found_blank:
|
||||
is_multiple = True
|
||||
break
|
||||
else: found_blank = False
|
||||
|
||||
#set up BranchLengths file
|
||||
BranchLengths_filename = hyphy_util.get_filled_temp_filename(hyphy_util.BranchLengths)
|
||||
if is_multiple:
|
||||
os.unlink(BranchLengths_filename)
|
||||
BranchLengths_filename = hyphy_util.get_filled_temp_filename(hyphy_util.BranchLengthsMF)
|
||||
print "Multiple Alignment Analyses"
|
||||
else: print "Single Alignment Analyses"
|
||||
|
||||
#setup Config file
|
||||
config_filename = hyphy_util.get_branch_lengths_config_filename(input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename)
|
||||
|
||||
#Run Hyphy
|
||||
hyphy_cmd = "%s BASEPATH=%s USEPATH=/dev/null %s" % (HYPHY_EXECUTABLE, HYPHY_PATH, config_filename)
|
||||
hyphy = os.popen(hyphy_cmd, 'r')
|
||||
#print hyphy.read()
|
||||
hyphy.close()
|
||||
|
||||
#remove temporary files
|
||||
os.unlink(BranchLengths_filename)
|
||||
os.unlink(tree_filename)
|
||||
os.unlink(config_filename)
|
||||
@@ -0,0 +1,93 @@
|
||||
<?xml version="1.0"?>
|
||||
<tool name="Branch Lengths" id="hyphy_branch_lengths_wrapper1">
|
||||
|
||||
<description>Estimation</description>
|
||||
|
||||
<command interpreter="python2.4">hyphy_branch_lengths_wrapper.py $input1 $out_file1 "$tree" "$model" "$base_freq" "Global"</command>
|
||||
|
||||
<inputs>
|
||||
<page>
|
||||
<param format="fasta" name="input1" type="data" label="Fasta file"/>
|
||||
<param name="tree" type="text" label="Tree Definition" size="20" help="For example: ((hg17,panTro1),(mm5,rn3),canFam1)"/>
|
||||
<param name="model" type="select" label="Substitution Model">
|
||||
<option value="000000">F81</option>
|
||||
<option value="010010">HKY85</option>
|
||||
<option value="012345">REV</option>
|
||||
</param>
|
||||
<!-- <param name="model_options" type="select" label="Model Options">
|
||||
<option value="Local">All model parameters are estimated independently for each branch</option>
|
||||
<option value="Global">Model parameters are shared by all branches, branch lengths are estimated independently</option>
|
||||
<option value="Global w/variation">Model parameters are shared by all branches, branch lengths come from a user-chosen distribution, whose parameters are estimated</option>
|
||||
<option value="Global w/variation+HM">Model parameters are shared by all branches, branch lengths come from a user-chosen distribution, whose parameters is estimated; rates at adjacent sites are correlated via a simple Hidden Markov model with an autocorrelation parameter lambda</option>
|
||||
</param> -->
|
||||
<param name="base_freq" type="select" label="Base Frequencies">
|
||||
<option value="Observed">Nucleotide frequencies collected from the data file will be used as equilibrium frequencies</option>
|
||||
<option value="Equal">Equal (.25) frequencies are used as equilibrium frequencies</option>
|
||||
</param>
|
||||
</page>
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data name="out_file1" format="tabular" />
|
||||
</outputs>
|
||||
<test>
|
||||
<param name="input1" value="branchlength_fasta_in.dat"/>
|
||||
<param name="tree" values="((hg17,panTro1),(mm5,rn3),canFam1)"/>
|
||||
<param name="model" values="012345"/>
|
||||
<param name="base_freq" values="Observed"/>
|
||||
<output name="out_file1" file="branchlength_out.dat"/>
|
||||
</test>
|
||||
<help>
|
||||
This tool takes a single or multiple FASTA alignment file and estimates branch lengths using HYPHY_, a maximum likelihood analyses package.
|
||||
|
||||
For the tree definition, you only need to specify the species build names. For example, you could use the tree *((hg17,panTro1),(mm5,rn3),canFam1)*, if your FASTA file looks like this::
|
||||
|
||||
>hg17.chr7(+):26907301-26907310|hg17_0
|
||||
GTGGGAGGT
|
||||
>panTro1.chr6(+):28037319-28037328|panTro1_0
|
||||
GTGGGAGGT
|
||||
>mm5.chr6(+):52104022-52104031|mm5_0
|
||||
GTGGGAGGT
|
||||
>rn3.chr4(+):80734395-80734404|rn3_0
|
||||
GTGGGAGGT
|
||||
>canFam1.chr14(+):42826409-42826418|canFam1_0
|
||||
GTGGGAGGT
|
||||
|
||||
>hg17.chr7(+):26907310-26907326|hg17_1
|
||||
AGTCAGAGTGTCTGAG
|
||||
>panTro1.chr6(+):28037328-28037344|panTro1_1
|
||||
AGTCAGAGTGTCTGAG
|
||||
>mm5.chr6(+):52104031-52104047|mm5_1
|
||||
AGTCAGAGTGTCTGAG
|
||||
>rn3.chr4(+):80734404-80734420|rn3_1
|
||||
AGTCAGAGTATCTGAG
|
||||
>canFam1.chr14(+):42826418-42826434|canFam1_1
|
||||
AGTCAGAGTGTCTGAG
|
||||
|
||||
>hg17.chr7(+):26907326-26907338|hg17_2
|
||||
GTAGAAGACCCC
|
||||
>panTro1.chr6(+):28037344-28037356|panTro1_2
|
||||
GTAGAAGACCCC
|
||||
>mm5.chr6(+):52104047-52104059|mm5_2
|
||||
GTAGACGATGCC
|
||||
>rn3.chr4(+):80734420-80734432|rn3_2
|
||||
GTAGATGATGCG
|
||||
>canFam1.chr14(+):42826434-42826446|canFam1_2
|
||||
GTAGAAGACCCC
|
||||
|
||||
>hg17.chr7(+):26907338-26907654|hg17_3
|
||||
GGGGAAGGAACGCAGGGCGAAGAGCTGGACTTCTCTGAGGAT---TCCTCGGCCTTCTCGT-----CGTTTCCTGG----CGGGGTGGCCGGAGAGATGGGCAAGAGACCCTCCTTCTCACGTTTCTTTTGCTTCATTCGGCGGTTCTGGAACCAGATCTTCACTTGGGTCTCGTTGAGCTGCAGGGATGCAGCGATCTCCACCCTGCGGGCGCGCGTCAGGTACTTGTTGAAGTGGAACTCCTTCTCCAGTTCCGTGAGCTGCTTGGTAGTGAAGTTGGTGCGCACCGCGTTGGGTTGACCCAGGTAGCCGTACTCTCCAACTTTCC
|
||||
>panTro1.chr6(+):28037356-28037672|panTro1_3
|
||||
GGGGAAGGAACGCAGGGCGAAGAGCTGGACTTCTCTGAGGAT---TCCTCGGCCTTCTCGT-----CGTTTCCTGG----CGGGGTGGCCGGAGAGATGGGCAAGAGACCCTCCTTCTCACGTTTCTTTTGCTTCATTCGGCGGTTCTGGAACCAGATCTTCACTTGGGTCTCGTTGAGCTGCAGGGATGCAGCGATCTCCACCCTGCGGGCGCGCGTCAGGTACTTGTTGAAGTGGAACTCCTTCTCCAGTTCCGTGAGCTGCTTGGTAGTGAAGTTGGTGCGCACCGCGTTGGGTTGACCCAGGTAGCCGTACTCTCCAACTTTCC
|
||||
>mm5.chr6(+):52104059-52104375|mm5_3
|
||||
GGAGAAGGGGCACTGGGCGAGGGGCTAGATTTCTCAGATGAT---TCTTCCGTTTTCTCAT-----CGCTGCCAGG----AGGAGTGGCAGGGGAGATGGGCAGGAGCCCCTCCTTCTCACGCTTCTTCTGCTTCATGCGGCGATTCTGGAACCAGATCTTCACCTGGGTCTCATTGAGCTGTAGGGACGCGGCAATCTCCACCCTGCGCGCTCGTGTAAGGTACTTGTTGAAGTGGAACTCCTTCTCCAGCTCTGTGAGCTGCTTGGTGGTGAAATTGGTGCGCACTGCGTTGGGTTGACCCACGTAGCCGTACTCTCCAACTTTCC
|
||||
>rn3.chr4(+):80734432-80734748|rn3_3
|
||||
GGAGAAGGGGCGCTGGGCGAGGAGCTGGATTTCTCAGATGAT---TCTTCAGTTTTCTCAT-----CGCTTCCAGG----AGGGGTGGCGGGTGAAATGGGCAAGAGCCCCTCTTTCTCGCGCTTCTTCTGCTTCATGCGGCGATTCTGGAACCAGATCTTCACCTGGGTCTCATTGAGTTGCAGGGACGCGGCTATCTCCACCCTGCGGGCTCTTGTTAGGTACTTGTTGAAGTGGAACTCCTTCTCCAGCTCTGTGAGCTGCTTGGTGGTGAAGTTGGTGCGCACTGCGTTGGGTTGACCCACGTAGCCATACTCTCCAACTTTCC
|
||||
>canFam1.chr14(+):42826446-42826762|canFam1_3
|
||||
GGAGACGGAATGCAGGGCGAGGAGCTGGATTTCTCTGAAGAT---TCCTCCGCCTTCTCCT-----CACTTCCTGG----CGGGGTGGCAGGGGAGATGGGCAAAAGGCCCTCTTTCTCTCGTTTCTTCTGCTTCATCCGGCGGTTCTGGAACCAGATCTTCACCTGGGTCTCGTTGAGCTGCAGGGATGCTGCGATCTCCACCCTGCGGGCGCGGGTCAGATACTTATTGAAGTGGAACTCCTTTTCCAGCTCGGTGAGCTGCTTGGTGGTGAAGTTGGTACGCACTGCATTCGGTTGACCCACGTAGCCGTACTCTCCAACTTTCC
|
||||
|
||||
|
||||
|
||||
.. _HYPHY: http://www.hyphy.org
|
||||
</help>
|
||||
</tool>
|
||||
|
||||
@@ -0,0 +1,516 @@
|
||||
#Dan Blankenberg
|
||||
#Contains file contents and helper methods for HYPHY configurations
|
||||
import tempfile, os
|
||||
|
||||
def get_filled_temp_filename(contents):
|
||||
fh = tempfile.NamedTemporaryFile('w')
|
||||
filename = fh.name
|
||||
fh.close()
|
||||
fh = open(filename, 'w')
|
||||
fh.write(contents)
|
||||
fh.close()
|
||||
return filename
|
||||
|
||||
#Hard Coded hyphy path, this will need to be the same across the cluster
|
||||
HYPHY_PATH = "/home/universe/linux-i686/HYPHY"
|
||||
HYPHY_EXECUTABLE = os.path.join(HYPHY_PATH,"HYPHY")
|
||||
|
||||
BranchLengthsMF = """
|
||||
VERBOSITY_LEVEL = -1;
|
||||
fscanf (PROMPT_FOR_FILE, "Lines", inLines);
|
||||
|
||||
_linesIn = Columns (inLines);
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
_currentGene = 1;
|
||||
_currentState = 0;
|
||||
geneSeqs = "";
|
||||
geneSeqs * 128;
|
||||
|
||||
for (l=0; l<_linesIn; l=l+1)
|
||||
{
|
||||
if (Abs(inLines[l]) == 0)
|
||||
{
|
||||
if (_currentState == 1)
|
||||
{
|
||||
geneSeqs * 0;
|
||||
DataSet ds = ReadFromString (geneSeqs);
|
||||
_processAGene (_currentGene);
|
||||
geneSeqs * 128;
|
||||
_currentGene = _currentGene + 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (_currentState == 0)
|
||||
{
|
||||
_currentState = 1;
|
||||
}
|
||||
geneSeqs * inLines[l];
|
||||
geneSeqs * "\\n";
|
||||
}
|
||||
}
|
||||
|
||||
if (_currentState == 1)
|
||||
{
|
||||
geneSeqs * 0;
|
||||
if (Abs(geneSeqs))
|
||||
{
|
||||
DataSet ds = ReadFromString (geneSeqs);
|
||||
_processAGene (_currentGene);
|
||||
}
|
||||
}
|
||||
|
||||
fprintf (resultFile,CLOSE_FILE);
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
function _processAGene (_geneID)
|
||||
{
|
||||
DataSetFilter filteredData = CreateFilter (ds,1);
|
||||
if (_currentGene == 1)
|
||||
{
|
||||
SelectTemplateModel (filteredData);
|
||||
|
||||
SetDialogPrompt ("Tree file");
|
||||
fscanf (PROMPT_FOR_FILE, "Tree", givenTree);
|
||||
fscanf (stdin, "String", resultFile);
|
||||
|
||||
/* do sequence to branch map */
|
||||
|
||||
validNames = {};
|
||||
taxonNameMap = {};
|
||||
|
||||
for (k=0; k<TipCount(givenTree); k=k+1)
|
||||
{
|
||||
validNames[TipName(givenTree,k)&&1] = 1;
|
||||
}
|
||||
|
||||
for (k=0; k<BranchCount(givenTree); k=k+1)
|
||||
{
|
||||
thisName = BranchName(givenTree,k);
|
||||
taxonNameMap[thisName&&1] = thisName;
|
||||
}
|
||||
|
||||
storeValidNames = validNames;
|
||||
fprintf (resultFile,CLEAR_FILE,KEEP_OPEN,"Block\\tBranch\\tLength\\tLowerBound\\tUpperBound\\n");
|
||||
}
|
||||
else
|
||||
{
|
||||
validNames = storeValidNames;
|
||||
}
|
||||
|
||||
for (k=0; k<ds.species; k=k+1)
|
||||
{
|
||||
GetString (thisName, ds,k);
|
||||
shortName = (thisName^{{"\\\\..+",""}})&&1;
|
||||
if (validNames[shortName])
|
||||
{
|
||||
taxonNameMap[shortName] = thisName;
|
||||
validNames - (shortName);
|
||||
SetParameter (ds,k,shortName);
|
||||
}
|
||||
else
|
||||
{
|
||||
fprintf (resultFile,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree,"\\n");
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
/* */
|
||||
|
||||
LikelihoodFunction lf = (filteredData,givenTree);
|
||||
|
||||
Optimize (res,lf);
|
||||
|
||||
timer = Time(0)-timer;
|
||||
|
||||
branchNames = BranchName (givenTree,-1);
|
||||
branchLengths = BranchLength (givenTree,-1);
|
||||
|
||||
|
||||
for (k=0; k<Columns(branchNames)-1; k=k+1)
|
||||
{
|
||||
COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
|
||||
COVARIANCE_PRECISION = 0.95;
|
||||
CovarianceMatrix (cmx,lf);
|
||||
if (k==0)
|
||||
{
|
||||
/* compute a scaling factor */
|
||||
ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
|
||||
scaleFactor = BranchLength (givenTree,0);
|
||||
ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
|
||||
}
|
||||
fprintf (resultFile,_geneID,"\\t",taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
|
||||
}
|
||||
|
||||
ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
|
||||
global treeScaler = 1;
|
||||
ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
|
||||
COVARIANCE_PARAMETER = "treeScaler";
|
||||
COVARIANCE_PRECISION = 0.95;
|
||||
CovarianceMatrix (cmx,lf);
|
||||
fprintf (resultFile,_geneID,"\\tTotal Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
|
||||
ClearConstraints (givenTree);
|
||||
return 0;
|
||||
}
|
||||
"""
|
||||
|
||||
BranchLengths = """
|
||||
DataSet ds = ReadDataFile (PROMPT_FOR_FILE);
|
||||
DataSetFilter filteredData = CreateFilter (ds,1);
|
||||
|
||||
SelectTemplateModel (filteredData);
|
||||
|
||||
SetDialogPrompt ("Tree file");
|
||||
fscanf (PROMPT_FOR_FILE, "Tree", givenTree);
|
||||
fscanf (stdin, "String", resultFile);
|
||||
|
||||
/* do sequence to branch map */
|
||||
|
||||
validNames = {};
|
||||
taxonNameMap = {};
|
||||
|
||||
for (k=0; k<TipCount(givenTree); k=k+1)
|
||||
{
|
||||
validNames[TipName(givenTree,k)&&1] = 1;
|
||||
}
|
||||
|
||||
for (k=0; k<BranchCount(givenTree); k=k+1)
|
||||
{
|
||||
thisName = BranchName(givenTree,k);
|
||||
taxonNameMap[thisName&&1] = thisName;
|
||||
}
|
||||
|
||||
for (k=0; k<ds.species; k=k+1)
|
||||
{
|
||||
GetString (thisName, ds,k);
|
||||
shortName = (thisName^{{"\\\\..+",""}})&&1;
|
||||
if (validNames[shortName])
|
||||
{
|
||||
taxonNameMap[shortName] = thisName;
|
||||
validNames - (shortName);
|
||||
SetParameter (ds,k,shortName);
|
||||
}
|
||||
else
|
||||
{
|
||||
fprintf (resultFile,CLEAR_FILE,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree);
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
/* */
|
||||
|
||||
LikelihoodFunction lf = (filteredData,givenTree);
|
||||
|
||||
Optimize (res,lf);
|
||||
|
||||
timer = Time(0)-timer;
|
||||
|
||||
branchNames = BranchName (givenTree,-1);
|
||||
branchLengths = BranchLength (givenTree,-1);
|
||||
|
||||
fprintf (resultFile,CLEAR_FILE,KEEP_OPEN,"Branch\\tLength\\tLowerBound\\tUpperBound\\n");
|
||||
|
||||
for (k=0; k<Columns(branchNames)-1; k=k+1)
|
||||
{
|
||||
COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
|
||||
COVARIANCE_PRECISION = 0.95;
|
||||
CovarianceMatrix (cmx,lf);
|
||||
if (k==0)
|
||||
{
|
||||
/* compute a scaling factor */
|
||||
ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
|
||||
scaleFactor = BranchLength (givenTree,0);
|
||||
ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
|
||||
}
|
||||
fprintf (resultFile,taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
|
||||
}
|
||||
|
||||
ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
|
||||
global treeScaler = 1;
|
||||
ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
|
||||
COVARIANCE_PARAMETER = "treeScaler";
|
||||
COVARIANCE_PRECISION = 0.95;
|
||||
CovarianceMatrix (cmx,lf);
|
||||
ClearConstraints (givenTree);
|
||||
fprintf (resultFile,"Total Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
|
||||
fprintf (resultFile,CLOSE_FILE);
|
||||
"""
|
||||
|
||||
|
||||
SimpleLocalFitter = """
|
||||
VERBOSITY_LEVEL = -1;
|
||||
COUNT_GAPS_IN_FREQUENCIES = 0;
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
function returnResultHeaders (dummy)
|
||||
{
|
||||
_analysisHeaders = {};
|
||||
_analysisHeaders[0] = "GENE";
|
||||
_analysisHeaders[1] = "BP";
|
||||
_analysisHeaders[2] = "S_sites";
|
||||
_analysisHeaders[3] = "NS_sites";
|
||||
_analysisHeaders[4] = "LogL";
|
||||
_analysisHeaders[5] = "AC";
|
||||
_analysisHeaders[6] = "AT";
|
||||
_analysisHeaders[7] = "CG";
|
||||
_analysisHeaders[8] = "CT";
|
||||
_analysisHeaders[9] = "GT";
|
||||
|
||||
for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
|
||||
{
|
||||
branchName = treeBranchNames[_biterator];
|
||||
|
||||
_analysisHeaders [Abs(_analysisHeaders)] = "length("+branchName+")";
|
||||
_analysisHeaders [Abs(_analysisHeaders)] = "dS("+branchName+")";
|
||||
_analysisHeaders [Abs(_analysisHeaders)] = "dN("+branchName+")";
|
||||
_analysisHeaders [Abs(_analysisHeaders)] = "omega("+branchName+")";
|
||||
}
|
||||
|
||||
return _analysisHeaders;
|
||||
}
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
function runAGeneFit (myID)
|
||||
{
|
||||
DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions);
|
||||
|
||||
if (_currentGene==1)
|
||||
{
|
||||
_MG94stdinOverload = {};
|
||||
_MG94stdinOverload ["0"] = "Local";
|
||||
_MG94stdinOverload ["1"] = modelSpecString;
|
||||
|
||||
ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"MG94custom.mdl",
|
||||
_MG94stdinOverload);
|
||||
|
||||
Tree codonTree = treeString;
|
||||
}
|
||||
else
|
||||
{
|
||||
HarvestFrequencies (observedFreq,filteredData,3,1,1);
|
||||
MULTIPLY_BY_FREQS = PopulateModelMatrix ("MG94custom", observedFreq);
|
||||
vectorOfFrequencies = BuildCodonFrequencies (observedFreq);
|
||||
Model MG94customModel = (MG94custom,vectorOfFrequencies,0);
|
||||
|
||||
Tree codonTree = treeString;
|
||||
}
|
||||
|
||||
LikelihoodFunction lf = (filteredData,codonTree);
|
||||
|
||||
Optimize (res,lf);
|
||||
|
||||
_snsAVL = _computeSNSSites ("filteredData", _Genetic_Code, vectorOfFrequencies, 0);
|
||||
_cL = ReturnVectorsOfCodonLengths (ComputeScalingStencils (0), "codonTree");
|
||||
|
||||
|
||||
_returnMe = {};
|
||||
_returnMe ["GENE"] = myID;
|
||||
_returnMe ["LogL"] = res[1][0];
|
||||
_returnMe ["BP"] = _snsAVL ["Sites"];
|
||||
_returnMe ["S_sites"] = _snsAVL ["SSites"];
|
||||
_returnMe ["NS_sites"] = _snsAVL ["NSSites"];
|
||||
_returnMe ["AC"] = AC;
|
||||
_returnMe ["AT"] = AT;
|
||||
_returnMe ["CG"] = CG;
|
||||
_returnMe ["CT"] = CT;
|
||||
_returnMe ["GT"] = GT;
|
||||
|
||||
|
||||
for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
|
||||
{
|
||||
branchName = treeBranchNames[_biterator];
|
||||
|
||||
_returnMe ["length("+branchName+")"] = (_cL["Total"])[_biterator];
|
||||
_returnMe ["dS("+branchName+")"] = (_cL["Syn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["S_sites"]);
|
||||
_returnMe ["dN("+branchName+")"] = (_cL["NonSyn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["NS_sites"]);
|
||||
|
||||
ExecuteCommands ("_lom = _standardizeRatio(codonTree."+treeBranchNames[_biterator]+".nonSynRate,codonTree."+treeBranchNames[_biterator]+".synRate);");
|
||||
_returnMe ["omega("+branchName+")"] = _lom;
|
||||
}
|
||||
|
||||
return _returnMe;
|
||||
}
|
||||
"""
|
||||
|
||||
FastaReader = """
|
||||
fscanf (stdin, "String", _coreAnalysis);
|
||||
fscanf (stdin, "String", _outputDriver);
|
||||
|
||||
ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"chooseGeneticCode.def");
|
||||
ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"dSdNTreeTools.ibf");
|
||||
ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"CodonTools.bf");
|
||||
ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"GrabBag.bf");
|
||||
|
||||
SetDialogPrompt ("Tree file");
|
||||
fscanf (PROMPT_FOR_FILE, "Tree", givenTree);
|
||||
|
||||
treeBranchNames = BranchName (givenTree,-1);
|
||||
treeBranchCount = Columns (treeBranchNames)-1;
|
||||
treeString = Format (givenTree,1,1);
|
||||
|
||||
SetDialogPrompt ("Multiple gene FASTA file");
|
||||
fscanf (PROMPT_FOR_FILE, "Lines", inLines);
|
||||
fscanf (stdin, "String", modelSpecString);
|
||||
fscanf (stdin, "String", _outPath);
|
||||
|
||||
ExecuteAFile (_outputDriver);
|
||||
ExecuteAFile (_coreAnalysis);
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
_linesIn = Columns (inLines);
|
||||
_currentGene = 1;
|
||||
_currentState = 0;
|
||||
/* 0 - waiting for a non-empty line */
|
||||
/* 1 - reading files */
|
||||
|
||||
geneSeqs = "";
|
||||
geneSeqs * 0;
|
||||
|
||||
_prepareFileOutput (_outPath);
|
||||
|
||||
for (l=0; l<_linesIn; l=l+1)
|
||||
{
|
||||
if (Abs(inLines[l]) == 0)
|
||||
{
|
||||
if (_currentState == 1)
|
||||
{
|
||||
geneSeqs * 0;
|
||||
DataSet ds = ReadFromString (geneSeqs);
|
||||
_processAGene (ds.species == treeBranchCount,_currentGene);
|
||||
geneSeqs * 128;
|
||||
_currentGene = _currentGene + 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (_currentState == 0)
|
||||
{
|
||||
_currentState = 1;
|
||||
}
|
||||
geneSeqs * inLines[l];
|
||||
geneSeqs * "\\n";
|
||||
}
|
||||
}
|
||||
|
||||
if (_currentState == 1)
|
||||
{
|
||||
geneSeqs * 0;
|
||||
DataSet ds = ReadFromString (geneSeqs);
|
||||
_processAGene (ds.species == treeBranchCount,_currentGene);
|
||||
}
|
||||
|
||||
_finishFileOutput (0);
|
||||
"""
|
||||
|
||||
TabWriter = """
|
||||
/*---------------------------------------------------------*/
|
||||
function _prepareFileOutput (_outPath)
|
||||
{
|
||||
_outputFilePath = _outPath;
|
||||
|
||||
_returnHeaders = returnResultHeaders(0);
|
||||
|
||||
fprintf (_outputFilePath, CLEAR_FILE, KEEP_OPEN, _returnHeaders[0]);
|
||||
for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
|
||||
{
|
||||
fprintf (_outputFilePath,"\\t",_returnHeaders[_biterator]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
fprintf (_outputFilePath,"\\n");
|
||||
return 0;
|
||||
}
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
|
||||
function _processAGene (valid, _geneID)
|
||||
{
|
||||
if (valid)
|
||||
{
|
||||
returnValue = runAGeneFit (_geneID);
|
||||
fprintf (_outputFilePath, returnValue[_returnHeaders[0]]);
|
||||
for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
|
||||
{
|
||||
fprintf (_outputFilePath,"\\t",returnValue[_returnHeaders[_biterator]]);
|
||||
}
|
||||
fprintf (_outputFilePath, "\\n");
|
||||
}
|
||||
else
|
||||
{
|
||||
fprintf (_outputFilePath,
|
||||
_geneID, ", Incorrect number of sequences\\n");
|
||||
}
|
||||
_currentState = 0;
|
||||
return 0;
|
||||
}
|
||||
|
||||
/*---------------------------------------------------------*/
|
||||
function _finishFileOutput (dummy)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
"""
|
||||
|
||||
def get_dN_dS_config_filename(SimpleLocalFitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename ):
|
||||
contents = """
|
||||
_genomeScreenOptions = {};
|
||||
|
||||
/* all paths are either absolute or relative
|
||||
to the DATA READER */
|
||||
|
||||
_genomeScreenOptions ["0"] = "%s";
|
||||
/* which analysis to run on each gene; */
|
||||
_genomeScreenOptions ["1"] = "%s";
|
||||
/* what output to produce; */
|
||||
|
||||
_genomeScreenOptions ["2"] = "%s";
|
||||
/* genetic code */
|
||||
_genomeScreenOptions ["3"] = "%s";
|
||||
/* tree file */
|
||||
_genomeScreenOptions ["4"] = "%s";
|
||||
/* alignment file */
|
||||
_genomeScreenOptions ["5"] = "%s";
|
||||
/* nucleotide bias string; can define any of the 203 models */
|
||||
_genomeScreenOptions ["6"] = "%s";
|
||||
/* output csv file */
|
||||
|
||||
ExecuteAFile ("%s", _genomeScreenOptions);
|
||||
""" % (SimpleLocalFitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename )
|
||||
return get_filled_temp_filename(contents)
|
||||
|
||||
|
||||
def get_branch_lengths_config_filename(input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename):
|
||||
contents = """
|
||||
_genomeScreenOptions = {};
|
||||
|
||||
/* all paths are either absolute or relative
|
||||
to the NucDataBranchLengths.bf */
|
||||
|
||||
_genomeScreenOptions ["0"] = "%s";
|
||||
/* the file to analyze; */
|
||||
_genomeScreenOptions ["1"] = "CUSTOM";
|
||||
/* use an arbitrary nucleotide model */
|
||||
_genomeScreenOptions ["2"] = "%s";
|
||||
/* which model to use */
|
||||
_genomeScreenOptions ["3"] = "%s";
|
||||
/* model options */
|
||||
_genomeScreenOptions ["4"] = "Estimated";
|
||||
/* rate parameters */
|
||||
_genomeScreenOptions ["5"] = "%s";
|
||||
/* base frequencies */
|
||||
_genomeScreenOptions ["6"] = "%s";
|
||||
/* the tree to use; */
|
||||
_genomeScreenOptions ["7"] = "%s";
|
||||
/* write .csv output to; */
|
||||
|
||||
ExecuteAFile ("%s", _genomeScreenOptions);
|
||||
""" % (input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename)
|
||||
return get_filled_temp_filename(contents)
|
||||
|
||||
Reference in New Issue
Block a user