Changed the location of normalization procedures for independent and epistatic models!

computePred.R

@@ -75,8 +75,13 @@ nbGaps = N[1] - apply(nbSeqs,2,sum)
 # output the conservation values
 dat=rbind(trace,KL=res[[1]][[2]],SE=res[[1]][[1]],gap=nbGaps/N[1],KL60=res[[2]][[2]],SE60=res[[2]][[1]],KL80=res[[3]][[2]],SE80=res[[3]][[1]])
 write.table(dat,paste0(prot,"_conservation.txt"))
+
 # compute log-odd ratios between mutated and wt sequence counts
 predInd = computePredNbSeqs(wt,nbSeqs)
+# output the sequence counts log-odd ratios
+write.table(predInd,paste0(prot,"_pred_evolInd.txt"))
+
+# Do the normalization for the independent component here!
 rownames(predInd)=aa
 if(simple){
   normPredInd = normalizePredWithNbSeqsZero(predInd,trace,wt)
@@ -85,14 +90,20 @@ if(simple){
   normPredInd = normalizePredWithNbSeqsZeroSelMult(predInd, trace, wt, list(pos,subsaa))
   names(normPredInd)=rawMut
 }
-# output the sequence counts log-odd ratios
-write.table(predInd,paste0(prot,"_pred_evolInd.txt"))
+
 # output the predicted mutational effects based on sequence counts (conservation at the bottom)
 write.table(normPredInd,paste0(prot,"_normPred_evolInd.txt"))
 print("done")
 
 print("running global epistatic model...")
 pred=computePredSimple(ali,distTrace,wt,5)
+
+# output the evolutionary distances between the query and the closest variants
+evolDist = pred/sum(trace^2)
+evolDist[is.na(evolDist)] = 1
+write.table(evolDist,paste0(prot,"_pred_evolEpi.txt"))
+
+# Do the normalization for the epistatic component here!
 rownames(pred)=aa
 if(simple){
   normPred=normalizePred(pred, trace, wt)
@@ -101,10 +112,6 @@ if(simple){
   normPred=normalizePredSelMult(pred, trace, wt, list(pos,subsaa))
   names(normPred)=rawMut
 }
-# output the evolutionary distances between the query and the closest variants
-evolDist = pred/sum(trace^2)
-evolDist[is.na(evolDist)] = 1
-write.table(evolDist,paste0(prot,"_pred_evolEpi.txt"))
 # output the predicted mutational effects based on evolutionary distance (conservation at the bottom)
 write.table(normPred,paste0(prot,"_normPred_evolEpi.txt"))
 print("done")

example/example-gemme-script.sh

+#!/bin/bash
+
+if [ "$1" == "clean" ]
+then
+    echo "Cleaning the folder."
+    rm -rf BLAT*
+    rm -f default.conf caracTest.dat
+else
+    echo "Running GEMME with a user-provided alignment file."
+    python $GEMME_PATH/gemme.py aliBLAT.fasta -r input -f aliBLAT.fasta
+fi