Hello,
  I'm am trying to create a codon model using batch commands. The transition matrix is fine, but I am having difficulties harvesting the frequencies. I can "harvestFrequencies(freq, df, 3, 3, 0)", and this returns a column of dimension 64.
The transition matrix is 61x61, the unallowed stop codons account for the difference. How can i harvest the freqs without counting stop codons?

Thank you,
-Nathan Clark
Swanson Lab
U of Washington

Dear Nathan,

In order to get rid of stop codons you'll need to post process the output of HarvestFrequencies - that function will always return 64 (4^3) counts if called with the (3,3,0) argument - HyPhy has no intrinsic knowledge of which codons are stop and which are not. The following example shows how to do this (assuming the use of a standard module to choose a genetic code). The example can also be downloaded from Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login


Code (]
/*1. Get an example file from the web */

GetURL (theData,"Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login;
DataSet ds = ReadFromString(theData);
fprintf (stdout, "Read these data:\n", ds,"\n");

/*2. Include a module to define a genetic code. This will define vector _Genetic_Code (see the comments in chooseGeneticCode.def for details)
      All stop codons (in the chosen code) will have index 10, and all sense codons will have indices ranging from 0-9 and 11-20
*/

pathToModule = HYPHY_BASE_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "chooseGeneticCode.def";
ExecuteCommands ("#include \""+ pathToModule + "\";");


/*3. Filter and get codon frequencies */

DataSetFilter                  filteredData = CreateFilter (ds, 3, "", "", GeneticCodeExclusions);
HarvestFrequencies            (efv, filteredData,3,3,1);

/* make a table of 0's and 1's: 1 for sense codons. See Examples/BatchLanguage/MatrixIndexing.bf for examples of matrix indexing tricks */

senseCodons     = Transpose(_Genetic_Code ["_MATRIX_ELEMENT_VALUE_ != 10"):

); senseCodonFreqs = Transpose(efv[senseCodons]); /* this vector will now hold the frequencies for all non-stop codons */ /*4. Print some diagnostic messages */ chars      = "ACGT"; allCodonIndex   = 0; senseCodonIndex = 0; fprintf (stdout, "\nGathered frequencies:\n"); for (p1=0;p1<4;p1=p1+1) {      for (p2=0;p2<4;p2=p2+1)      {            for (p3=0;p3<4;p3=p3+1)            {                        fprintf (stdout, chars[p1], chars[p2], chars[p3], " : ");                  if (senseCodons[allCodonIndex] == 1)                  {                        fprintf (stdout, senseCodonFreqs[senseCodonIndex], "\n");                                                senseCodonIndex = senseCodonIndex + 1;                  }                  else                  {                        fprintf (stdout, "stop codon\n");                  }                  allCodonIndex = allCodonIndex + 1;            }      } }

Cheers,
Sergei

Thanks!
It worked briliantly.
-Nathan
:)