File: balibase_test.org

package info (click to toggle)
kalign 1%3A3.4.0-1
links: PTS, VCS
area: main
in suites: forky, sid, trixie
size: 2,160 kB
sloc: ansic: 10,480; sh: 65; makefile: 57; cpp: 30
file content (1308 lines) | stat: -rw-r--r-- 48,555 bytes
parent folder | download | duplicates (2)
#+TITLE:  Test kalign on BAliBASE
#+AUTHOR: Timo Lassmann 
#+EMAIL:  timo.lassmann@telethonkids.org.au
#+DATE:   2019-04-23
#+LATEX_CLASS: report
#+OPTIONS:  toc:nil
#+OPTIONS: H:4
#+LATEX_CMD: pdflatex
* Introduction 
  Let's run kalign on the core BAliBASE test alignments. 
  
* Method 

** Step 1: install balibase 

   #+BEGIN_SRC sh
     cd 
     mkdir -p data
     cd data
     if [ ! -f BAliBASE_R1-5.tar.gz ]; then
         wget https://www.lbgi.fr/balibase/BalibaseDownload/BAliBASE_R1-5.tar.gz

     fi
     tar -zxvf  BAliBASE_R1-5.tar.gz
   #+END_SRC

   #+RESULTS:

   BAliBASE comes with a C program to score alignments. Unfortunately, something is wrong with the bundled expat library. To fix this install expat system and remove the links to the bundled library by editing the makefile: 

   #+BEGIN_EXAMPLE makefile
   CC	= cc
   CFLAGS  = -c -O #-I$(EXPAT_INC)
   LFLAGS	= -O -lm -lexpat #-L$(EXPAT_LIB) -lexpat
   EXPAT_LIB	= expat-1.95.2/lib
   EXPAT_INC	= expat-1.95.2/include
   #+END_EXAMPLE

   Also, by default bali_score will use all columns for the TC score. In the original clustal omega paper only core blocks are used to calculate the TC scores supercite:sievers-2014-fast-scalab. To change this modify =bali_score.c= line 104 from: 

   #+BEGIN_EXAMPLE C 
   method='C' 
   #+END_EXAMPLE

   to: 

   #+BEGIN_EXAMPLE C 
   method='B' 
   #+END_EXAMPLE

   Then compile:

   #+BEGIN_SRC sh
     cd ~/data/bb3_release/bali_score_src
     make 
   #+END_SRC

   #+RESULTS:
   | cc | -c | -O         | readxml.c    |        |        |              |    |     |         |
   | cc | -c | -O         | init.c       |        |        |              |    |     |         |
   | cc | -c | -O         | util.c       |        |        |              |    |     |         |
   | cc | -c | -O         | bali_score.c |        |        |              |    |     |         |
   | cc | -o | bali_score | readxml.o    | init.o | util.o | bali_score.o | -O | -lm | -lexpat |

** Step 2: run tests

   #+BEGIN_SRC sh :session onesh
     cd ~/data/bb3_release
     cd .. 

     mkdir -p bb3_release_tmp_aln
     cd bb3_release_tmp_aln
     KALIGNOUTDIR=$PWD 
     echo $KALIGNOUTDIR


   #+END_SRC

   #+RESULTS:
   |                                                                             |
   | sh-4.4$ sh-4.4$ sh-4.4$ sh-4.4$ sh-4.4$ /home/user/data/bb3_release_tmp_aln |
   #+BEGIN_SRC sh :session onesh :results raw
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -alnp param_25.txt %s -o %s/%s_%skalign25.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_kalign.sh


     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -alnp param_50.txt  %s -o %s/%s_%skalign50.msf\n", $1,outdir,a[n-2],a[n-1] }' >> run_kalign.sh

     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -alnp param_75.txt  %s -o %s/%s_%skalign75.msf\n", $1,outdir,a[n-2],a[n-1] }' >> run_kalign.sh

   #+END_SRC 

   #+BEGIN_SRC sh :session onesh :results raw 
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -alnp testparam  %s -o %s/%s_%skalign.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_kalign.sh

     find ~/data/bb3_release -name "*.tfa"  | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); printf "clustalo --dealign -i  %s --outfmt=msf -o %s/%s_%sclustalo.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_clustalo.sh

     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/);
          printf "kalign %s -set 0 -o %s/%s_%skalign_1.msf\n", $1,outdir,a[n-2],a[n-1] ;
          printf "kalign %s -set 1 -o %s/%s_%skalign_2.msf\n", $1,outdir,a[n-2],a[n-1] ;
          printf "kalign %s -set 2 -o %s/%s_%skalign_3.msf\n", $1,outdir,a[n-2],a[n-1] ;
          printf "kalign %s -set 3 -o %s/%s_%skalign_4.msf\n", $1,outdir,a[n-2],a[n-1] ;
          }' > run_kalign_cmsa.sh


     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/);
          printf "cmsa -a  %s/%s_%skalign_*.msf -out %s/%s_%skalign_cmsa.msf -f msf \n",outdir,a[n-2],a[n-1],outdir,a[n-2],a[n-1] ;
          }' > run_cmsa.sh

     chmod 755 run_kalign_cmsa.sh 


     chmod 755 run_kalign.sh 
     chmod 755 run_clustalo.sh

   #+END_SRC

   #+RESULTS:

   #+BEGIN_SRC sh :session onesh :results raw 

     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores.sh

     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_1.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_set1.sh
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_2.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_set2.sh
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_3.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_set3.sh
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_4.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_set4.sh
     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_5.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_set5.sh


     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%skalign_cmsa.msf\n", $1,outdir,a[n-2],a[n-1] }' > alignment_scores_cmsa.sh

     find ~/data/bb3_release -name "*.xml" | awk -v outdir="$KALIGNOUTDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%sclustalo.msf\n", $1,outdir,a[n-2],a[n-1] }' >> alignment_scores.sh


     chmod 755 alignment_scores.sh
   #+END_SRC


   #+RESULTS:

   run the tests 

   #+BEGIN_SRC sh :session onesh
     parallel --jobs 5 < ./run_kalign.sh 
     parallel --jobs 5 < ./run_clustalo.sh 
     ./alignment_scores.sh | grep auto > scores2.csv

   #+END_SRC
   
** Step 3: plot scores 

*** Check for dependencies
   
    Here I define the list of libraries I'll be using. 
   



    #+NAME: liblist
    #+BEGIN_SRC R -n :exports code :results none
      libraries <- c("tidyverse","ggplot2","stringi","cowplot")

    #+END_SRC
   
    Script to test if libraries are present.
    #+BEGIN_SRC R -n :tangle test_for_libraries.R :shebang #!/usr/bin/env Rscript :noweb yes :exports code :results none
      <<liblist>>
      Sys.info()["nodename"]
      for(library in libraries) 
      { 
        f = is.element(library, installed.packages()[,1])
        print(paste("Library",library, "is installed?", f))
        if(!f)
        {
          message("Missing library:",library )
          quit(status=1)
        }
      }
      quit(status=0)
    #+END_SRC
   
   
    #+BEGIN_SRC sh -n :results output :exports both
      ./test_for_libraries.R
    #+END_SRC

    #+RESULTS:
    : nodename 
    :   "work" 
    : [1] "Library tidyverse is installed? TRUE"
    : [1] "Library ggplot2 is installed? TRUE"
    : [1] "Library stringi is installed? TRUE"
    : [1] "Library cowplot is installed? TRUE"

    install.packages("tidyverse")
    Code block to load the libraries in R code.

    #+NAME: Rlibraries
    #+BEGIN_SRC R -n :exports code :results none :noweb yes
      <<liblist>>
      lapply(libraries, FUN = function(X) {
        do.call("library", list(X)) 
      })

    #+END_SRC

*** Plotting functions 
    #+BEGIN_SRC R :session  one :results none :export none :noweb yes 
      <<Rlibraries>> 
      readBaliscores <-function(file,name){
        mat  <- read.table(file);
        colnames(mat) <- c("Type","Name","SP","TC")
        mat$Type <- name
        mat$Name <- sub(".*/" ,"", mat$Name)
        mat$Name <- sub("[A-Z,_]*[.]{1}[A-Z]*$" ,"", mat$Name,ignore.case = TRUE)
        mat <- as.tibble(mat)
        x = str_split(mat$Name, "_", n = Inf, simplify = TRUE)
        mat$Group <- x[,1] 
        return(mat)

      }


    #+END_SRC


   

    #+BEGIN_SRC R :session one :results output graphics :file BalibaseSP_scores.jpeg :exports both :width 160 :height 80 :noweb yes 

      mat <- readBaliscores("scores_kalign_old.csv","kalign old");
      mat <- rbind(mat,readBaliscores("scores_kalign_new.csv","kalign new"));
      mat <- rbind(mat,readBaliscores("scores_kalign_pbil2.csv","kalign pbil 2"));
      mat <- rbind(mat,readBaliscores("scores_kalign_16seed.csv","kalign 16 seed"));
      mat <- rbind(mat,readBaliscores("scores_kalign_test.csv","rel"))
      mat <- rbind(mat,readBaliscores("scores_kalign_optimized.csv","Opt"))
      ##      mat <- rbind(mat,readBaliscores("scores_kalign_3.csv","kalign 3"));

      ## mata
      ## <- rbind(mat,readBaliscores("scores_kalign_bibpm.csv","bibpm"));
      ## mat 
      ## <- rbind(mat,readBaliscores("scores_kalign_bibpm_zero.csv","bibpm_zero"));


      p1 <- ggplot(mat, aes(Group, SP))
      p1 <- p1 + geom_boxplot(aes(colour = Type))

      means <- aggregate(SP ~  Type, mat, median)
      means$SP <- round(means$SP,digits = 4)
      p2 <- ggplot(mat, aes(Type, SP))
      p2 <- p2 + geom_boxplot(aes(colour = Type))
      p2 <- p2 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
      p2 <- p2 + geom_text(data = means, aes(label = SP, y = SP + 0.08))

      p3 <- ggplot(mat, aes(Group, TC))
      p3 <- p3 + geom_boxplot(aes(colour = Type))

      means <- aggregate(TC ~  Type, mat, median)
      means$TC <- round(means$TC,digits = 4)

      p4 <- ggplot(mat, aes(Type, TC))
      p4 <- p4 + geom_boxplot(aes(colour = Type))
      p4 <- p4 + geom_boxplot(aes(colour = Type))
      p4 <- p4 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
      p4 <- p4 + geom_text(data = means, aes(label = TC, y = TnC + 0.08))



      p  = plot_grid(p1,p2,p3,p4, labels=c("SP", "TC"), ncol = 2, nrow= 2)




    #+END_SRC

    #+RESULTS:
    [[file:BalibaseSP_scores.jpeg]]

* Benchmarking 

  Programs to test: 

  #+NAME: Benchprograms
  | Name         | Download                                                                     |
  |--------------+------------------------------------------------------------------------------|
  | Muscle       | https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux32.tar.gz |
  | Clustalomega | http://clustal.org/omega/clustalo-1.2.4-Ubuntu-x86_64                        |
  | q-score      | https://www.drive5.com/qscore/qscore_src.tar.gz                              |
  | Kalign2      | http://msa.sbc.su.se/downloads/kalign/current.tar.gz                         |
  
  Since the clustal papers provide a great overview of alignment programs I think these two should suffice. Note I am not comparing kalign against consistency based methods as these are a) all more accurate and b) slower. Kalign is targeted at the accurate and fast niche. 

  Benchmarks: 
  #+NAME: Benchmarkdata 
  | Name          | Download                                                          |
  |---------------+-------------------------------------------------------------------|
  | Balibase1-5   | http://www.lbgi.fr/balibase/BalibaseDownload/BAliBASE_R1-5.tar.gz |
  | Balifam       | http://clustal.org/omega/bali3fam-26.tar.gz                       |
  | Bralibaseset1 | http://projects.binf.ku.dk/pgardner/bralibase/data-set1.tar.gz    |
  | Bralibaseset2 | http://projects.binf.ku.dk/pgardner/bralibase/data-set2.tar.gz    |
  | HomFam        | http://www.clustal.org/omega/homfam-20110613-25.tar.gz            |
  | QuanTest2     | http://bioinf.ucd.ie/quantest2.tar                                |
  


** Step one: create test directory and download all data

   #+BEGIN_SRC bash :exports both :results none 
     cd 
     mkdir -p kalignbenchmark
     cd  kalignbenchmark
     mkdir -p data
     mkdir -p programs 
     mkdir -p scratch
   #+END_SRC

*** Download benchmark datasets 
     
    #+BEGIN_SRC bash -n :results raw :exports both :var tbl=Benchmarkdata :colnames yes
      cd
      cd ~/kalignbenchmark/data
      for idx in ${!tbl[*]}; do
          #echo  ${tbl[$idx]} 
          if [[ ${tbl[$idx]} =~ ^http* ]];
          then
              echo "wget ${tbl[$idx]}"
              wget ${tbl[$idx]}
          fi
      done
    #+END_SRC



    inter all data sets 
    #+BEGIN_SRC bash  :results none :exports code 
      cd
      cd ~/kalignbenchmark/data
      mkdir -p homfam 
      mv homfam-20110613-25.tar.gz homfam 
      cd homfam 
      tar -zxvf homfam-20110613-25.tar.gz

      cd ~/kalignbenchmark/data
      for filename in *.tar.gz; do
          tar -zxvf  $filename
      done
      tar -xvf quantest2.tar
    #+END_SRC

*** Compile baliscore program and place in program directory 

    #+BEGIN_SRC bash

      cd ~/kalignbenchmark/data/bb3_release/bali_score_src

      gcc *.c -lm -lexpat -o bali_score 
      cp bali_score ~/kalignbenchmark/programs

    #+END_SRC

    #+RESULTS:

*** Download  programs 

    #+BEGIN_SRC bash -n :results none :exports none :var tbl=Benchprograms :colnames yes
      cd
      cd ~/kalignbenchmark/programs
      for idx in ${!tbl[*]}; do
          #echo  ${tbl[$idx]} 
          if [[ ${tbl[$idx]} =~ ^http* ]];
          then
              echo "wget ${tbl[$idx]}"
              wget ${tbl[$idx]}
          fi
      done

      chmod 755 clustalo-1.2.4-Ubuntu-x86_64
      tar -zxvf muscle3.8.31_i86linux32.tar.gz 
      tar -zxvf qscore_src.tar.gz

      mkdir -p kalign2_src
      cd kalign2_src
      mv ../current.tar.gz .
      tar -zxvf current.tar.gz
      ./configure
      make 
      cp kalign ../kalign2
    #+END_SRC




*** Compile q-score 

    NOTE: this is manual! 

    Add:

    #+BEGIN_EXAMPLE C
    #include <limits.h> 
    #+END_EXAMPLE

    to qscore header: =qscore.h=  

    then: 
    #+BEGIN_EXAMPLE bash
    make
    #+END_EXAMPLE



** Step 2: test  if all programs run 

   Let's make sure individual programs run on all the test sets. 

   small code chunk to set path:

   #+NAME: path
   #+BEGIN_SRC bash :results none :exports code :noweb yes
     export PATH=~/kalignbenchmark/programs:$PATH
   #+END_SRC

*** muscle 

    Can I run muscle? 
    
    #+BEGIN_SRC bash :results value :exports code :noweb yes 
      <<path>>

      muscle3.8.31_i86linux32 -version
    #+END_SRC

    #+RESULTS:
    : MUSCLE v3.8.31 by Robert C. Edgar
    
    Let's try on a balibase case: 

    #+BEGIN_SRC bash :results raw :exports code :noweb yes
      <<path>> 
      IN=~/kalignbenchmark/data/bb3_release/RV11/BB11002.tfa
      REF=~/kalignbenchmark/data/bb3_release/RV11/BB11002.msf
      OUT=~/kalignbenchmark/scratch/BB11002_muscle.msf 
      muscle3.8.31_i86linux32 -msf -in $IN -out $OUT
      bali_score $REF $OUT

    #+END_SRC



    Let's try on a balifam 

    #+BEGIN_SRC bash :results raw :exports both :noweb yes
      <<path>> 

      IN=~/kalignbenchmark/data/BB20021-package.vie
      ADD=~/kalignbenchmark/data/bb3_release/RV20/BB20021.tfa
      REF=~/kalignbenchmark/data/bb3_release/RV20/BB20021.msf

      OUT=~/kalignbenchmark/scratch/BB20021_muscle.fa

      cd ~/kalignbenchmark/scratch

      cat $IN $ADD > test.fa
      time muscle3.8.31_i86linux32 -maxiters 2 -in test.fa -out $OUT
      kalign $REF -r -f fasta -o ref.fa
      qscore -test $OUT -ref ref.fa

    #+END_SRC

    #+RESULTS:
    [2019-06-11 20:55:42] :     LOG : kalign -r -f fasta -o ref.fa /home/user/kalignbenchmark/data/bb3_release/RV20/BB20021.msf 
    [2019-06-11 20:55:42] :     LOG : Detected protein sequences.
    [2019-06-11 20:55:42] :     LOG : Detected: 53 sequences.
    [2019-06-11 20:55:42] :     LOG : Output in fasta format





*** clustalo 
    Can I run muscle? 
    
    #+BEGIN_SRC bash :results value :exports code :noweb yes 
      <<path>>

      clustalo-1.2.4-Ubuntu-x86_64 --version 
    #+END_SRC
    Let's try on a balibase case: 

    #+BEGIN_SRC bash :results raw :exports code :noweb yes
      <<path>> 
      IN=~/kalignbenchmark/data/bb3_release/RV11/BB11002.tfa
      REF=~/kalignbenchmark/data/bb3_release/RV11/BB11002.msf
      OUT=~/kalignbenchmark/scratch/BB11002_clustalo.msf 

      clustalo-1.2.4-Ubuntu-x86_64 --outfmt=msf  --in $IN  --out $OUT

      bali_score $REF $OUT 

    #+END_SRC

    #+RESULTS:
    Using GCG format reference alignment

    Comparing test alignment in /home/user/kalignbenchmark/scratch/BB11002_clustalo.msf
    with reference alignment in /home/user/kalignbenchmark/data/bb3_release/RV11/BB11002.msf

	  SP score= 0.483

	  TC score= 0.000
    auto /home/user/kalignbenchmark/scratch/BB11002_clustalo.msf 0.483 0.000


** Step 3: scripts to run benchmarks  

   #+NAME: scratchdir
   #+BEGIN_SRC bash :results none :exports code :noweb yes
     SCRATCHDIR="~/kalignbenchmark/scratch"
   #+END_SRC


*** Balibase 

    #+BEGIN_SRC bash -n :tangle gen_run_balibase.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 
      <<scratchdir>>
    #+END_SRC

    #+BEGIN_SRC bash -n :tangle gen_run_balibase.sh :noweb yes :exports code :results none :noweb yes 

      DIR=`pwd`
      function usage()
      {
          printf "This script will generate scripts to run kalign, clustal omega and muscle on the balibase MSA benchmark data set.\n\n" ;
          printf "usage: $0\n\n" ;
          exit 1;
      }

      while getopts h  opt
      do
          case ${opt} in
              h) usage;;
              ,*) usage;;
          esac
      done

      printf "Generating kalign run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.tfa" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -i %s -f msf -o %s/%s_%s_kalign.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_kalign.sh
      chmod 755 run_kalign.sh

      printf "Generating muscle run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.tfa" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "muscle3.8.31_i86linux32 -msf -in %s -out  %s/%s_%s_muscle.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_muscle.sh
      chmod 755 run_muscle.sh

      printf "Generating clustalo run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.tfa" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "clustalo-1.2.4-Ubuntu-x86_64 --outfmt=msf  --in %s  --out %s/%s_%s_clustalo.msf\n", $1,outdir,a[n-2],a[n-1] }' > run_clustalo.sh
      chmod 755 run_clustalo.sh


      printf "Generating kalign scoring run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.xml" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/programs/bali_score %s %s/%s_%s_kalign.msf | grep auto\n", $1,outdir,a[n-2],a[n-1] }' > score_bb3_kalign.sh
      chmod 755  score_bb3_kalign.sh

      printf "Generating muscle scoring run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.xml" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/programs/bali_score %s %s/%s_%s_muscle.msf | grep auto\n", $1,outdir,a[n-2],a[n-1] }' > score_bb3_muscle.sh
      chmod 755  score_bb3_muscle.sh

      printf "Generating clustalo scoring run script\n";
      find ~/kalignbenchmark/data/bb3_release -name "*.xml" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/programs/bali_score %s %s/%s_%s_clustalo.msf | grep auto\n", $1,outdir,a[n-2],a[n-1] }' > score_bb3_clustalo.sh
      chmod 755  score_bb3_clustalo.sh


    #+END_SRC


**** Kick off jobs 


     #+BEGIN_SRC bash -n :tangle run_balibase.sh :noweb yes :exports code :results none :noweb yes :shebang #!/usr/bin/env bash
       <<path>>
       ./run_muscle.sh
       ./run_kalign.sh
       ./run_clustalo.sh 
     #+END_SRC

**** score alignments 

     #+BEGIN_SRC bash 
       ./score_bb3_kalign.sh > scores_kalign.csv 
       ./score_bb3_muscle.sh > scores_muscle.csv 
       ./score_bb3_clustalo.sh > scores_clustalo.csv 
     #+END_SRC

     #+RESULTS:

**** Plot 

     #+BEGIN_SRC R :session  one :results none :export none :noweb yes 
       <<Rlibraries>> 
       readBaliscores <-function(file,name){
         mat  <- read.table(file);
         colnames(mat) <- c("Type","Name","SP","TC")
         mat$Type <- name
         mat$Name <- sub(".*/" ,"", mat$Name)
         mat$Name <- sub("[A-Z,_]*[.]{1}[A-Z]*$" ,"", mat$Name,ignore.case = TRUE)
         mat <- as.tibble(mat)
         x = str_split(mat$Name, "_", n = Inf, simplify = TRUE)
         mat$Group <- x[,1] 
         return(mat)

       }


     #+END_SRC


     #+BEGIN_SRC R :session one :results output graphics :file Balibase_scores.jpeg :exports both :noweb yes 

       mat <- readBaliscores("scores_kalign.csv","kalign");
       mat <- rbind(mat,readBaliscores("scores_muscle.csv","muscle"));
       mat <- rbind(mat,readBaliscores("scores_clustalo.csv","clustalo"));
       ##      mat <- rbind(mat,readBaliscores("scores_kalign_3.csv","kalign 3"));

       ## mata
       ## <- rbind(mat,readBaliscores("scores_kalign_bibpm.csv","bibpm"));
       ## mat 
       ## <- rbind(mat,readBaliscores("scores_kalign_bibpm_zero.csv","bibpm_zero"));


       p1 <- ggplot(mat, aes(Group, SP))
       p1 <- p1 + geom_boxplot(aes(colour = Type))

       means <- aggregate(SP ~  Type, mat, median)
       means$SP <- round(means$SP,digits = 4)
       p2 <- ggplot(mat, aes(Type, SP))
       p2 <- p2 + geom_boxplot(aes(colour = Type))
       p2 <- p2 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
       p2 <- p2 + geom_text(data = means, aes(label = SP, y = SP + 0.08))

       p3 <- ggplot(mat, aes(Group, TC))
       p3 <- p3 + geom_boxplot(aes(colour = Type))

       means <- aggregate(TC ~  Type, mat, median)
       means$TC <- round(means$TC,digits = 4)

       p4 <- ggplot(mat, aes(Type, TC))
       p4 <- p4 + geom_boxplot(aes(colour = Type))
       p4 <- p4 + geom_boxplot(aes(colour = Type))
       p4 <- p4 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
       p4 <- p4 + geom_text(data = means, aes(label = TC, y = TC + 0.08))

       p1 <- p1 + theme(legend.position="none")
       p2 <- p2 + theme(legend.position="none")
       p3 <- p3 + theme(legend.position="none")
       p4 <- p4 + theme(legend.position="none")

       p  = plot_grid(p1,p2,p3,p4, labels=c("SP", "TC"), ncol = 2, nrow= 2)

       ggsave("Balibase_scores.jpeg", p,width = 16,  height = 12,  dpi = 300, units = c( "cm"))
       p1 <- ggplot(mat, aes(Group, SP))
       p1 <- p1 + geom_boxplot(aes(colour = Type))
       p1 <- p1 + scale_color_discrete(name = "Program")      
       p1 <- p1 + theme(legend.position = c(0.25, 0.15),legend.text=element_text(size=16),  legend.direction = "horizontal")
       p1 <- p1 + xlab(label="")
       p1 <- p1 + theme(plot.margin = margin(t = 0, b = -0.4, unit  = "cm"))# unit(c(1,2,3,4), "cm")) # t = 0, r = 0, b = 0, l = 0 ("left", "right", "bottom", "top")
       ggsave("Figure1.jpeg",p1,width = 16,  height = 12,  dpi = 300, units = c( "cm"))
     #+END_SRC

     #+RESULTS:
     [[file:Balibase_scores.jpeg]]


*** Bralibase 


    #+BEGIN_SRC bash -n :tangle gen_run_bralibase.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 
      <<scratchdir>>
    #+END_SRC

    #+BEGIN_SRC bash -n :tangle gen_run_bralibase.sh :noweb yes :exports code :results none :noweb yes 

      DIR=`pwd`
      function usage()
      {
          printf "This script will generate scripts to run kalign, clustal omega and muscle on the balibase MSA benchmark data set.\n\n" ;
          printf "usage: $0\n\n" ;
          exit 1;
      }

      while getopts h  opt
      do
          case ${opt} in
              h) usage;;
              ,*) usage;;
          esac
      done

      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2  -name "*_test.fa" -exec rm -rf {} \;
      printf "Generate commands to reformat the reference alignments\n\n";

      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*.fa"  -o -name '*.fasta' | grep structural |  awk -v outdir="$SCRATCHDIR" '{printf "kalign %s --rename -reformat -f msf -o %s.msf\n", $1,$1 }' > run_kalign_reformat.sh
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*.fa"  -o -name '*.fasta' | grep unaligned |  awk -v outdir="$SCRATCHDIR" '{printf "kalign %s --rename -reformat -f fasta -o %s_test.fa\n", $1,$1 }' >> run_kalign_reformat.sh
      chmod 755 run_kalign_reformat.sh
      ./run_kalign_reformat.sh

      printf "Generating kalign run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2  -name "*_test.fa" | grep unaligned |  awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign %s -f msf -o %s/%s_%s_%s_kalign.msf\n", $1,outdir,a[n-4],a[n-3],a[n-2] }' > run_kalign_bralibase.sh
      chmod 755 run_kalign_bralibase.sh

      printf "Generating muscle run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*_test.fa" -o -name '*.fasta' | grep unaligned | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "muscle3.8.31_i86linux32 -msf -in %s -out  %s/%s_%s_%s_muscle.msf\n", $1,outdir,a[n-4],a[n-3],a[n-2] }' > run_muscle_bralibase.sh
      chmod 755 run_muscle_bralibase.sh

      printf "Generating clustalo run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*_test.fa" -o -name '*.fasta' | grep unaligned | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "clustalo-1.2.4-Ubuntu-x86_64 --outfmt=msf  --in %s --force --out %s/%s_%s_%s_clustalo.msf\n", $1,outdir,a[n-4],a[n-3],a[n-2] }' > run_clustalo_bralibase.sh
      chmod 755 run_clustalo_bralibase.sh

      printf "Generating kalign scoring run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*.msf" | grep structural | awk -v outdir="$SCRATCHDIR" '{name=$1; gsub(/structural/, "unaligned"); n=split ($1,a,/[\/,.]/); printf "~/kalignbenchmark/programs/bali_score  %s %s/%s_%s_%s_kalign.msf | grep auto\n", name,outdir,a[n-4],a[n-3],a[n-2] }' > score_kalign_bralibase.sh
      chmod 755  score_kalign_bralibase.sh

      printf "Generating muscle scoring run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*.msf" | grep structural | awk -v outdir="$SCRATCHDIR" '{name=$1; gsub(/structural/, "unaligned"); n=split ($1,a,/[\/,.]/); printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%s_%s_muscle.msf | grep auto\n", name,outdir,a[n-4],a[n-3],a[n-2] }' > score_muscle_bralibase.sh
      chmod 755  score_muscle_bralibase.sh

      printf "Generating clustalo scoring run script\n";
      find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*.msf" | grep structural | awk -v outdir="$SCRATCHDIR" '{name=$1; gsub(/structural/, "unaligned"); n=split ($1,a,/[\/,.]/); printf "~/data/bb3_release/bali_score_src/bali_score %s %s/%s_%s_%s_clustalo.msf | grep auto\n", name,outdir,a[n-4],a[n-3],a[n-2] }' > score_clustalo_bralibase.sh
      chmod 755  score_clustalo_bralibase.sh


    #+END_SRC

**** Kick off jobs 


     #+BEGIN_SRC bash -n :tangle run_bralibase.sh :noweb yes :exports code :results none :noweb yes :shebang  #!/usr/bin/env bash
       <<path>>

       ./run_kalign_bralibase.sh   
       ./run_muscle_bralibase.sh
       ./run_clustalo_bralibase.sh       


     #+END_SRC




**** score alignments 


     #+BEGIN_SRC bash 
       ./score_kalign_bralibase.sh > scores_kalign_bralibase.csv 
       ./score_muscle_bralibase.sh > scores_muscle_bralibase.csv 
       ./score_clustalo_bralibase.sh > scores_clustalo_bralibase.csv 

     #+END_SRC

     #+RESULTS:

**** Plot 

     #+BEGIN_SRC R :session  one :results none :export none :noweb yes 
       <<Rlibraries>> 
       readBaliscores <-function(file,name){
         mat  <- read.table(file);
         colnames(mat) <- c("Type","Name","SP","TC")
         mat$Type <- name
         mat$Name <- sub(".*/" ,"", mat$Name)
         mat$Name <- sub("[A-Z,_]*[.]{1}[A-Z]*$" ,"", mat$Name,ignore.case = TRUE)
         mat <- as.tibble(mat)
         x = str_split(mat$Name, "_", n = Inf, simplify = TRUE)
         mat$Group <- x[,1] 
         return(mat)

       }


     #+END_SRC


     #+BEGIN_SRC R :session one :results output graphics :file Bralibase_scores.jpeg :exports both :noweb yes 

       mat <- readBaliscores("scores_kalign_bralibase.csv","kalign");
       mat <- rbind(mat,readBaliscores("scores_muscle_bralibase.csv","muscle"));
       mat <- rbind(mat,readBaliscores("scores_clustalo_bralibase.csv","clustalo"));
       ##      mat <- rbind(mat,readBaliscores("scores_kalign_3.csv","kalign 3"));

       ## mata
       ## <- rbind(mat,readBaliscores("scores_kalign_bibpm.csv","bibpm"));
       ## mat 
       ## <- rbind(mat,readBaliscores("scores_kalign_bibpm_zero.csv","bibpm_zero"));


       p1 <- ggplot(mat, aes(Group, SP))
       p1 <- p1 + geom_boxplot(aes(colour = Type))

       means <- aggregate(SP ~  Type, mat, median)
       means$SP <- round(means$SP,digits = 4)
       p2 <- ggplot(mat, aes(Type, SP))
       p2 <- p2 + geom_boxplot(aes(colour = Type))
       p2 <- p2 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
       p2 <- p2 + geom_text(data = means, aes(label = SP, y = SP + 0.08))

       p3 <- ggplot(mat, aes(Group, TC))
       p3 <- p3 + geom_boxplot(aes(colour = Type)) 


       means <- aggregate(TC ~  Type, mat, median)
       means$TC <- round(means$TC,digits = 4)


       p4 <- ggplot(mat, aes(Type, TC))
       p4 <- p4 + geom_boxplot(aes(colour = Type))
       p4 <- p4 + geom_boxplot(aes(colour = Type))
       p4 <- p4 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
       p4 <- p4 + geom_text(data = means, aes(label = TC, y = TC + 0.08))

       p1 <- p1 + theme(legend.position="none")
       p2 <- p2 + theme(legend.position="none")
       p3 <- p3 + theme(legend.position="none")
       p4 <- p4 + theme(legend.position="none")


       p  = plot_grid(p1,p2,p3,p4, labels=c("SP", "TC"), ncol = 2, nrow= 2)

       ggsave("Bralibase_scores.jpeg", p,width = 16,  height = 12,  dpi = 300, units = c( "cm"))

       p1 <- ggplot(mat, aes(Group, SP))
       p1 <- p1 + geom_boxplot(aes(colour = Type))
       p1 <- p1 + scale_color_discrete(name = "Program")      
       p1 <- p1 + theme(legend.position = c(0.25, 0.15),legend.text=element_text(size=16),  legend.direction = "horizontal")
       p1 <- p1 + xlab(label="")
       p1 <- p1 + theme(plot.margin = margin(t = 0, b = -0.4, unit  = "cm"))# unit(c(1,2,3,4), "cm")) # t = 0, r = 0, b = 0, l = 0 ("left", "right", "bottom", "top")
       ggsave("Figure2.jpeg",p1,width = 16,  height = 12,  dpi = 300, units = c( "cm"))
     #+END_SRC

     #+RESULTS:
     [[file:Bralibase_scores.jpeg]]

*** Homfam


****  Create test alignments 
     create script to merge =*_ref.vie= and =*_test-only.vie= alignments.
     #+BEGIN_SRC bash -n :tangle create_homfam_tests.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 
       <<scratchdir>>
     #+END_SRC


     #+BEGIN_SRC bash -n :tangle create_homfam_tests.sh :noweb yes :exports code :results none :noweb yes


       DIR=`pwd`
       function usage()
       {
           printf "This script will merge *_ref.vie and *_test-only.vie alignments and convert them into msf format.\n\n" ;
           printf "usage: $0\n\n" ;
           exit 1;
       }

       while getopts h  opt
       do
           case ${opt} in
               h) usage;;
               ,*) usage;;
           esac
       done
       find ~/kalignbenchmark/data/homfam  -name "*_test-only.vie" | awk '{test_aln=$1;ref=$1 ;gsub("_test-only.vie", "_ref.vie", ref);gsub("_test-only.vie", "_testaln.fa", test_aln); printf "kalign %s %s --reformat --rename -f fasta -o %s\n",ref, $1, test_aln }' > homfam_reformat.sh
       find ~/kalignbenchmark/data/homfam  -name "*_ref.vie" | awk '{ref=$1;gsub("_ref.vie", "_ref.msf", ref);printf "kalign %s --reformat --rename  -f msf -o %s\n", $1,ref}' >> homfam_reformat.sh

       chmod 755 homfam_reformat.sh
     #+END_SRC     


**** Run alignment programs on Homfam 


     #+BEGIN_SRC bash -n :tangle gen_run_homfam.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 
       <<scratchdir>>
     #+END_SRC

     #+BEGIN_SRC bash -n :tangle gen_run_homfam.sh :noweb yes :exports code :results none :noweb yes

       DIR=`pwd`
       function usage()
       {
           printf "This script will generate scripts to run kalign, clustal omega and muscle on the balibase MSA benchmark data set.\n\n" ;
           printf "usage: $0\n\n" ;
           exit 1;
       }

       while getopts h  opt
       do
           case ${opt} in
               h) usage;;
               ,*) usage;;
           esac
       done



       printf "Generating kalign run script\n"
       

printf "export PATH=~/kalignbenchmark/programs:$PATH\n\n" > run_homfam.sh
       find ~/kalignbenchmark/data/homfam  -name "*_testaln.fa" | awk -v outdir="$SCRATCHDIR" '{
       ref_aln=$1;gsub("_testaln.fa", "_ref.msf",ref_aln);
       printf "timescorealn -test %s -ref %s -program kalign --scratch ~/kalignbenchmark/scratch -out scores_homfam.csv

\n", $1,ref_aln;
       printf "timescorealn -test %s -ref %s -program muscle --scratch ~/kalignbenchmark/scratch -out scores_homfam.csv\n", $1,ref_aln;
       printf "timescorealn -test %s -ref %s -program clustalo --scratch ~/kalignbenchmark/scratch -out scores_homfam.csv\n", $1,ref_aln;
       }' >> run_homfam.sh

       chmod 755 run_homfam.sh




     #+END_SRC

     #+BEGIN_EXAMPLE bash -n :tangle run_homfam.sh :noweb yes :exports code :results none :noweb yes :shebang #!/usr/bin/env bash
       <<path>>
       ./run_homfam.sh
     #+END_EXAMPLE






     /home/user/kalignbenchmark/data/homfam

     printf "Generating clustalo run script\n";
     find ~/kalignbenchmark/data/data-set1  ~/kalignbenchmark/data/data-set2 -name "*_test.fa" -o -name '*.fasta' | grep unaligned | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "clustalo-1.2.4-Ubuntu-x86_64 --outfmt=msf  --in %s --force --out %s/%s_%s_%s_clustalo.msf\n", $1,outdir,a[n-4],a[n-3],a[n-2] }' > run_clustalo_bralibase.sh
     chmod 755 run_clustalo_bralibase.sh
**** plot scores

     #+BEGIN_SRC R :session one :results output graphics :file homfam_runtime.jpeg :exports both :noweb yes
       <<Rlibraries>> 
       mat  <- read.csv("scores_homfam.csv");
       p <- ggplot(mat, aes(NUMSEQ,Time))
       p <- p + geom_point(aes(colour = Program,size = AVGLEN))
       p <- p + scale_x_continuous(trans = 'log10', breaks=c(100,3000,10000,100000) ,labels = scales::comma)
       p <- p + scale_y_continuous(trans = 'log10',breaks=c(0.01, 0.1, 1 , 10,100,1000,10000), labels=scales::comma )

       ggsave("homfam_runtime.jpeg",p,width = 18,  height = 12,  dpi = 300, units = c( "cm"))
     #+END_SRC

     #+RESULTS:
     [[file:homfam_runtime.jpeg]]
     q[[file:homfam_runtime.jpeg]]





     #+BEGIN_SRC R :session one :results output graphics :file homfam_scores.jpeg :exports both :noweb yes


       p <- ggplot(mat, aes(Program, SP))
       p <- p + geom_boxplot(aes(colour = Program))
       ggsave("homfam_scores.jpeg",p,width = 16,  height = 12,  dpi = 300, units = c( "cm"))
     #+END_SRC

     #+RESULTS:
     [[file:homfam_scores.jpeg]]

*** QuantTest 2 

**** script to generate run and score commands

    #+BEGIN_SRC bash -n :tangle gen_run_quanttest2.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 
      <<scratchdir>>

    #+END_SRC

    #+BEGIN_SRC bash -n :tangle gen_run_quanttest2.sh :noweb yes :exports code :results none :noweb yes 

      DIR=`pwd`
      function usage()
      {
          printf "This script will generate scripts to run kalign, clustal omega and muscle on the quanttest2 MSA benchmark data set.\n\n" ;
          printf "usage: $0\n\n" ;
          exit 1;
      }

      while getopts h  opt
      do
          case ${opt} in
              h) usage;;
              ,*) usage;;
          esac
      done

      printf "Generating kalign run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "kalign -i %s -o %s/%s_kalign.afa\n", $1,outdir,a[n-1] }' > run_kalign_quanttest.sh
      chmod 755 run_kalign_quanttest.sh

      printf "Generating muscle run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "muscle3.8.31_i86linux32 -in %s -out  %s/%s_muscle.afa\n", $1,outdir,a[n-1] }' > run_muscle_quanttest.sh
      chmod 755 run_muscle_quanttest.sh

      printf "Generating clustalo run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); printf "clustalo-1.2.4-Ubuntu-x86_64 --threads 8 --outfmt=a2m  --in %s  --out %s/%s_clustalo.afa\n", $1,outdir,a[n-1] }' > run_clustalo_quanttest.sh
      chmod 755 run_clustalo_quanttest.sh


      printf "Generating kalign scoring run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/data/QuanTest2/quantest2.py test@test  %s/%s_kalign.afa ~/kalignbenchmark/data/QuanTest2/SS/%s.ss  | grep kalignbenchmark \n", outdir,a[n-1],a[n-1] }' > score_kalign_quanttest.sh
      chmod 755  score_kalign_quanttest.sh

      printf "Generating muscle scoring run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/data/QuanTest2/quantest2.py test@test  %s/%s_muscle.afa ~/kalignbenchmark/data/QuanTest2/SS/%s.ss  | grep kalignbenchmark \n", outdir,a[n-1],a[n-1] }' > score_muscle_quanttest.sh
      chmod 755  score_muscle_quanttest.sh


      printf "Generating clustalo scoring run script\n";
      find ~/kalignbenchmark/data/QuanTest2/Test -name "*.vie" | awk -v outdir="$SCRATCHDIR" '{n=split ($1,a,/[\/,.]/); ;printf "~/kalignbenchmark/data/QuanTest2/quantest2.py test@test  %s/%s_clustalo.afa ~/kalignbenchmark/data/QuanTest2/SS/%s.ss  | grep kalignbenchmark \n", outdir,a[n-1],a[n-1] }' > score_clustalo_quanttest.sh
      chmod 755  score_clustalo_quanttest.sh


    #+END_SRC

**** Run jobs and score alignments 

     #+BEGIN_SRC bash -n :tangle run_quanttest.sh :noweb yes :exports code :results none :noweb yes :shebang #!/usr/bin/env bash
       <<path>>

       ./run_muscle_quanttest.sh
       ./run_kalign_quanttest.sh
       ./run_clustalo_quanttest.sh

       ./score_kalign_quanttest.sh > scores_quanttest_kalign.csv
       ./score_muscle_quanttest.sh > scores_quanttest_muscle.csv
       ./score_clustalo_quanttest.sh > scores_quanttest_clustalo.csv
     #+END_SRC

**** Plot

  #+BEGIN_SRC R :session  one :results none :export none :noweb yes 
    <<Rlibraries>> 
    readquanttestscores <-function(file,name){
      mat  <- read.table(file);
      colnames(mat) <- c("Name","Score")
      mat$Type <- name
      mat$Name <- sub(".*/" ,"", mat$Name)
      mat$Name <- sub("[_]{1}[A-Z]*[.]{1}[A-Z]*$" ,"", mat$Name,ignore.case = TRUE)
      mat <- as_tibble(mat)
      
      return(mat)

    }


  #+END_SRC


  #+BEGIN_SRC R :session one :results output graphics :file Quantest.jpeg :exports both :noweb yes 


    mat <- readquanttestscores("scores_quanttest_clustalo.csv","clustalo")
    mat <- rbind(mat, readquanttestscores("scores_quanttest_kalign.csv","kalign"));
    mat <- rbind(mat, readquanttestscores("scores_quanttest_muscle.csv","muscle"));
#mat <- mat  %>% mutate(Score = as.dbl(Score))      

       
       p4 <- ggplot(mat, aes(Type, Score))
       p4 <- p4 + geom_boxplot(aes(colour = Type))

       means <- aggregate(Score ~  Type, mat, median)
       means$Score <- round(means$Score,digits = 4)
       
       p4 <- p4 + stat_summary(fun.y=mean, colour="darkred", geom="point", shape=18, size=3,show.legend = FALSE) 
       p4 <- p4 + geom_text(data = means, aes(label = Score, y = Score + 0.08))

       ggsave("Quantest.jpeg",p4,width = 16,  height = 12,  dpi = 300, units = c( "cm"))


t = t.test(filter(mat,Type == "kalign")$Score,filter(mat,Type == "muscle")$Score, paired = FALSE, alternative = "greater")

      t = t.test(filter(mat,Type == "kalign")$Score,filter(mat,Type == "clustalo")$Score, paired = FALSE, alternative = "greater")

  #+END_SRC

  #+RESULTS:
  [[file:Quantest.jpeg]]

*** Rose
    :PROPERTIES:
    :ORDERED:  t
    :END:


    #+BEGIN_SRC bash -n :tangle create_rose_tests.sh :shebang #!/usr/bin/env bash :noweb yes :exports code :results none :noweb yes 


      DIR=`pwd`
      function usage()
      {
          printf "This script will merge *_ref.vie and *_test-only.vie alignments and convert them into msf format.\n\n" ;
          printf "usage: $0\n\n" ;
          exit 1;
      }

      while getopts h  opt
      do
          case ${opt} in
              h) usage;;
              ,*) usage;;
          esac
      done


      numseq=(100 1000 10000 100000 1000000 10000000)
      for item in ${numseq[*]}
      do
          printf "# kalignbenchmark sample 1\n" > rose.in;
          printf "\n" >> rose.in;
          printf "%%include protein-defaults\n"  >> rose.in;
          printf "\n"  >> rose.in;
          printf "StdOut= False\n" >> rose.in
          printf "OutputFilebase = \"rosetest%d\"\n" $item  >> rose.in;
          printf "AlignmentFormat= \"FASTA\"\n" >> rose.in;
          printf "SequenceLen = 500\n" >> rose.in;
          printf "SequenceNum = %d\n" $item >> rose.in;
          printf "ChooseFromLeaves = False\n" >> rose.in;
          printf "TreeWithSequences = True\n" >> rose.in;
          rose rose.in


      done




    #+END_SRC



    # rose sample 1

    %include protein-defaults

    OutputFilename = "TIMO"
    OutputFilebase = "timotest"
    AlignmentFormat= "FASTA"
    SequenceLen = 500
    SequenceNum = 1000
    ChooseFromLeaves = False
    TreeWithSequences = True
* Paper figure 


  #+BEGIN_SRC R :session  one :results none :export none :noweb yes 
    <<Rlibraries>> 
    readBaliscores <-function(file,name){
      mat  <- read.table(file);
      colnames(mat) <- c("Type","Name","SP","TC")
      mat$Type <- name
      mat$Name <- sub(".*/" ,"", mat$Name)
      mat$Name <- sub("[A-Z,_]*[.]{1}[A-Z]*$" ,"", mat$Name,ignore.case = TRUE)
      mat <- as.tibble(mat)
      x = str_split(mat$Name, "_", n = Inf, simplify = TRUE)
      mat$Group <- x[,1] 
      return(mat)

    }


  #+END_SRC


  #+BEGIN_SRC R :session one :results output graphics :file Paper_figure.jpeg :exports both :noweb yes 

    library(ggplot2)
    library(ggsignif)
    my_comparisons = list( c("clustalo", "kalign"), c("kalign", "muscle"), c("clustalo", "muscle") )

    mat <- readBaliscores("scores_kalign.csv","kalign");
    mat <- rbind(mat,readBaliscores("scores_muscle.csv","muscle"));
    mat <- rbind(mat,readBaliscores("scores_clustalo.csv","clustalo"));
    p1 <- ggplot(mat, aes(Group, SP))
    p1 <- p1 + geom_boxplot(aes(colour = Type))
    p1 <- p1 + scale_color_discrete(name = "Program")      
    p1 <- p1 + theme(legend.position = c(0.25, 0.15),legend.text=element_text(size=16),  legend.direction = "horizontal")
    p1 <- p1 + xlab(label="")
    ##  p1 <- p1 +  geom_signif(comparisons = my_comparisons ,map_signif_level=TRUE)
    mat <- readBaliscores("scores_kalign_bralibase.csv","kalign");
    mat <- rbind(mat,readBaliscores("scores_muscle_bralibase.csv","muscle"));
    mat <- rbind(mat,readBaliscores("scores_clustalo_bralibase.csv","clustalo"));
    ##      mat <- rbind(mat,readBaliscores("scores_kalign_3.csv","kalign 3"));

    ## mata
    ## <- rbind(mat,readBaliscores("scores_kalign_bibpm.csv","bibpm"));
    ## mat 
    ## <- rbind(mat,readBaliscores("scores_kalign_bibpm_zero.csv","bibpm_zero"));


    p2 <- ggplot(mat, aes(Group, SP))
    p2 <- p2 + geom_boxplot(aes(colour = Type))
    p2 <- p2 + scale_color_discrete(name = "Program")      
    p2 <- p2 + theme(legend.position = c(0.25, 0.15),legend.text=element_text(size=16),  legend.direction = "horizontal")
    p2 <- p2 + xlab(label="")
                                            #      p2 <- p2 + theme(plot.margin = margin(t = 0, b = -0.2, unit  = "cm"))# unit(c(1,2,3,4), "cm")) # t = 0, r 

    mat  <- read.csv("scores_homfam.csv");
    p3 <- ggplot(mat, aes(NUMSEQ,Time))
    p3 <- p3 + geom_point(aes(colour = Program,size = AVGLEN))
    p3 <- p3 + scale_x_continuous(trans = 'log10', breaks=c(100,3000,10000,100000) ,labels = scales::comma)
    p3 <- p3 + scale_y_continuous(trans = 'log10',breaks=c(0.01, 0.1, 1 , 10,100,1000,10000),labels=scales::comma )

    p3 <- p3 + xlab("Number of sequences") 
    p3 <- p3 + ylab("Time (s)") 
    p = plot_grid(p1, p2, p3, labels = c('A','B','C'), ncol = 1, nrow = 3, rel_heights = c(0.7,0.7,1.0))
    ggsave("Paper_figure.jpeg",p,width = 18,  height = 24,  dpi = 300, units = c( "cm")) 
  #+END_SRC

  #+RESULTS:
  [[file:Paper_figure.jpeg]]



  Test for statistical significant differences: 

  #+BEGIN_SRC R :session one 
    assess_diff_alignment_scores <-function(mat, groupname,tests){
      f <- filter(mat, Group == groupname)  %>% select(-TC,-Group)  %>% spread(Type, SP) %>% column_to_rownames("Name") 

      t = t.test(f$kalign,f$clustalo, paired = TRUE, alternative = "greater")
      if(t$p.value * tests <= 0.05){
        message("sig diff:  kalign, clustalo ", groupname," ", t$p.value);
        message("Clustalo ", mean(f$clustalo));
        message("Kalign ", mean(f$kalign));
        message("Muscle ", mean(f$muscle));
      }

      t = t.test(f$kalign,f$muscle, paired = TRUE, alternative = "greater")
      if(t$p.value * tests  <= 0.05){
        message("sig diff: kalign, muscle ", groupname," ", t$p.value);
        message("Clustalo ", mean(f$clustalo));
        message("Kalign ", mean(f$kalign));
        message("Muscle ", mean(f$muscle));
      }

      
    }

    mat <- readBaliscores("scores_kalign.csv","kalign");
    mat <- rbind(mat,readBaliscores("scores_muscle.csv","muscle"));
    mat <- rbind(mat,readBaliscores("scores_clustalo.csv","clustalo"));

    g <- unique(mat$Group) 
    tests = length(g) * length(unique(mat$Type))
    for (i in g){
      assess_diff_alignment_scores(mat,i,tests);
    }

    mat <- readBaliscores("scores_kalign_bralibase.csv","kalign");
    mat <- rbind(mat,readBaliscores("scores_muscle_bralibase.csv","muscle"));
    mat <- rbind(mat,readBaliscores("scores_clustalo_bralibase.csv","clustalo"));

    g <- unique(mat$Group) 
    tests = length(g) * length(unique(mat$Type))
    for (i in g){
      assess_diff_alignment_scores(mat,i,tests);
    }

  #+END_SRC

  #+RESULTS:
* Slurm scripts 



  #+BEGIN_SRC bash -n :tangle srun_kalign.sh :shebang #!/usr/bin/env bash :exports code :results none :noweb yes
    #SBATCH -N 1 # number of nodes
    #SBATCH -n 1 # number of cores
    #SBATCH -t 0-10:00 # time (D-HH:MM)
    #SBATCH -o slurm.%N.%j.out # STDOUT
    #SBATCH -e slurm.%N.%j.err # STDERR
    DIR=`pwd`
    OUTPUT=
    INPUT=


    function usage()
    {

        printf "Slurm script to run kalign.\n\n" ;
        printf "usage: $0 -i input -o output\n\n" ;
        exit 1;
    }

    while getopts i:o:  opt
    do
        case ${opt} in
            i) INPUT=${OPTARG};;
            o) OUTPUT=${OPTARG};;            
            ,*) usage;;
        esac
    done
    if [ "${INPUT}" == "" ]; then usage; fi
    if [ "${OUTPUT}" == "" ]; then usage; fi

    kalign -i $INPUT -f msf -o $OUTPUT
  #+END_SRC