diff baseline/Baseline_Functions.r @ 0:8a5a2abbb870 draft default tip

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/baseline/Baseline_Functions.r	Mon Aug 29 05:36:10 2016 -0400
@@ -0,0 +1,2287 @@
+#########################################################################################
+# License Agreement
+# 
+# THIS WORK IS PROVIDED UNDER THE TERMS OF THIS CREATIVE COMMONS PUBLIC LICENSE 
+# ("CCPL" OR "LICENSE"). THE WORK IS PROTECTED BY COPYRIGHT AND/OR OTHER 
+# APPLICABLE LAW. ANY USE OF THE WORK OTHER THAN AS AUTHORIZED UNDER THIS LICENSE 
+# OR COPYRIGHT LAW IS PROHIBITED.
+# 
+# BY EXERCISING ANY RIGHTS TO THE WORK PROVIDED HERE, YOU ACCEPT AND AGREE TO BE 
+# BOUND BY THE TERMS OF THIS LICENSE. TO THE EXTENT THIS LICENSE MAY BE CONSIDERED 
+# TO BE A CONTRACT, THE LICENSOR GRANTS YOU THE RIGHTS CONTAINED HERE IN 
+# CONSIDERATION OF YOUR ACCEPTANCE OF SUCH TERMS AND CONDITIONS.
+#
+# BASELIne: Bayesian Estimation of Antigen-Driven Selection in Immunoglobulin Sequences
+# Coded by: Mohamed Uduman & Gur Yaari
+# Copyright 2012 Kleinstein Lab
+# Version: 1.3 (01/23/2014)
+#########################################################################################
+
+# Global variables  
+  
+  FILTER_BY_MUTATIONS = 1000
+
+  # Nucleotides
+  NUCLEOTIDES = c("A","C","G","T")
+  
+  # Amino Acids
+  AMINO_ACIDS <- c("F", "F", "L", "L", "S", "S", "S", "S", "Y", "Y", "*", "*", "C", "C", "*", "W", "L", "L", "L", "L", "P", "P", "P", "P", "H", "H", "Q", "Q", "R", "R", "R", "R", "I", "I", "I", "M", "T", "T", "T", "T", "N", "N", "K", "K", "S", "S", "R", "R", "V", "V", "V", "V", "A", "A", "A", "A", "D", "D", "E", "E", "G", "G", "G", "G")
+  names(AMINO_ACIDS) <- c("TTT", "TTC", "TTA", "TTG", "TCT", "TCC", "TCA", "TCG", "TAT", "TAC", "TAA", "TAG", "TGT", "TGC", "TGA", "TGG", "CTT", "CTC", "CTA", "CTG", "CCT", "CCC", "CCA", "CCG", "CAT", "CAC", "CAA", "CAG", "CGT", "CGC", "CGA", "CGG", "ATT", "ATC", "ATA", "ATG", "ACT", "ACC", "ACA", "ACG", "AAT", "AAC", "AAA", "AAG", "AGT", "AGC", "AGA", "AGG", "GTT", "GTC", "GTA", "GTG", "GCT", "GCC", "GCA", "GCG", "GAT", "GAC", "GAA", "GAG", "GGT", "GGC", "GGA", "GGG")
+  names(AMINO_ACIDS) <- names(AMINO_ACIDS)
+
+  #Amino Acid Traits
+  #"*" "A" "C" "D" "E" "F" "G" "H" "I" "K" "L" "M" "N" "P" "Q" "R" "S" "T" "V" "W" "Y"
+  #B = "Hydrophobic/Burried"  N = "Intermediate/Neutral"  S="Hydrophilic/Surface") 
+  TRAITS_AMINO_ACIDS_CHOTHIA98 <- c("*","N","B","S","S","B","N","N","B","S","B","B","S","N","S","S","N","N","B","B","N")
+  names(TRAITS_AMINO_ACIDS_CHOTHIA98) <- sort(unique(AMINO_ACIDS))
+  TRAITS_AMINO_ACIDS <- array(NA,21)
+  
+  # Codon Table
+  CODON_TABLE <- as.data.frame(matrix(NA,ncol=64,nrow=12))
+
+  # Substitution Model: Smith DS et al. 1996
+  substitution_Literature_Mouse <- matrix(c(0, 0.156222928, 0.601501588, 0.242275484, 0.172506739, 0, 0.241239892, 0.586253369, 0.54636291, 0.255795364, 0, 0.197841727, 0.290240811, 0.467680608, 0.24207858, 0),nrow=4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
+  substitution_Flu_Human <- matrix(c(0,0.2795596,0.5026927,0.2177477,0.1693210,0,0.3264723,0.5042067,0.4983549,0.3328321,0,0.1688130,0.2021079,0.4696077,0.3282844,0),4,4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
+  substitution_Flu25_Human <- matrix(c(0,0.2580641,0.5163685,0.2255674,0.1541125,0,0.3210224,0.5248651,0.5239281,0.3101292,0,0.1659427,0.1997207,0.4579444,0.3423350,0),4,4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
+  load("FiveS_Substitution.RData")
+
+  # Mutability Models: Shapiro GS et al. 2002
+  triMutability_Literature_Human <- matrix(c(0.24, 1.2, 0.96, 0.43, 2.14, 2, 1.11, 1.9, 0.85, 1.83, 2.36, 1.31, 0.82, 0.52, 0.89, 1.33, 1.4, 0.82, 1.83, 0.73, 1.83, 1.62, 1.53, 0.57, 0.92, 0.42, 0.42, 1.47, 3.44, 2.58, 1.18, 0.47, 0.39, 1.12, 1.8, 0.68, 0.47, 2.19, 2.35, 2.19, 1.05, 1.84, 1.26, 0.28, 0.98, 2.37, 0.66, 1.58, 0.67, 0.92, 1.76, 0.83, 0.97, 0.56, 0.75, 0.62, 2.26, 0.62, 0.74, 1.11, 1.16, 0.61, 0.88, 0.67, 0.37, 0.07, 1.08, 0.46, 0.31, 0.94, 0.62, 0.57, 0.29, NA, 1.44, 0.46, 0.69, 0.57, 0.24, 0.37, 1.1, 0.99, 1.39, 0.6, 2.26, 1.24, 1.36, 0.52, 0.33, 0.26, 1.25, 0.37, 0.58, 1.03, 1.2, 0.34, 0.49, 0.33, 2.62, 0.16, 0.4, 0.16, 0.35, 0.75, 1.85, 0.94, 1.61, 0.85, 2.09, 1.39, 0.3, 0.52, 1.33, 0.29, 0.51, 0.26, 0.51, 3.83, 2.01, 0.71, 0.58, 0.62, 1.07, 0.28, 1.2, 0.74, 0.25, 0.59, 1.09, 0.91, 1.36, 0.45, 2.89, 1.27, 3.7, 0.69, 0.28, 0.41, 1.17, 0.56, 0.93, 3.41, 1, 1, NA, 5.9, 0.74, 2.51, 2.24, 2.24, 1.95, 3.32, 2.34, 1.3, 2.3, 1, 0.66, 0.73, 0.93, 0.41, 0.65, 0.89, 0.65, 0.32, NA, 0.43, 0.85, 0.43, 0.31, 0.31, 0.23, 0.29, 0.57, 0.71, 0.48, 0.44, 0.76, 0.51, 1.7, 0.85, 0.74, 2.23, 2.08, 1.16, 0.51, 0.51, 1, 0.5, NA, NA, 0.71, 2.14), nrow=64,byrow=T)
+  triMutability_Literature_Mouse <- matrix(c(1.31, 1.35, 1.42, 1.18, 2.02, 2.02, 1.02, 1.61, 1.99, 1.42, 2.01, 1.03, 2.02, 0.97, 0.53, 0.71, 1.19, 0.83, 0.96, 0.96, 0, 1.7, 2.22, 0.59, 1.24, 1.07, 0.51, 1.68, 3.36, 3.36, 1.14, 0.29, 0.33, 0.9, 1.11, 0.63, 1.08, 2.07, 2.27, 1.74, 0.22, 1.19, 2.37, 1.15, 1.15, 1.56, 0.81, 0.34, 0.87, 0.79, 2.13, 0.49, 0.85, 0.97, 0.36, 0.82, 0.66, 0.63, 1.15, 0.94, 0.85, 0.25, 0.93, 1.19, 0.4, 0.2, 0.44, 0.44, 0.88, 1.06, 0.77, 0.39, 0, 0, 0, 0, 0, 0, 0.43, 0.43, 0.86, 0.59, 0.59, 0, 1.18, 0.86, 2.9, 1.66, 0.4, 0.2, 1.54, 0.43, 0.69, 1.71, 0.68, 0.55, 0.91, 0.7, 1.71, 0.09, 0.27, 0.63, 0.2, 0.45, 1.01, 1.63, 0.96, 1.48, 2.18, 1.2, 1.31, 0.66, 2.13, 0.49, 0, 0, 0, 2.97, 2.8, 0.79, 0.4, 0.5, 0.4, 0.11, 1.68, 0.42, 0.13, 0.44, 0.93, 0.71, 1.11, 1.19, 2.71, 1.08, 3.43, 0.4, 0.67, 0.47, 1.02, 0.14, 1.56, 1.98, 0.53, 0.33, 0.63, 2.06, 1.77, 1.46, 3.74, 2.93, 2.1, 2.18, 0.78, 0.73, 2.93, 0.63, 0.57, 0.17, 0.85, 0.52, 0.31, 0.31, 0, 0, 0.51, 0.29, 0.83, 0.54, 0.28, 0.47, 0.9, 0.99, 1.24, 2.47, 0.73, 0.23, 1.13, 0.24, 2.12, 0.24, 0.33, 0.83, 1.41, 0.62, 0.28, 0.35, 0.77, 0.17, 0.72, 0.58, 0.45, 0.41), nrow=64,byrow=T)
+  triMutability_Names <- c("AAA", "AAC", "AAG", "AAT", "ACA", "ACC", "ACG", "ACT", "AGA", "AGC", "AGG", "AGT", "ATA", "ATC", "ATG", "ATT", "CAA", "CAC", "CAG", "CAT", "CCA", "CCC", "CCG", "CCT", "CGA", "CGC", "CGG", "CGT", "CTA", "CTC", "CTG", "CTT", "GAA", "GAC", "GAG", "GAT", "GCA", "GCC", "GCG", "GCT", "GGA", "GGC", "GGG", "GGT", "GTA", "GTC", "GTG", "GTT", "TAA", "TAC", "TAG", "TAT", "TCA", "TCC", "TCG", "TCT", "TGA", "TGC", "TGG", "TGT", "TTA", "TTC", "TTG", "TTT")
+  load("FiveS_Mutability.RData")
+
+# Functions
+  
+  # Translate codon to amino acid
+  translateCodonToAminoAcid<-function(Codon){
+     return(AMINO_ACIDS[Codon])
+  }
+
+  # Translate amino acid to trait change
+  translateAminoAcidToTraitChange<-function(AminoAcid){
+     return(TRAITS_AMINO_ACIDS[AminoAcid])
+  }
+    
+  # Initialize Amino Acid Trait Changes
+  initializeTraitChange <- function(traitChangeModel=1,species=1,traitChangeFileName=NULL){
+    if(!is.null(traitChangeFileName)){
+      tryCatch(
+          traitChange <- read.delim(traitChangeFileName,sep="\t",header=T)
+          , error = function(ex){
+            cat("Error|Error reading trait changes. Please check file name/path and format.\n")
+            q()
+          }
+        )
+    }else{
+      traitChange <- TRAITS_AMINO_ACIDS_CHOTHIA98
+    }
+    TRAITS_AMINO_ACIDS <<- traitChange
+ } 
+  
+  # Read in formatted nucleotide substitution matrix
+  initializeSubstitutionMatrix <- function(substitutionModel,species,subsMatFileName=NULL){
+    if(!is.null(subsMatFileName)){
+      tryCatch(
+          subsMat <- read.delim(subsMatFileName,sep="\t",header=T)
+          , error = function(ex){
+            cat("Error|Error reading substitution matrix. Please check file name/path and format.\n")
+            q()
+          }
+        )
+      if(sum(apply(subsMat,1,sum)==1)!=4) subsMat = t(apply(subsMat,1,function(x)x/sum(x)))
+    }else{
+      if(substitutionModel==1)subsMat <- substitution_Literature_Mouse
+      if(substitutionModel==2)subsMat <- substitution_Flu_Human      
+      if(substitutionModel==3)subsMat <- substitution_Flu25_Human      
+       
+    }
+
+    if(substitutionModel==0){
+      subsMat <- matrix(1,4,4)
+      subsMat[,] = 1/3
+      subsMat[1,1] = 0
+      subsMat[2,2] = 0
+      subsMat[3,3] = 0
+      subsMat[4,4] = 0
+    }
+    
+    
+    NUCLEOTIDESN = c(NUCLEOTIDES,"N", "-")
+    if(substitutionModel==5){
+      subsMat <- FiveS_Substitution
+      return(subsMat)
+    }else{
+      subsMat <- rbind(subsMat,rep(NA,4),rep(NA,4))
+      return( matrix(data.matrix(subsMat),6,4,dimnames=list(NUCLEOTIDESN,NUCLEOTIDES) ) )
+    }
+  }
+
+   
+  # Read in formatted Mutability file
+  initializeMutabilityMatrix <- function(mutabilityModel=1, species=1,mutabilityMatFileName=NULL){
+    if(!is.null(mutabilityMatFileName)){
+        tryCatch(
+            mutabilityMat <- read.delim(mutabilityMatFileName,sep="\t",header=T)
+            , error = function(ex){
+              cat("Error|Error reading mutability matrix. Please check file name/path and format.\n")
+              q()
+            }
+          )
+    }else{
+      mutabilityMat <- triMutability_Literature_Human
+      if(species==2) mutabilityMat <- triMutability_Literature_Mouse
+    }
+
+  if(mutabilityModel==0){ mutabilityMat <- matrix(1,64,3)}
+  
+    if(mutabilityModel==5){
+      mutabilityMat <- FiveS_Mutability
+      return(mutabilityMat)
+    }else{
+      return( matrix( data.matrix(mutabilityMat), 64, 3, dimnames=list(triMutability_Names,1:3)) )
+    }
+  }
+
+  # Read FASTA file formats
+  # Modified from read.fasta from the seqinR package
+  baseline.read.fasta <-
+  function (file = system.file("sequences/sample.fasta", package = "seqinr"), 
+      seqtype = c("DNA", "AA"), as.string = FALSE, forceDNAtolower = TRUE, 
+      set.attributes = TRUE, legacy.mode = TRUE, seqonly = FALSE, 
+      strip.desc = FALSE,  sizeof.longlong = .Machine$sizeof.longlong, 
+      endian = .Platform$endian, apply.mask = TRUE) 
+  {
+      seqtype <- match.arg(seqtype)
+  
+          lines <- readLines(file)
+          
+          if (legacy.mode) {
+              comments <- grep("^;", lines)
+              if (length(comments) > 0) 
+                  lines <- lines[-comments]
+          }
+          
+          
+          ind_groups<-which(substr(lines, 1L, 3L) == ">>>")
+          lines_mod<-lines
+  
+          if(!length(ind_groups)){
+              lines_mod<-c(">>>All sequences combined",lines)            
+          }
+          
+          ind_groups<-which(substr(lines_mod, 1L, 3L) == ">>>")
+  
+          lines <- array("BLA",dim=(length(ind_groups)+length(lines_mod)))
+          id<-sapply(1:length(ind_groups),function(i)ind_groups[i]+i-1)+1
+          lines[id] <- "THIS IS A FAKE SEQUENCE"
+          lines[-id] <- lines_mod
+          rm(lines_mod)
+  
+  		ind <- which(substr(lines, 1L, 1L) == ">")
+          nseq <- length(ind)
+          if (nseq == 0) {
+               stop("no line starting with a > character found")
+          }        
+          start <- ind + 1
+          end <- ind - 1
+  
+          while( any(which(ind%in%end)) ){
+            ind=ind[-which(ind%in%end)]
+            nseq <- length(ind)
+            if (nseq == 0) {
+                stop("no line starting with a > character found")
+            }        
+            start <- ind + 1
+            end <- ind - 1        
+          }
+          
+          end <- c(end[-1], length(lines))
+          sequences <- lapply(seq_len(nseq), function(i) paste(lines[start[i]:end[i]], collapse = ""))
+          if (seqonly) 
+              return(sequences)
+          nomseq <- lapply(seq_len(nseq), function(i) {
+          
+              #firstword <- strsplit(lines[ind[i]], " ")[[1]][1]
+              substr(lines[ind[i]], 2, nchar(lines[ind[i]]))
+          
+          })
+          if (seqtype == "DNA") {
+              if (forceDNAtolower) {
+                  sequences <- as.list(tolower(chartr(".","-",sequences)))
+              }else{
+                  sequences <- as.list(toupper(chartr(".","-",sequences)))
+              }
+          }
+          if (as.string == FALSE) 
+              sequences <- lapply(sequences, s2c)
+          if (set.attributes) {
+              for (i in seq_len(nseq)) {
+                  Annot <- lines[ind[i]]
+                  if (strip.desc) 
+                    Annot <- substr(Annot, 2L, nchar(Annot))
+                  attributes(sequences[[i]]) <- list(name = nomseq[[i]], 
+                    Annot = Annot, class = switch(seqtype, AA = "SeqFastaAA", 
+                      DNA = "SeqFastadna"))
+              }
+          }
+          names(sequences) <- nomseq
+          return(sequences)
+  }
+
+  
+  # Replaces non FASTA characters in input files with N  
+  replaceNonFASTAChars <-function(inSeq="ACGTN-AApA"){
+    gsub('[^ACGTNacgt[:punct:]-[:punct:].]','N',inSeq,perl=TRUE)
+  }    
+  
+  # Find the germlines in the FASTA list
+  germlinesInFile <- function(seqIDs){
+    firstChar = sapply(seqIDs,function(x){substr(x,1,1)})
+    secondChar = sapply(seqIDs,function(x){substr(x,2,2)})
+    return(firstChar==">" & secondChar!=">")
+  }
+  
+  # Find the groups in the FASTA list
+  groupsInFile <- function(seqIDs){
+    sapply(seqIDs,function(x){substr(x,1,2)})==">>"
+  }
+
+  # In the process of finding germlines/groups, expand from the start to end of the group
+  expandTillNext <- function(vecPosToID){    
+    IDs = names(vecPosToID)
+    posOfInterests =  which(vecPosToID)
+  
+    expandedID = rep(NA,length(IDs))
+    expandedIDNames = gsub(">","",IDs[posOfInterests])
+    startIndexes = c(1,posOfInterests[-1])
+    stopIndexes = c(posOfInterests[-1]-1,length(IDs))
+    expandedID  = unlist(sapply(1:length(startIndexes),function(i){
+                                    rep(i,stopIndexes[i]-startIndexes[i]+1)
+                                  }))
+    names(expandedID) = unlist(sapply(1:length(startIndexes),function(i){
+                                    rep(expandedIDNames[i],stopIndexes[i]-startIndexes[i]+1)
+                                  }))  
+    return(expandedID)                                                                                                  
+  }
+    
+  # Process FASTA (list) to return a matrix[input, germline)
+  processInputAdvanced <- function(inputFASTA){
+  
+    seqIDs = names(inputFASTA)
+    numbSeqs = length(seqIDs)
+    posGermlines1 = germlinesInFile(seqIDs)
+    numbGermlines = sum(posGermlines1)
+    posGroups1 = groupsInFile(seqIDs)
+    numbGroups = sum(posGroups1)
+    consDef = NA
+    
+    if(numbGermlines==0){
+      posGermlines = 2
+      numbGermlines = 1  
+    }
+  
+      glPositionsSum = cumsum(posGermlines1)
+      glPositions = table(glPositionsSum)
+      #Find the position of the conservation row
+      consDefPos = as.numeric(names(glPositions[names(glPositions)!=0 & glPositions==1]))+1  
+    if( length(consDefPos)> 0 ){
+      consDefID =  match(consDefPos, glPositionsSum) 
+      #The coservation rows need to be pulled out and stores seperately 
+      consDef =  inputFASTA[consDefID]
+      inputFASTA =  inputFASTA[-consDefID]
+  
+      seqIDs = names(inputFASTA)
+      numbSeqs = length(seqIDs)
+      posGermlines1 = germlinesInFile(seqIDs)
+      numbGermlines = sum(posGermlines1)
+      posGroups1 = groupsInFile(seqIDs)
+      numbGroups = sum(posGroups1)
+      if(numbGermlines==0){
+        posGermlines = 2
+        numbGermlines = 1  
+      }    
+    }
+    
+    posGroups <- expandTillNext(posGroups1)
+    posGermlines <- expandTillNext(posGermlines1)
+    posGermlines[posGroups1] = 0
+    names(posGermlines)[posGroups1] = names(posGroups)[posGroups1]
+    posInput = rep(TRUE,numbSeqs)
+    posInput[posGroups1 | posGermlines1] = FALSE
+    
+    matInput = matrix(NA, nrow=sum(posInput), ncol=2)
+    rownames(matInput) = seqIDs[posInput]
+    colnames(matInput) = c("Input","Germline")
+    
+    vecInputFASTA = unlist(inputFASTA)  
+    matInput[,1] = vecInputFASTA[posInput]
+    matInput[,2] = vecInputFASTA[ which( names(inputFASTA)%in%paste(">",names(posGermlines)[posInput],sep="") )[ posGermlines[posInput]] ]
+    
+    germlines = posGermlines[posInput]
+    groups = posGroups[posInput]
+    
+    return( list("matInput"=matInput, "germlines"=germlines, "groups"=groups, "conservationDefinition"=consDef ))      
+  }
+
+
+  # Replace leading and trailing dashes in the sequence
+  replaceLeadingTrailingDashes <- function(x,readEnd){
+    iiGap = unlist(gregexpr("-",x[1]))
+    ggGap = unlist(gregexpr("-",x[2]))  
+    #posToChange = intersect(iiGap,ggGap)
+    
+    
+    seqIn = replaceLeadingTrailingDashesHelper(x[1])
+    seqGL = replaceLeadingTrailingDashesHelper(x[2])
+    seqTemplate = rep('N',readEnd)
+    seqIn <- c(seqIn,seqTemplate[(length(seqIn)+1):readEnd])
+    seqGL <- c(seqGL,seqTemplate[(length(seqGL)+1):readEnd])
+#    if(posToChange!=-1){
+#      seqIn[posToChange] = "-"
+#      seqGL[posToChange] = "-"
+#    }
+  
+    seqIn = c2s(seqIn[1:readEnd])
+    seqGL = c2s(seqGL[1:readEnd])
+  
+    lenGL = nchar(seqGL)
+    if(lenGL<readEnd){
+      seqGL = paste(seqGL,c2s(rep("N",readEnd-lenGL)),sep="")
+    }
+  
+    lenInput = nchar(seqIn)
+    if(lenInput<readEnd){
+      seqIn = paste(seqIn,c2s(rep("N",readEnd-lenInput)),sep="")
+    }    
+    return( c(seqIn,seqGL) )
+  }  
+
+  replaceLeadingTrailingDashesHelper <- function(x){
+    grepResults = gregexpr("-*",x)
+    grepResultsPos = unlist(grepResults)
+    grepResultsLen =  attr(grepResults[[1]],"match.length")   
+    #print(paste("x = '", x, "'", sep=""))
+    x = s2c(x)
+    if(x[1]=="-"){
+      x[1:grepResultsLen[1]] = "N"      
+    }
+    if(x[length(x)]=="-"){
+      x[(length(x)-grepResultsLen[length(grepResultsLen)]+1):length(x)] = "N"      
+    }
+    return(x)
+  }
+
+
+
+  
+  # Check sequences for indels
+  checkForInDels <- function(matInputP){
+    insPos <- checkInsertion(matInputP)
+    delPos <- checkDeletions(matInputP)
+    return(list("Insertions"=insPos, "Deletions"=delPos))
+  }
+
+  # Check sequences for insertions
+  checkInsertion <- function(matInputP){
+    insertionCheck = apply( matInputP,1, function(x){
+                                          inputGaps <- as.vector( gregexpr("-",x[1])[[1]] )
+                                          glGaps <- as.vector( gregexpr("-",x[2])[[1]] )                                          
+                                          return( is.finite( match(FALSE, glGaps%in%inputGaps ) ) )
+                                        })   
+    return(as.vector(insertionCheck))
+  }
+  # Fix inserstions
+  fixInsertions <- function(matInputP){
+    insPos <- checkInsertion(matInputP)
+    sapply((1:nrow(matInputP))[insPos],function(rowIndex){
+                                                x <- matInputP[rowIndex,]
+                                                inputGaps <- gregexpr("-",x[1])[[1]]
+                                                glGaps <- gregexpr("-",x[2])[[1]]
+                                                posInsertions <- glGaps[!(glGaps%in%inputGaps)]
+                                                inputInsertionToN <- s2c(x[2])
+                                                inputInsertionToN[posInsertions]!="-"
+                                                inputInsertionToN[posInsertions] <- "N"
+                                                inputInsertionToN <- c2s(inputInsertionToN)
+                                                matInput[rowIndex,2] <<- inputInsertionToN 
+                                              })                                                               
+    return(insPos)
+  } 
+    
+  # Check sequences for deletions
+  checkDeletions <-function(matInputP){
+    deletionCheck = apply( matInputP,1, function(x){
+                                          inputGaps <- as.vector( gregexpr("-",x[1])[[1]] )
+                                          glGaps <- as.vector( gregexpr("-",x[2])[[1]] )
+                                          return( is.finite( match(FALSE, inputGaps%in%glGaps ) ) )
+                                      })
+    return(as.vector(deletionCheck))                                      
+  }
+  # Fix sequences with deletions
+  fixDeletions <- function(matInputP){
+    delPos <- checkDeletions(matInputP)    
+    sapply((1:nrow(matInputP))[delPos],function(rowIndex){
+                                                x <- matInputP[rowIndex,]
+                                                inputGaps <- gregexpr("-",x[1])[[1]]
+                                                glGaps <- gregexpr("-",x[2])[[1]]
+                                                posDeletions <- inputGaps[!(inputGaps%in%glGaps)]
+                                                inputDeletionToN <- s2c(x[1])
+                                                inputDeletionToN[posDeletions] <- "N"
+                                                inputDeletionToN <- c2s(inputDeletionToN)
+                                                matInput[rowIndex,1] <<- inputDeletionToN 
+                                              })                                                                   
+    return(delPos)
+  }  
+    
+
+  # Trim DNA sequence to the last codon
+  trimToLastCodon <- function(seqToTrim){
+    seqLen = nchar(seqToTrim)  
+    trimmedSeq = s2c(seqToTrim)
+    poi = seqLen
+    tailLen = 0
+    
+    while(trimmedSeq[poi]=="-" || trimmedSeq[poi]=="."){
+      tailLen = tailLen + 1
+      poi = poi - 1   
+    }
+    
+    trimmedSeq = c2s(trimmedSeq[1:(seqLen-tailLen)])
+    seqLen = nchar(trimmedSeq)
+    # Trim sequence to last codon
+  	if( getCodonPos(seqLen)[3] > seqLen )
+  	  trimmedSeq = substr(seqToTrim,1, ( (getCodonPos(seqLen)[1])-1 ) )
+    
+    return(trimmedSeq)
+  }
+  
+  # Given a nuclotide position, returns the pos of the 3 nucs that made the codon
+  # e.g. nuc 86 is part of nucs 85,86,87
+  getCodonPos <- function(nucPos){
+    codonNum =  (ceiling(nucPos/3))*3
+    return( (codonNum-2):codonNum)
+  }
+  
+  # Given a nuclotide position, returns the codon number
+  # e.g. nuc 86  = codon 29
+  getCodonNumb <- function(nucPos){
+    return( ceiling(nucPos/3) )
+  }
+  
+  # Given a codon, returns all the nuc positions that make the codon
+  getCodonNucs <- function(codonNumb){
+    getCodonPos(codonNumb*3)
+  }  
+
+  computeCodonTable <- function(testID=1){
+                  
+    if(testID<=4){    
+      # Pre-compute every codons
+      intCounter = 1
+      for(pOne in NUCLEOTIDES){
+        for(pTwo in NUCLEOTIDES){
+          for(pThree in NUCLEOTIDES){
+            codon = paste(pOne,pTwo,pThree,sep="")
+            colnames(CODON_TABLE)[intCounter] =  codon
+            intCounter = intCounter + 1
+            CODON_TABLE[,codon] = mutationTypeOptimized(cbind(permutateAllCodon(codon),rep(codon,12)))
+          }  
+        }
+      }
+      chars = c("N","A","C","G","T", "-")
+      for(a in chars){
+        for(b in chars){
+          for(c in chars){
+            if(a=="N" | b=="N" | c=="N"){ 
+              #cat(paste(a,b,c),sep="","\n") 
+              CODON_TABLE[,paste(a,b,c,sep="")] = rep(NA,12)
+            }
+          }  
+        }
+      }
+      
+      chars = c("-","A","C","G","T")
+      for(a in chars){
+        for(b in chars){
+          for(c in chars){
+            if(a=="-" | b=="-" | c=="-"){ 
+              #cat(paste(a,b,c),sep="","\n") 
+              CODON_TABLE[,paste(a,b,c,sep="")] = rep(NA,12)
+            }
+          }  
+        }
+      }
+      CODON_TABLE <<- as.matrix(CODON_TABLE)
+    }
+  }
+  
+  collapseClone <- function(vecInputSeqs,glSeq,readEnd,nonTerminalOnly=0){
+  #print(length(vecInputSeqs))
+    vecInputSeqs = unique(vecInputSeqs) 
+    if(length(vecInputSeqs)==1){
+      return( list( c(vecInputSeqs,glSeq), F) )
+    }else{
+      charInputSeqs <- sapply(vecInputSeqs, function(x){
+                                              s2c(x)[1:readEnd]
+                                            })
+      charGLSeq <- s2c(glSeq)
+      matClone <- sapply(1:readEnd, function(i){
+                                            posNucs = unique(charInputSeqs[i,])
+                                            posGL = charGLSeq[i]
+                                            error = FALSE                                            
+                                            if(posGL=="-" & sum(!(posNucs%in%c("-","N")))==0 ){
+                                              return(c("-",error))
+                                            }
+                                            if(length(posNucs)==1)
+                                              return(c(posNucs[1],error))
+                                            else{
+                                              if("N"%in%posNucs){
+                                                error=TRUE
+                                              }
+                                              if(sum(!posNucs[posNucs!="N"]%in%posGL)==0){
+                                                return( c(posGL,error) )  
+                                              }else{
+                                                #return( c(sample(posNucs[posNucs!="N"],1),error) )  
+                                                if(nonTerminalOnly==0){
+                                                  return( c(sample(charInputSeqs[i,charInputSeqs[i,]!="N" & charInputSeqs[i,]!=posGL],1),error) )  
+                                                }else{
+                                                  posNucs = charInputSeqs[i,charInputSeqs[i,]!="N" & charInputSeqs[i,]!=posGL]
+                                                  posNucsTable = table(posNucs)
+                                                  if(sum(posNucsTable>1)==0){
+                                                    return( c(posGL,error) )
+                                                  }else{
+                                                    return( c(sample( posNucs[posNucs%in%names(posNucsTable)[posNucsTable>1]],1),error) )
+                                                  }
+                                                }
+                                                
+                                              }
+                                            } 
+                                          })
+      
+                                          
+      #print(length(vecInputSeqs))                                        
+      return(list(c(c2s(matClone[1,]),glSeq),"TRUE"%in%matClone[2,]))
+    }
+  }
+
+  # Compute the expected for each sequence-germline pair
+  getExpectedIndividual <- function(matInput){
+  if( any(grep("multicore",search())) ){ 
+    facGL <- factor(matInput[,2])
+    facLevels = levels(facGL)
+    LisGLs_MutabilityU = mclapply(1:length(facLevels),  function(x){
+                                                      computeMutabilities(facLevels[x])
+                                                    })
+    facIndex = match(facGL,facLevels)
+    
+    LisGLs_Mutability = mclapply(1:nrow(matInput),  function(x){
+                                                      cInput = rep(NA,nchar(matInput[x,1]))
+                                                      cInput[s2c(matInput[x,1])!="N"] = 1
+                                                      LisGLs_MutabilityU[[facIndex[x]]] * cInput                                                   
+                                                    })
+                                                    
+    LisGLs_Targeting =  mclapply(1:dim(matInput)[1],  function(x){
+                                                      computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
+                                                    })
+                                                    
+    LisGLs_MutationTypes  = mclapply(1:length(matInput[,2]),function(x){
+                                                    #print(x)
+                                                    computeMutationTypes(matInput[x,2])
+                                                })
+    
+    LisGLs_Exp = mclapply(1:dim(matInput)[1],  function(x){
+                                                  computeExpected(LisGLs_Targeting[[x]],LisGLs_MutationTypes[[x]])
+                                                })
+    
+    ul_LisGLs_Exp =  unlist(LisGLs_Exp)                                            
+    return(matrix(ul_LisGLs_Exp,ncol=4,nrow=(length(ul_LisGLs_Exp)/4),byrow=T))
+  }else{
+    facGL <- factor(matInput[,2])
+    facLevels = levels(facGL)
+    LisGLs_MutabilityU = lapply(1:length(facLevels),  function(x){
+      computeMutabilities(facLevels[x])
+    })
+    facIndex = match(facGL,facLevels)
+    
+    LisGLs_Mutability = lapply(1:nrow(matInput),  function(x){
+      cInput = rep(NA,nchar(matInput[x,1]))
+      cInput[s2c(matInput[x,1])!="N"] = 1
+      LisGLs_MutabilityU[[facIndex[x]]] * cInput                                                   
+    })
+    
+    LisGLs_Targeting =  lapply(1:dim(matInput)[1],  function(x){
+      computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
+    })
+    
+    LisGLs_MutationTypes  = lapply(1:length(matInput[,2]),function(x){
+      #print(x)
+      computeMutationTypes(matInput[x,2])
+    })
+    
+    LisGLs_Exp = lapply(1:dim(matInput)[1],  function(x){
+      computeExpected(LisGLs_Targeting[[x]],LisGLs_MutationTypes[[x]])
+    })
+    
+    ul_LisGLs_Exp =  unlist(LisGLs_Exp)                                            
+    return(matrix(ul_LisGLs_Exp,ncol=4,nrow=(length(ul_LisGLs_Exp)/4),byrow=T))
+    
+  }
+  }
+
+  # Compute mutabilities of sequence based on the tri-nucleotide model
+  computeMutabilities <- function(paramSeq){
+    seqLen = nchar(paramSeq)
+    seqMutabilites = rep(NA,seqLen)
+  
+    gaplessSeq = gsub("-", "", paramSeq)
+    gaplessSeqLen = nchar(gaplessSeq)
+    gaplessSeqMutabilites = rep(NA,gaplessSeqLen)
+    
+    if(mutabilityModel!=5){
+      pos<- 3:(gaplessSeqLen)
+      subSeq =  substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))    
+      gaplessSeqMutabilites[pos] =      
+        tapply( c(
+                                        getMutability( substr(subSeq,1,3), 3) , 
+                                        getMutability( substr(subSeq,2,4), 2), 
+                                        getMutability( substr(subSeq,3,5), 1) 
+                                        ),rep(1:(gaplessSeqLen-2),3),mean,na.rm=TRUE
+                                      )
+      #Pos 1
+      subSeq =  substr(gaplessSeq,1,3)
+      gaplessSeqMutabilites[1] =  getMutability(subSeq , 1)
+      #Pos 2
+      subSeq =  substr(gaplessSeq,1,4)
+      gaplessSeqMutabilites[2] =  mean( c(
+                                            getMutability( substr(subSeq,1,3), 2) , 
+                                            getMutability( substr(subSeq,2,4), 1) 
+                                          ),na.rm=T
+                                      ) 
+      seqMutabilites[which(s2c(paramSeq)!="-")]<- gaplessSeqMutabilites
+      return(seqMutabilites)
+    }else{
+      
+      pos<- 3:(gaplessSeqLen)
+      subSeq =  substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))    
+      gaplessSeqMutabilites[pos] = sapply(subSeq,function(x){ getMutability5(x) }, simplify=T)
+      seqMutabilites[which(s2c(paramSeq)!="-")]<- gaplessSeqMutabilites
+      return(seqMutabilites)
+    }
+
+  }
+
+  # Returns the mutability of a triplet at a given position
+  getMutability <- function(codon, pos=1:3){
+    triplets <- rownames(mutability)
+    mutability[  match(codon,triplets) ,pos]
+  }
+
+  getMutability5 <- function(fivemer){
+    return(mutability[fivemer])
+  }
+
+  # Returns the substitution probabilty
+  getTransistionProb <- function(nuc){
+    substitution[nuc,]
+  }
+
+  getTransistionProb5 <- function(fivemer){    
+    if(any(which(fivemer==colnames(substitution)))){
+      return(substitution[,fivemer])
+    }else{
+      return(array(NA,4))
+    }
+  }
+
+  # Given a nuc, returns the other 3 nucs it can mutate to
+  canMutateTo <- function(nuc){
+    NUCLEOTIDES[- which(NUCLEOTIDES==nuc)]
+  }
+  
+  # Given a nucleotide, returns the probabilty of other nucleotide it can mutate to 
+  canMutateToProb <- function(nuc){
+    substitution[nuc,canMutateTo(nuc)]
+  }
+
+  # Compute targeting, based on precomputed mutatbility & substitution  
+  computeTargeting <- function(param_strSeq,param_vecMutabilities){
+
+    if(substitutionModel!=5){
+      vecSeq = s2c(param_strSeq)
+      matTargeting = sapply( 1:length(vecSeq), function(x) { param_vecMutabilities[x] * getTransistionProb(vecSeq[x]) } )  
+      #matTargeting = apply( rbind(vecSeq,param_vecMutabilities),2, function(x) { as.vector(as.numeric(x[2]) * getTransistionProb(x[1])) } )
+      dimnames( matTargeting ) =  list(NUCLEOTIDES,1:(length(vecSeq))) 
+      return (matTargeting)
+    }else{
+      
+      seqLen = nchar(param_strSeq)
+      seqsubstitution = matrix(NA,ncol=seqLen,nrow=4)
+      paramSeq <- param_strSeq
+      gaplessSeq = gsub("-", "", paramSeq)
+      gaplessSeqLen = nchar(gaplessSeq)
+      gaplessSeqSubstitution  = matrix(NA,ncol=gaplessSeqLen,nrow=4) 
+      
+      pos<- 3:(gaplessSeqLen)
+      subSeq =  substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))    
+      gaplessSeqSubstitution[,pos] = sapply(subSeq,function(x){ getTransistionProb5(x) }, simplify=T)
+      seqsubstitution[,which(s2c(paramSeq)!="-")]<- gaplessSeqSubstitution
+      #matTargeting <- param_vecMutabilities  %*% seqsubstitution
+      matTargeting <- sweep(seqsubstitution,2,param_vecMutabilities,`*`)
+      dimnames( matTargeting ) =  list(NUCLEOTIDES,1:(seqLen)) 
+      return (matTargeting)      
+    }
+  }  
+
+  # Compute the mutations types   
+  computeMutationTypes <- function(param_strSeq){
+  #cat(param_strSeq,"\n")
+    #vecSeq = trimToLastCodon(param_strSeq)
+    lenSeq = nchar(param_strSeq)
+    vecCodons = sapply({1:(lenSeq/3)}*3-2,function(x){substr(param_strSeq,x,x+2)})
+    matMutationTypes = matrix( unlist(CODON_TABLE[,vecCodons]) ,ncol=lenSeq,nrow=4, byrow=F)
+    dimnames( matMutationTypes ) =  list(NUCLEOTIDES,1:(ncol(matMutationTypes)))
+    return(matMutationTypes)   
+  }  
+  computeMutationTypesFast <- function(param_strSeq){
+    matMutationTypes = matrix( CODON_TABLE[,param_strSeq] ,ncol=3,nrow=4, byrow=F)
+    #dimnames( matMutationTypes ) =  list(NUCLEOTIDES,1:(length(vecSeq)))
+    return(matMutationTypes)   
+  }  
+  mutationTypeOptimized <- function( matOfCodons ){
+   apply( matOfCodons,1,function(x){ mutationType(x[2],x[1]) } ) 
+  }  
+
+  # Returns a vector of codons 1 mutation away from the given codon
+  permutateAllCodon <- function(codon){
+    cCodon = s2c(codon)
+    matCodons = t(array(cCodon,dim=c(3,12)))
+    matCodons[1:4,1] = NUCLEOTIDES
+    matCodons[5:8,2] = NUCLEOTIDES
+    matCodons[9:12,3] = NUCLEOTIDES
+    apply(matCodons,1,c2s)
+  }
+
+  # Given two codons, tells you if the mutation is R or S (based on your definition)
+  mutationType <- function(codonFrom,codonTo){
+    if(testID==4){
+      if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
+        return(NA)
+      }else{
+        mutationType = "S"
+        if( translateAminoAcidToTraitChange(translateCodonToAminoAcid(codonFrom)) != translateAminoAcidToTraitChange(translateCodonToAminoAcid(codonTo)) ){
+          mutationType = "R"                                                              
+        }
+        if(translateCodonToAminoAcid(codonTo)=="*" | translateCodonToAminoAcid(codonFrom)=="*"){
+          mutationType = "Stop"
+        }
+        return(mutationType)
+      }  
+    }else if(testID==5){  
+      if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
+        return(NA)
+      }else{
+        if(codonFrom==codonTo){
+          mutationType = "S"
+        }else{
+          codonFrom = s2c(codonFrom)
+          codonTo = s2c(codonTo)  
+          mutationType = "Stop"
+          nucOfI = codonFrom[which(codonTo!=codonFrom)]
+          if(nucOfI=="C"){
+            mutationType = "R"  
+          }else if(nucOfI=="G"){
+            mutationType = "S"
+          }
+        }
+        return(mutationType)
+      }
+    }else{
+      if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
+        return(NA)
+      }else{
+        mutationType = "S"
+        if( translateCodonToAminoAcid(codonFrom) != translateCodonToAminoAcid(codonTo) ){
+          mutationType = "R"                                                              
+        }
+        if(translateCodonToAminoAcid(codonTo)=="*" | translateCodonToAminoAcid(codonFrom)=="*"){
+          mutationType = "Stop"
+        }
+        return(mutationType)
+      }  
+    }    
+  }
+
+  
+  #given a mat of targeting & it's corresponding mutationtypes returns 
+  #a vector of Exp_RCDR,Exp_SCDR,Exp_RFWR,Exp_RFWR
+  computeExpected <- function(paramTargeting,paramMutationTypes){
+    # Replacements
+    RPos = which(paramMutationTypes=="R")  
+      #FWR
+      Exp_R_FWR = sum(paramTargeting[ RPos[which(FWR_Nuc_Mat[RPos]==T)] ],na.rm=T)
+      #CDR
+      Exp_R_CDR = sum(paramTargeting[ RPos[which(CDR_Nuc_Mat[RPos]==T)] ],na.rm=T)
+    # Silents
+    SPos = which(paramMutationTypes=="S")  
+      #FWR
+      Exp_S_FWR = sum(paramTargeting[ SPos[which(FWR_Nuc_Mat[SPos]==T)] ],na.rm=T)
+      #CDR
+      Exp_S_CDR = sum(paramTargeting[ SPos[which(CDR_Nuc_Mat[SPos]==T)] ],na.rm=T)
+  
+      return(c(Exp_R_CDR,Exp_S_CDR,Exp_R_FWR,Exp_S_FWR))
+  }
+  
+  # Count the mutations in a sequence
+  # each mutation is treated independently 
+  analyzeMutations2NucUri_website <- function( rev_in_matrix ){
+    paramGL = rev_in_matrix[2,]
+    paramSeq = rev_in_matrix[1,]  
+    
+    #Fill seq with GL seq if gapped
+    #if( any(paramSeq=="-") ){
+    #  gapPos_Seq =  which(paramSeq=="-")
+    #  gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "-"]
+    #  paramSeq[gapPos_Seq_ToReplace] =  paramGL[gapPos_Seq_ToReplace]
+    #}
+  
+  
+    #if( any(paramSeq=="N") ){
+    #  gapPos_Seq =  which(paramSeq=="N")
+    #  gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
+    #  paramSeq[gapPos_Seq_ToReplace] =  paramGL[gapPos_Seq_ToReplace]
+    #}  
+      
+    analyzeMutations2NucUri(  matrix(c( paramGL, paramSeq  ),2,length(paramGL),byrow=T)  )
+    
+  }
+
+  #1 = GL 
+  #2 = Seq
+  analyzeMutations2NucUri <- function( in_matrix=matrix(c(c("A","A","A","C","C","C"),c("A","G","G","C","C","A")),2,6,byrow=T) ){
+    paramGL = in_matrix[2,]
+    paramSeq = in_matrix[1,]
+    paramSeqUri = paramGL
+    #mutations = apply(rbind(paramGL,paramSeq), 2, function(x){!x[1]==x[2]})
+    mutations_val = paramGL != paramSeq   
+    if(any(mutations_val)){
+      mutationPos = {1:length(mutations_val)}[mutations_val]  
+      mutationPos = mutationPos[sapply(mutationPos, function(x){!any(paramSeq[getCodonPos(x)]=="N")})]
+      length_mutations =length(mutationPos)
+      mutationInfo = rep(NA,length_mutations)
+      if(any(mutationPos)){  
+
+        pos<- mutationPos
+        pos_array<-array(sapply(pos,getCodonPos))
+        codonGL =  paramGL[pos_array]
+        
+        codonSeq = sapply(pos,function(x){
+                                  seqP = paramGL[getCodonPos(x)]
+                                  muCodonPos = {x-1}%%3+1 
+                                  seqP[muCodonPos] = paramSeq[x]
+                                  return(seqP)
+                                })      
+        GLcodons =  apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
+        Seqcodons =   apply(codonSeq,2,c2s)
+        mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})     
+        names(mutationInfo) = mutationPos
+    }
+    if(any(!is.na(mutationInfo))){
+      return(mutationInfo[!is.na(mutationInfo)])    
+    }else{
+      return(NA)
+    }
+    
+    
+    }else{
+      return (NA)
+    }
+  }
+  
+  processNucMutations2 <- function(mu){
+    if(!is.na(mu)){
+      #R
+      if(any(mu=="R")){
+        Rs = mu[mu=="R"]
+        nucNumbs = as.numeric(names(Rs))
+        R_CDR = sum(as.integer(CDR_Nuc[nucNumbs]),na.rm=T)
+        R_FWR = sum(as.integer(FWR_Nuc[nucNumbs]),na.rm=T)      
+      }else{
+        R_CDR = 0
+        R_FWR = 0
+      }    
+      
+      #S
+      if(any(mu=="S")){
+        Ss = mu[mu=="S"]
+        nucNumbs = as.numeric(names(Ss))
+        S_CDR = sum(as.integer(CDR_Nuc[nucNumbs]),na.rm=T)
+        S_FWR = sum(as.integer(FWR_Nuc[nucNumbs]),na.rm=T)      
+      }else{
+        S_CDR = 0
+        S_FWR = 0
+      }    
+      
+      
+      retVec = c(R_CDR,S_CDR,R_FWR,S_FWR)
+      retVec[is.na(retVec)]=0
+      return(retVec)
+    }else{
+      return(rep(0,4))
+    }
+  }        
+  
+  
+  ## Z-score Test
+  computeZScore <- function(mat, test="Focused"){
+    matRes <- matrix(NA,ncol=2,nrow=(nrow(mat)))
+    if(test=="Focused"){
+      #Z_Focused_CDR
+      #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
+      P = apply(mat[,c(5,6,8)],1,function(x){(x[1]/sum(x))})
+      R_mean = apply(cbind(mat[,c(1,2,4)],P),1,function(x){x[4]*(sum(x[1:3]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,1] = (mat[,1]-R_mean)/R_sd
+    
+      #Z_Focused_FWR
+      #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
+      P = apply(mat[,c(7,6,8)],1,function(x){(x[1]/sum(x))})
+      R_mean = apply(cbind(mat[,c(3,2,4)],P),1,function(x){x[4]*(sum(x[1:3]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,2] = (mat[,3]-R_mean)/R_sd
+    }
+  
+    if(test=="Local"){
+      #Z_Focused_CDR
+      #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
+      P = apply(mat[,c(5,6)],1,function(x){(x[1]/sum(x))})
+      R_mean = apply(cbind(mat[,c(1,2)],P),1,function(x){x[3]*(sum(x[1:2]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,1] = (mat[,1]-R_mean)/R_sd
+    
+      #Z_Focused_FWR
+      #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
+      P = apply(mat[,c(7,8)],1,function(x){(x[1]/sum(x))})
+      R_mean = apply(cbind(mat[,c(3,4)],P),1,function(x){x[3]*(sum(x[1:2]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,2] = (mat[,3]-R_mean)/R_sd
+    }
+    
+    if(test=="Imbalanced"){
+      #Z_Focused_CDR
+      #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
+      P = apply(mat[,5:8],1,function(x){((x[1]+x[2])/sum(x))})
+      R_mean = apply(cbind(mat[,1:4],P),1,function(x){x[5]*(sum(x[1:4]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,1] = (mat[,1]-R_mean)/R_sd
+    
+      #Z_Focused_FWR
+      #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
+      P = apply(mat[,5:8],1,function(x){((x[3]+x[4])/sum(x))})
+      R_mean = apply(cbind(mat[,1:4],P),1,function(x){x[5]*(sum(x[1:4]))})
+      R_sd=sqrt(R_mean*(1-P))
+      matRes[,2] = (mat[,3]-R_mean)/R_sd
+    }    
+      
+    matRes[is.nan(matRes)] = NA
+    return(matRes)
+  }
+
+  # Return a p-value for a z-score
+  z2p <- function(z){
+    p=NA
+    if( !is.nan(z) && !is.na(z)){   
+      if(z>0){
+        p = (1 - pnorm(z,0,1))
+      } else if(z<0){
+        p = (-1 * pnorm(z,0,1))
+      } else{
+        p = 0.5
+      }
+    }else{
+      p = NA
+    }
+    return(p)
+  }    
+  
+  
+  ## Bayesian  Test
+
+  # Fitted parameter for the bayesian framework
+BAYESIAN_FITTED<-c(0.407277142798302, 0.554007336744485, 0.63777155771234, 0.693989162719009, 0.735450014674917, 0.767972534429806, 0.794557287143399, 0.816906816601605, 0.83606796225341, 0.852729446430296, 0.867370424541641, 0.880339760590323, 0.891900995024999, 0.902259181289864, 0.911577919359,0.919990301665853, 0.927606458124537, 0.934518806350661, 0.940805863754375, 0.946534836475715, 0.951763691199255, 0.95654428191308, 0.960920179487397, 0.964930893680829, 0.968611312149038, 0.971992459313836, 0.975102110004818, 0.977964943023096, 0.980603428208439, 0.983037660179428, 0.985285800977406, 0.987364285326685, 0.989288037855441, 0.991070478823525, 0.992723699729969, 0.994259575477392, 0.995687688867975, 0.997017365051493, 0.998257085153047, 0.999414558305388, 1.00049681357804, 1.00151036237481, 1.00246080204981, 1.00335370751909, 1.0041939329768, 1.0049859393417, 1.00573382091263, 1.00644127217376, 1.00711179729107, 1.00774845526417, 1.00835412715854, 1.00893143010366, 1.00948275846309, 1.01001030293661, 1.01051606798079, 1.01100188771288, 1.01146944044216, 1.01192026195449, 1.01235575766094, 1.01277721370986)
+  CONST_i <- sort(c(((2^(seq(-39,0,length.out=201)))/2)[1:200],(c(0:11,13:99)+0.5)/100,1-(2^(seq(-39,0,length.out=201)))/2))
+  
+  # Given x, M & p, returns a pdf 
+  calculate_bayes <- function ( x=3, N=10, p=0.33,
+                                i=CONST_i,
+                                max_sigma=20,length_sigma=4001
+                              ){
+    if(!0%in%N){
+      G <- max(length(x),length(N),length(p))
+      x=array(x,dim=G)
+      N=array(N,dim=G)
+      p=array(p,dim=G)
+      sigma_s<-seq(-max_sigma,max_sigma,length.out=length_sigma)
+      sigma_1<-log({i/{1-i}}/{p/{1-p}})
+      index<-min(N,60)
+      y<-dbeta(i,x+BAYESIAN_FITTED[index],N+BAYESIAN_FITTED[index]-x)*(1-p)*p*exp(sigma_1)/({1-p}^2+2*p*{1-p}*exp(sigma_1)+{p^2}*exp(2*sigma_1))
+      if(!sum(is.na(y))){
+        tmp<-approx(sigma_1,y,sigma_s)$y
+        tmp/sum(tmp)/{2*max_sigma/{length_sigma-1}}
+      }else{
+        return(NA)
+      }
+    }else{
+      return(NA)
+    }
+  }  
+  # Given a mat of observed & expected, return a list of CDR & FWR pdf for selection
+  computeBayesianScore <- function(mat, test="Focused", max_sigma=20,length_sigma=4001){
+    flagOneSeq = F
+    if(nrow(mat)==1){
+      mat=rbind(mat,mat)
+      flagOneSeq = T
+    }
+    if(test=="Focused"){
+      #CDR
+      P = c(apply(mat[,c(5,6,8)],1,function(x){(x[1]/sum(x))}),0.5)
+      N = c(apply(mat[,c(1,2,4)],1,function(x){(sum(x))}),0)
+      X = c(mat[,1],0)
+      bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesCDR = bayesCDR[-length(bayesCDR)]
+  
+      #FWR
+      P = c(apply(mat[,c(7,6,8)],1,function(x){(x[1]/sum(x))}),0.5)
+      N = c(apply(mat[,c(3,2,4)],1,function(x){(sum(x))}),0)
+      X = c(mat[,3],0)
+      bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesFWR = bayesFWR[-length(bayesFWR)]     
+    }
+    
+    if(test=="Local"){
+      #CDR
+      P = c(apply(mat[,c(5,6)],1,function(x){(x[1]/sum(x))}),0.5)
+      N = c(apply(mat[,c(1,2)],1,function(x){(sum(x))}),0)
+      X = c(mat[,1],0)
+      bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesCDR = bayesCDR[-length(bayesCDR)]
+  
+      #FWR
+      P = c(apply(mat[,c(7,8)],1,function(x){(x[1]/sum(x))}),0.5)
+      N = c(apply(mat[,c(3,4)],1,function(x){(sum(x))}),0)
+      X = c(mat[,3],0)
+      bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesFWR = bayesFWR[-length(bayesFWR)]     
+    } 
+     
+    if(test=="Imbalanced"){
+      #CDR
+      P = c(apply(mat[,c(5:8)],1,function(x){((x[1]+x[2])/sum(x))}),0.5)
+      N = c(apply(mat[,c(1:4)],1,function(x){(sum(x))}),0)
+      X = c(apply(mat[,c(1:2)],1,function(x){(sum(x))}),0)
+      bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesCDR = bayesCDR[-length(bayesCDR)]
+  
+      #FWR
+      P = c(apply(mat[,c(5:8)],1,function(x){((x[3]+x[4])/sum(x))}),0.5)
+      N = c(apply(mat[,c(1:4)],1,function(x){(sum(x))}),0)
+      X = c(apply(mat[,c(3:4)],1,function(x){(sum(x))}),0)
+      bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesFWR = bayesFWR[-length(bayesFWR)]     
+    }
+
+    if(test=="ImbalancedSilent"){
+      #CDR
+      P = c(apply(mat[,c(6,8)],1,function(x){((x[1])/sum(x))}),0.5)
+      N = c(apply(mat[,c(2,4)],1,function(x){(sum(x))}),0)
+      X = c(apply(mat[,c(2,4)],1,function(x){(x[1])}),0)
+      bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesCDR = bayesCDR[-length(bayesCDR)]
+  
+      #FWR
+      P = c(apply(mat[,c(6,8)],1,function(x){((x[2])/sum(x))}),0.5)
+      N = c(apply(mat[,c(2,4)],1,function(x){(sum(x))}),0)
+      X = c(apply(mat[,c(2,4)],1,function(x){(x[2])}),0)
+      bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})    
+      bayesFWR = bayesFWR[-length(bayesFWR)]     
+    }
+        
+    if(flagOneSeq==T){
+      bayesCDR = bayesCDR[1]  
+      bayesFWR = bayesFWR[1]
+    }
+    return( list("CDR"=bayesCDR, "FWR"=bayesFWR) )
+  }
+  
+  ##Covolution
+  break2chunks<-function(G=1000){
+  base<-2^round(log(sqrt(G),2),0)
+  return(c(rep(base,floor(G/base)-1),base+G-(floor(G/base)*base)))
+  }  
+  
+  PowersOfTwo <- function(G=100){
+    exponents <- array()
+    i = 0
+    while(G > 0){
+      i=i+1
+      exponents[i] <- floor( log2(G) )
+      G <- G-2^exponents[i]
+    }
+    return(exponents)
+  }
+  
+  convolutionPowersOfTwo <- function( cons, length_sigma=4001 ){
+    G = ncol(cons)
+    if(G>1){
+      for(gen in log(G,2):1){
+        ll<-seq(from=2,to=2^gen,by=2)
+        sapply(ll,function(l){cons[,l/2]<<-weighted_conv(cons[,l],cons[,l-1],length_sigma=length_sigma)})
+      }
+    }
+    return( cons[,1] )
+  }
+  
+  convolutionPowersOfTwoByTwos <- function( cons, length_sigma=4001,G=1 ){
+    if(length(ncol(cons))) G<-ncol(cons)
+    groups <- PowersOfTwo(G)
+    matG <- matrix(NA, ncol=length(groups), nrow=length(cons)/G )
+    startIndex = 1
+    for( i in 1:length(groups) ){
+      stopIndex <- 2^groups[i] + startIndex - 1
+      if(stopIndex!=startIndex){
+        matG[,i] <- convolutionPowersOfTwo( cons[,startIndex:stopIndex], length_sigma=length_sigma )
+        startIndex = stopIndex + 1
+      }
+      else {
+        if(G>1) matG[,i] <- cons[,startIndex:stopIndex]
+        else matG[,i] <- cons
+        #startIndex = stopIndex + 1
+      }
+    }
+    return( list( matG, groups ) )
+  }
+  
+  weighted_conv<-function(x,y,w=1,m=100,length_sigma=4001){
+    lx<-length(x)
+    ly<-length(y)
+    if({lx<m}| {{lx*w}<m}| {{ly}<m}| {{ly*w}<m}){
+      if(w<1){
+        y1<-approx(1:ly,y,seq(1,ly,length.out=m))$y
+        x1<-approx(1:lx,x,seq(1,lx,length.out=m/w))$y
+        lx<-length(x1)
+        ly<-length(y1)
+      }
+      else {
+        y1<-approx(1:ly,y,seq(1,ly,length.out=m*w))$y
+        x1<-approx(1:lx,x,seq(1,lx,length.out=m))$y
+        lx<-length(x1)
+        ly<-length(y1)
+      }
+    }
+    else{
+      x1<-x
+      y1<-approx(1:ly,y,seq(1,ly,length.out=floor(lx*w)))$y
+      ly<-length(y1)
+    }
+    tmp<-approx(x=1:(lx+ly-1),y=convolve(x1,rev(y1),type="open"),xout=seq(1,lx+ly-1,length.out=length_sigma))$y
+    tmp[tmp<=0] = 0
+    return(tmp/sum(tmp))
+  }
+  
+  calculate_bayesGHelper <- function( listMatG,length_sigma=4001 ){
+    matG <- listMatG[[1]]
+    groups <- listMatG[[2]]
+    i = 1
+    resConv <- matG[,i]
+    denom <- 2^groups[i]
+    if(length(groups)>1){
+      while( i<length(groups) ){
+        i = i + 1
+        resConv <- weighted_conv(resConv, matG[,i], w= {{2^groups[i]}/denom} ,length_sigma=length_sigma)
+        #cat({{2^groups[i]}/denom},"\n")
+        denom <- denom + 2^groups[i]
+      }
+    }
+    return(resConv)
+  }
+  
+  # Given a list of PDFs, returns a convoluted PDF    
+  groupPosteriors <- function( listPosteriors, max_sigma=20, length_sigma=4001 ,Threshold=2 ){  
+    listPosteriors = listPosteriors[ !is.na(listPosteriors) ]
+    Length_Postrior<-length(listPosteriors)
+    if(Length_Postrior>1 & Length_Postrior<=Threshold){
+      cons = matrix(unlist(listPosteriors),length(listPosteriors[[1]]),length(listPosteriors))
+      listMatG <- convolutionPowersOfTwoByTwos(cons,length_sigma=length_sigma)
+      y<-calculate_bayesGHelper(listMatG,length_sigma=length_sigma)
+      return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
+    }else if(Length_Postrior==1) return(listPosteriors[[1]])
+    else  if(Length_Postrior==0) return(NA)
+    else {
+      cons = matrix(unlist(listPosteriors),length(listPosteriors[[1]]),length(listPosteriors))
+      y = fastConv(cons,max_sigma=max_sigma, length_sigma=length_sigma )
+      return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
+    }
+  }
+
+  fastConv<-function(cons, max_sigma=20, length_sigma=4001){
+    chunks<-break2chunks(G=ncol(cons))
+    if(ncol(cons)==3) chunks<-2:1
+    index_chunks_end <- cumsum(chunks)
+    index_chunks_start <- c(1,index_chunks_end[-length(index_chunks_end)]+1)
+    index_chunks <- cbind(index_chunks_start,index_chunks_end)
+    
+    case <- sum(chunks!=chunks[1])
+    if(case==1) End <- max(1,((length(index_chunks)/2)-1))
+    else End <- max(1,((length(index_chunks)/2)))
+    
+    firsts <- sapply(1:End,function(i){
+          	    indexes<-index_chunks[i,1]:index_chunks[i,2]
+          	    convolutionPowersOfTwoByTwos(cons[ ,indexes])[[1]]
+          	  })
+    if(case==0){
+    	result<-calculate_bayesGHelper( convolutionPowersOfTwoByTwos(firsts) )
+    }else if(case==1){
+      last<-list(calculate_bayesGHelper(
+      convolutionPowersOfTwoByTwos( cons[ ,index_chunks[length(index_chunks)/2,1]:index_chunks[length(index_chunks)/2,2]] )
+                                      ),0)
+      result_first<-calculate_bayesGHelper(convolutionPowersOfTwoByTwos(firsts))
+      result<-calculate_bayesGHelper(
+        list(
+          cbind(
+          result_first,last[[1]]),
+          c(log(index_chunks_end[length(index_chunks)/2-1],2),log(index_chunks[length(index_chunks)/2,2]-index_chunks[length(index_chunks)/2,1]+1,2))
+        )
+      )
+    }
+    return(as.vector(result))
+  }
+    
+  # Computes the 95% CI for a pdf
+  calcBayesCI <- function(Pdf,low=0.025,up=0.975,max_sigma=20, length_sigma=4001){
+    if(length(Pdf)!=length_sigma) return(NA)
+    sigma_s=seq(-max_sigma,max_sigma,length.out=length_sigma)
+    cdf = cumsum(Pdf)
+    cdf = cdf/cdf[length(cdf)]  
+    return( c(sigma_s[findInterval(low,cdf)-1] , sigma_s[findInterval(up,cdf)]) ) 
+  }
+  
+  # Computes a mean for a pdf
+  calcBayesMean <- function(Pdf,max_sigma=20,length_sigma=4001){
+    if(length(Pdf)!=length_sigma) return(NA)
+    sigma_s=seq(-max_sigma,max_sigma,length.out=length_sigma)
+    norm = {length_sigma-1}/2/max_sigma
+    return( (Pdf%*%sigma_s/norm)  ) 
+  }
+  
+  # Returns the mean, and the 95% CI for a pdf
+  calcBayesOutputInfo <- function(Pdf,low=0.025,up=0.975,max_sigma=20, length_sigma=4001){
+    if(is.na(Pdf)) 
+     return(rep(NA,3))  
+    bCI = calcBayesCI(Pdf=Pdf,low=low,up=up,max_sigma=max_sigma,length_sigma=length_sigma)
+    bMean = calcBayesMean(Pdf=Pdf,max_sigma=max_sigma,length_sigma=length_sigma)
+    return(c(bMean, bCI))
+  }   
+
+  # Computes the p-value of a pdf
+  computeSigmaP <- function(Pdf, length_sigma=4001, max_sigma=20){
+    if(length(Pdf)>1){
+      norm = {length_sigma-1}/2/max_sigma
+      pVal = {sum(Pdf[1:{{length_sigma-1}/2}]) + Pdf[{{length_sigma+1}/2}]/2}/norm
+      if(pVal>0.5){
+        pVal = pVal-1
+      }
+      return(pVal)
+    }else{
+      return(NA)
+    }
+  }    
+  
+  # Compute p-value of two distributions
+  compareTwoDistsFaster <-function(sigma_S=seq(-20,20,length.out=4001), N=10000, dens1=runif(4001,0,1), dens2=runif(4001,0,1)){
+  #print(c(length(dens1),length(dens2)))
+  if(length(dens1)>1 & length(dens2)>1 ){
+    dens1<-dens1/sum(dens1)
+    dens2<-dens2/sum(dens2)
+    cum2 <- cumsum(dens2)-dens2/2
+    tmp<- sum(sapply(1:length(dens1),function(i)return(dens1[i]*cum2[i])))
+    #print(tmp)
+    if(tmp>0.5)tmp<-tmp-1
+    return( tmp )
+    }
+    else {
+    return(NA)
+    }
+    #return (sum(sapply(1:N,function(i)(sample(sigma_S,1,prob=dens1)>sample(sigma_S,1,prob=dens2))))/N)
+  }  
+  
+  # get number of seqeunces contributing to the sigma (i.e. seqeunces with mutations)
+  numberOfSeqsWithMutations <- function(matMutations,test=1){
+    if(test==4)test=2
+    cdrSeqs <- 0
+    fwrSeqs <- 0    
+    if(test==1){#focused
+      cdrMutations <- apply(matMutations, 1, function(x){ sum(x[c(1,2,4)]) })
+      fwrMutations <- apply(matMutations, 1, function(x){ sum(x[c(3,4,2)]) })
+      if( any(which(cdrMutations>0)) ) cdrSeqs <- sum(cdrMutations>0)
+      if( any(which(fwrMutations>0)) ) fwrSeqs <- sum(fwrMutations>0) 
+    }
+    if(test==2){#local
+      cdrMutations <- apply(matMutations, 1, function(x){ sum(x[c(1,2)]) })
+      fwrMutations <- apply(matMutations, 1, function(x){ sum(x[c(3,4)]) })
+      if( any(which(cdrMutations>0)) ) cdrSeqs <- sum(cdrMutations>0)
+      if( any(which(fwrMutations>0)) ) fwrSeqs <- sum(fwrMutations>0) 
+    }
+  return(c("CDR"=cdrSeqs, "FWR"=fwrSeqs))
+}  
+
+
+
+shadeColor <- function(sigmaVal=NA,pVal=NA){
+  if(is.na(sigmaVal) & is.na(pVal)) return(NA)
+  if(is.na(sigmaVal) & !is.na(pVal)) sigmaVal=sign(pVal)
+  if(is.na(pVal) || pVal==1 || pVal==0){
+    returnColor = "#FFFFFF";
+  }else{
+    colVal=abs(pVal);
+    
+    if(sigmaVal<0){      
+        if(colVal>0.1)
+          returnColor = "#CCFFCC";
+        if(colVal<=0.1)
+          returnColor = "#99FF99";
+        if(colVal<=0.050)
+          returnColor = "#66FF66";
+        if(colVal<=0.010)
+          returnColor = "#33FF33";
+        if(colVal<=0.005)
+          returnColor = "#00FF00";
+      
+    }else{
+      if(colVal>0.1)
+        returnColor = "#FFCCCC";
+      if(colVal<=0.1)
+        returnColor = "#FF9999";
+      if(colVal<=0.05)
+        returnColor = "#FF6666";
+      if(colVal<=0.01)
+        returnColor = "#FF3333";
+      if(colVal<0.005)
+        returnColor = "#FF0000";
+    }
+  }
+  
+  return(returnColor)
+}
+
+
+
+plotHelp <- function(xfrac=0.05,yfrac=0.05,log=FALSE){
+  if(!log){
+    x = par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac
+    y = par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac
+  }else {
+    if(log==2){
+      x = par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac
+      y = 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac)
+    }
+    if(log==1){
+      x = 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac)
+      y = par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac
+    }
+    if(log==3){
+      x = 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac)
+      y = 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac)
+    }
+  }
+  return(c("x"=x,"y"=y))
+}
+
+# SHMulation
+
+  # Based on targeting, introduce a single mutation & then update the targeting 
+  oneMutation <- function(){
+    # Pick a postion + mutation
+    posMutation = sample(1:(seqGermlineLen*4),1,replace=F,prob=as.vector(seqTargeting))
+    posNucNumb = ceiling(posMutation/4)                    # Nucleotide number
+    posNucKind = 4 - ( (posNucNumb*4) - posMutation )   # Nuc the position mutates to
+  
+    #mutate the simulation sequence
+    seqSimVec <-  s2c(seqSim)
+    seqSimVec[posNucNumb] <- NUCLEOTIDES[posNucKind]
+    seqSim <<-  c2s(seqSimVec)
+    
+    #update Mutability, Targeting & MutationsTypes
+    updateMutabilityNTargeting(posNucNumb)
+  
+    #return(c(posNucNumb,NUCLEOTIDES[posNucKind])) 
+    return(posNucNumb)
+  }  
+  
+  updateMutabilityNTargeting <- function(position){
+    min_i<-max((position-2),1)
+    max_i<-min((position+2),nchar(seqSim))
+    min_ii<-min(min_i,3)
+    
+    #mutability - update locally
+    seqMutability[(min_i):(max_i)] <<- computeMutabilities(substr(seqSim,position-4,position+4))[(min_ii):(max_i-min_i+min_ii)]
+    
+    
+    #targeting - compute locally
+    seqTargeting[,min_i:max_i] <<- computeTargeting(substr(seqSim,min_i,max_i),seqMutability[min_i:max_i])                 
+    seqTargeting[is.na(seqTargeting)] <<- 0
+    #mutCodonPos = getCodonPos(position) 
+    mutCodonPos = seq(getCodonPos(min_i)[1],getCodonPos(max_i)[3])
+    #cat(mutCodonPos,"\n")                                                  
+    mutTypeCodon = getCodonPos(position)
+    seqMutationTypes[,mutTypeCodon] <<- computeMutationTypesFast( substr(seqSim,mutTypeCodon[1],mutTypeCodon[3]) ) 
+    # Stop = 0
+    if(any(seqMutationTypes[,mutCodonPos]=="Stop",na.rm=T )){
+      seqTargeting[,mutCodonPos][seqMutationTypes[,mutCodonPos]=="Stop"] <<- 0
+    }
+    
+  
+    #Selection
+    selectedPos = (min_i*4-4)+(which(seqMutationTypes[,min_i:max_i]=="R"))  
+    # CDR
+    selectedCDR = selectedPos[which(matCDR[selectedPos]==T)]
+    seqTargeting[selectedCDR] <<-  seqTargeting[selectedCDR] *  exp(selCDR)
+    seqTargeting[selectedCDR] <<- seqTargeting[selectedCDR]/baseLineCDR_K
+        
+    # FWR
+    selectedFWR = selectedPos[which(matFWR[selectedPos]==T)]
+    seqTargeting[selectedFWR] <<-  seqTargeting[selectedFWR] *  exp(selFWR)
+    seqTargeting[selectedFWR] <<- seqTargeting[selectedFWR]/baseLineFWR_K      
+    
+  }  
+  
+
+
+  # Validate the mutation: if the mutation has not been sampled before validate it, else discard it.   
+  validateMutation <- function(){  
+    if( !(mutatedPos%in%mutatedPositions) ){ # if it's a new mutation
+      uniqueMutationsIntroduced <<- uniqueMutationsIntroduced + 1
+      mutatedPositions[uniqueMutationsIntroduced] <<-  mutatedPos  
+    }else{
+      if(substr(seqSim,mutatedPos,mutatedPos)==substr(seqGermline,mutatedPos,mutatedPos)){ # back to germline mutation
+        mutatedPositions <<-  mutatedPositions[-which(mutatedPositions==mutatedPos)]
+        uniqueMutationsIntroduced <<-  uniqueMutationsIntroduced - 1
+      }      
+    }
+  }  
+  
+  
+  
+  # Places text (labels) at normalized coordinates 
+  myaxis <- function(xfrac=0.05,yfrac=0.05,log=FALSE,w="text",cex=1,adj=1,thecol="black"){
+    par(xpd=TRUE)
+    if(!log)
+      text(par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac,par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac,w,cex=cex,adj=adj,col=thecol)
+    else {
+    if(log==2)
+    text(
+      par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac,
+      10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac),
+      w,cex=cex,adj=adj,col=thecol)
+    if(log==1)
+      text(
+      10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac),
+      par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac,
+      w,cex=cex,adj=adj,col=thecol)
+    if(log==3)
+      text(
+      10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac),
+      10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac),
+      w,cex=cex,adj=adj,col=thecol)
+    }
+    par(xpd=FALSE)
+  }
+  
+  
+  
+  # Count the mutations in a sequence
+  analyzeMutations <- function( inputMatrixIndex, model = 0 , multipleMutation=0, seqWithStops=0){
+
+    paramGL = s2c(matInput[inputMatrixIndex,2])
+    paramSeq = s2c(matInput[inputMatrixIndex,1])            
+    
+    #if( any(paramSeq=="N") ){
+    #  gapPos_Seq =  which(paramSeq=="N")
+    #  gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
+    #  paramSeq[gapPos_Seq_ToReplace] =  paramGL[gapPos_Seq_ToReplace]
+    #}        
+    mutations_val = paramGL != paramSeq   
+    
+    if(any(mutations_val)){
+      mutationPos = which(mutations_val)#{1:length(mutations_val)}[mutations_val]  
+      length_mutations =length(mutationPos)
+      mutationInfo = rep(NA,length_mutations)
+                          
+      pos<- mutationPos
+      pos_array<-array(sapply(pos,getCodonPos))
+      codonGL =  paramGL[pos_array]
+      codonSeqWhole =  paramSeq[pos_array]
+      codonSeq = sapply(pos,function(x){
+                                seqP = paramGL[getCodonPos(x)]
+                                muCodonPos = {x-1}%%3+1 
+                                seqP[muCodonPos] = paramSeq[x]
+                                return(seqP)
+                              })
+      GLcodons =  apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
+      SeqcodonsWhole =  apply(matrix(codonSeqWhole,length_mutations,3,byrow=TRUE),1,c2s)      
+      Seqcodons =   apply(codonSeq,2,c2s)
+      
+      mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})     
+      names(mutationInfo) = mutationPos     
+      
+      mutationInfoWhole = apply(rbind(GLcodons , SeqcodonsWhole),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})           
+      names(mutationInfoWhole) = mutationPos
+
+      mutationInfo <- mutationInfo[!is.na(mutationInfo)]
+      mutationInfoWhole <- mutationInfoWhole[!is.na(mutationInfoWhole)]
+      
+      if(any(!is.na(mutationInfo))){       
+  
+        #Filter based on Stop (at the codon level)
+        if(seqWithStops==1){
+          nucleotidesAtStopCodons = names(mutationInfoWhole[mutationInfoWhole!="Stop"])
+          mutationInfo = mutationInfo[nucleotidesAtStopCodons]
+          mutationInfoWhole = mutationInfo[nucleotidesAtStopCodons]
+        }else{
+          countStops = sum(mutationInfoWhole=="Stop")
+          if(seqWithStops==2 & countStops==0) mutationInfo = NA
+          if(seqWithStops==3 & countStops>0) mutationInfo = NA
+        }         
+        
+        if(any(!is.na(mutationInfo))){
+          #Filter mutations based on multipleMutation
+          if(multipleMutation==1 & !is.na(mutationInfo)){
+            mutationCodons = getCodonNumb(as.numeric(names(mutationInfoWhole)))
+            tableMutationCodons <- table(mutationCodons)
+            codonsWithMultipleMutations <- as.numeric(names(tableMutationCodons[tableMutationCodons>1]))
+            if(any(codonsWithMultipleMutations)){
+              #remove the nucleotide mutations in the codons with multiple mutations
+              mutationInfo <- mutationInfo[!(mutationCodons %in% codonsWithMultipleMutations)]
+              #replace those codons with Ns in the input sequence
+              paramSeq[unlist(lapply(codonsWithMultipleMutations, getCodonNucs))] = "N"
+              matInput[inputMatrixIndex,1] <<- c2s(paramSeq)
+            }
+          }
+
+          #Filter mutations based on the model
+          if(any(mutationInfo)==T | is.na(any(mutationInfo))){        
+            
+            if(model==1 & !is.na(mutationInfo)){
+              mutationInfo <- mutationInfo[mutationInfo=="S"]
+            }  
+            if(any(mutationInfo)==T | is.na(any(mutationInfo))) return(mutationInfo)
+            else return(NA)
+          }else{
+            return(NA)
+          }
+        }else{
+          return(NA)
+        }
+        
+        
+      }else{
+        return(NA)
+      }
+    
+    
+    }else{
+      return (NA)
+    }    
+  }  
+
+   analyzeMutationsFixed <- function( inputArray, model = 0 , multipleMutation=0, seqWithStops=0){
+
+    paramGL = s2c(inputArray[2])
+    paramSeq = s2c(inputArray[1])            
+    inputSeq <- inputArray[1]
+    #if( any(paramSeq=="N") ){
+    #  gapPos_Seq =  which(paramSeq=="N")
+    #  gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
+    #  paramSeq[gapPos_Seq_ToReplace] =  paramGL[gapPos_Seq_ToReplace]
+    #}        
+    mutations_val = paramGL != paramSeq   
+    
+    if(any(mutations_val)){
+      mutationPos = which(mutations_val)#{1:length(mutations_val)}[mutations_val]  
+      length_mutations =length(mutationPos)
+      mutationInfo = rep(NA,length_mutations)
+                          
+      pos<- mutationPos
+      pos_array<-array(sapply(pos,getCodonPos))
+      codonGL =  paramGL[pos_array]
+      codonSeqWhole =  paramSeq[pos_array]
+      codonSeq = sapply(pos,function(x){
+                                seqP = paramGL[getCodonPos(x)]
+                                muCodonPos = {x-1}%%3+1 
+                                seqP[muCodonPos] = paramSeq[x]
+                                return(seqP)
+                              })
+      GLcodons =  apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
+      SeqcodonsWhole =  apply(matrix(codonSeqWhole,length_mutations,3,byrow=TRUE),1,c2s)      
+      Seqcodons =   apply(codonSeq,2,c2s)
+      
+      mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})     
+      names(mutationInfo) = mutationPos     
+      
+      mutationInfoWhole = apply(rbind(GLcodons , SeqcodonsWhole),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})           
+      names(mutationInfoWhole) = mutationPos
+
+      mutationInfo <- mutationInfo[!is.na(mutationInfo)]
+      mutationInfoWhole <- mutationInfoWhole[!is.na(mutationInfoWhole)]
+      
+      if(any(!is.na(mutationInfo))){       
+  
+        #Filter based on Stop (at the codon level)
+        if(seqWithStops==1){
+          nucleotidesAtStopCodons = names(mutationInfoWhole[mutationInfoWhole!="Stop"])
+          mutationInfo = mutationInfo[nucleotidesAtStopCodons]
+          mutationInfoWhole = mutationInfo[nucleotidesAtStopCodons]
+        }else{
+          countStops = sum(mutationInfoWhole=="Stop")
+          if(seqWithStops==2 & countStops==0) mutationInfo = NA
+          if(seqWithStops==3 & countStops>0) mutationInfo = NA
+        }         
+        
+        if(any(!is.na(mutationInfo))){
+          #Filter mutations based on multipleMutation
+          if(multipleMutation==1 & !is.na(mutationInfo)){
+            mutationCodons = getCodonNumb(as.numeric(names(mutationInfoWhole)))
+            tableMutationCodons <- table(mutationCodons)
+            codonsWithMultipleMutations <- as.numeric(names(tableMutationCodons[tableMutationCodons>1]))
+            if(any(codonsWithMultipleMutations)){
+              #remove the nucleotide mutations in the codons with multiple mutations
+              mutationInfo <- mutationInfo[!(mutationCodons %in% codonsWithMultipleMutations)]
+              #replace those codons with Ns in the input sequence
+              paramSeq[unlist(lapply(codonsWithMultipleMutations, getCodonNucs))] = "N"
+              #matInput[inputMatrixIndex,1] <<- c2s(paramSeq)
+              inputSeq <- c2s(paramSeq)
+            }
+          }
+          
+          #Filter mutations based on the model
+          if(any(mutationInfo)==T | is.na(any(mutationInfo))){        
+            
+            if(model==1 & !is.na(mutationInfo)){
+              mutationInfo <- mutationInfo[mutationInfo=="S"]
+            }  
+            if(any(mutationInfo)==T | is.na(any(mutationInfo))) return(list(mutationInfo,inputSeq))
+            else return(list(NA,inputSeq))
+          }else{
+            return(list(NA,inputSeq))
+          }
+        }else{
+          return(list(NA,inputSeq))
+        }
+        
+        
+      }else{
+        return(list(NA,inputSeq))
+      }
+    
+    
+    }else{
+      return (list(NA,inputSeq))
+    }    
+  }  
+ 
+  # triMutability Background Count
+  buildMutabilityModel <- function( inputMatrixIndex, model=0 , multipleMutation=0, seqWithStops=0, stopMutations=0){
+    
+    #rowOrigMatInput = matInput[inputMatrixIndex,]    
+    seqGL =  gsub("-", "", matInput[inputMatrixIndex,2])
+    seqInput = gsub("-", "", matInput[inputMatrixIndex,1])    
+    #matInput[inputMatrixIndex,] <<- cbind(seqInput,seqGL)
+    tempInput <- cbind(seqInput,seqGL)
+    seqLength = nchar(seqGL)      
+    list_analyzeMutationsFixed<- analyzeMutationsFixed(tempInput, model, multipleMutation, seqWithStops)
+    mutationCount <- list_analyzeMutationsFixed[[1]]
+    seqInput <- list_analyzeMutationsFixed[[2]]
+    BackgroundMatrix = mutabilityMatrix
+    MutationMatrix = mutabilityMatrix    
+    MutationCountMatrix = mutabilityMatrix    
+    if(!is.na(mutationCount)){
+      if((stopMutations==0 & model==0) | (stopMutations==1 & (sum(mutationCount=="Stop")<length(mutationCount))) | (model==1 & (sum(mutationCount=="S")>0)) ){ 
+                  
+        fivermerStartPos = 1:(seqLength-4)
+        fivemerLength <- length(fivermerStartPos)
+        fivemerGL <- substr(rep(seqGL,length(fivermerStartPos)),(fivermerStartPos),(fivermerStartPos+4))
+        fivemerSeq <- substr(rep(seqInput,length(fivermerStartPos)),(fivermerStartPos),(fivermerStartPos+4))
+    
+        #Background
+        for(fivemerIndex in 1:fivemerLength){
+          fivemer = fivemerGL[fivemerIndex]
+          if(!any(grep("N",fivemer))){
+            fivemerCodonPos = fivemerCodon(fivemerIndex)
+            fivemerReadingFrameCodon = substr(fivemer,fivemerCodonPos[1],fivemerCodonPos[3]) 
+            fivemerReadingFrameCodonInputSeq = substr(fivemerSeq[fivemerIndex],fivemerCodonPos[1],fivemerCodonPos[3])          
+            
+            # All mutations model
+            #if(!any(grep("N",fivemerReadingFrameCodon))){
+              if(model==0){
+                if(stopMutations==0){
+                  if(!any(grep("N",fivemerReadingFrameCodonInputSeq)))
+                    BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + 1)              
+                }else{
+                  if( !any(grep("N",fivemerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(fivemerReadingFrameCodon)!="*" ){
+                    positionWithinCodon = which(fivemerCodonPos==3)#positionsWithinCodon[(fivemerCodonPos[1]%%3)+1]
+                    BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + probNonStopMutations[fivemerReadingFrameCodon,positionWithinCodon])
+                  }
+                }
+              }else{ # Only silent mutations
+                if( !any(grep("N",fivemerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(fivemerReadingFrameCodon)!="*" & translateCodonToAminoAcid(fivemerReadingFrameCodonInputSeq)==translateCodonToAminoAcid(fivemerReadingFrameCodon)){
+                  positionWithinCodon = which(fivemerCodonPos==3)
+                  BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + probSMutations[fivemerReadingFrameCodon,positionWithinCodon])
+                }
+              }
+            #}
+          }
+        }
+        
+        #Mutations
+        if(stopMutations==1) mutationCount = mutationCount[mutationCount!="Stop"]
+        if(model==1) mutationCount = mutationCount[mutationCount=="S"]  
+        mutationPositions = as.numeric(names(mutationCount))
+        mutationCount = mutationCount[mutationPositions>2 & mutationPositions<(seqLength-1)]
+        mutationPositions =  mutationPositions[mutationPositions>2 & mutationPositions<(seqLength-1)]
+        countMutations = 0 
+        for(mutationPosition in mutationPositions){
+          fivemerIndex = mutationPosition-2
+          fivemer = fivemerSeq[fivemerIndex]
+          GLfivemer = fivemerGL[fivemerIndex]
+          fivemerCodonPos = fivemerCodon(fivemerIndex)
+          fivemerReadingFrameCodon = substr(fivemer,fivemerCodonPos[1],fivemerCodonPos[3]) 
+          fivemerReadingFrameCodonGL = substr(GLfivemer,fivemerCodonPos[1],fivemerCodonPos[3])
+          if(!any(grep("N",fivemer)) & !any(grep("N",GLfivemer))){
+            if(model==0){
+                countMutations = countMutations + 1              
+                MutationMatrix[GLfivemer] <- (MutationMatrix[GLfivemer] + 1)
+                MutationCountMatrix[GLfivemer] <- (MutationCountMatrix[GLfivemer] + 1)             
+            }else{
+              if( translateCodonToAminoAcid(fivemerReadingFrameCodonGL)!="*" ){
+                  countMutations = countMutations + 1
+                  positionWithinCodon = which(fivemerCodonPos==3)
+                  glNuc =  substr(fivemerReadingFrameCodonGL,positionWithinCodon,positionWithinCodon)
+                  inputNuc =  substr(fivemerReadingFrameCodon,positionWithinCodon,positionWithinCodon)
+                  MutationMatrix[GLfivemer] <- (MutationMatrix[GLfivemer] + substitution[glNuc,inputNuc])
+                  MutationCountMatrix[GLfivemer] <- (MutationCountMatrix[GLfivemer] + 1)                                    
+              }                
+            }                  
+          }              
+        }
+        
+        seqMutability = MutationMatrix/BackgroundMatrix
+        seqMutability = seqMutability/sum(seqMutability,na.rm=TRUE)
+        #cat(inputMatrixIndex,"\t",countMutations,"\n")
+        return(list("seqMutability"  = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix, "BackgroundMatrix"=BackgroundMatrix))      
+        
+      }        
+    }
+  
+  }  
+  
+  #Returns the codon position containing the middle nucleotide
+  fivemerCodon <- function(fivemerIndex){
+    codonPos = list(2:4,1:3,3:5)
+    fivemerType = fivemerIndex%%3
+    return(codonPos[[fivemerType+1]])
+  }
+
+  #returns probability values for one mutation in codons resulting in R, S or Stop
+  probMutations <- function(typeOfMutation){    
+    matMutationProb <- matrix(0,ncol=3,nrow=125,dimnames=list(words(alphabet = c(NUCLEOTIDES,"N"), length=3),c(1:3)))   
+    for(codon in rownames(matMutationProb)){
+        if( !any(grep("N",codon)) ){
+        for(muPos in 1:3){
+          matCodon = matrix(rep(s2c(codon),3),nrow=3,ncol=3,byrow=T)
+          glNuc = matCodon[1,muPos]
+          matCodon[,muPos] = canMutateTo(glNuc) 
+          substitutionRate = substitution[glNuc,matCodon[,muPos]]
+          typeOfMutations = apply(rbind(rep(codon,3),apply(matCodon,1,c2s)),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})        
+          matMutationProb[codon,muPos] <- sum(substitutionRate[typeOfMutations==typeOfMutation])
+        }
+      }
+    }
+    
+    return(matMutationProb) 
+  }
+  
+  
+  
+  
+#Mapping Trinucleotides to fivemers
+mapTriToFivemer <- function(triMutability=triMutability_Literature_Human){
+  rownames(triMutability) <- triMutability_Names
+  Fivemer<-rep(NA,1024)
+  names(Fivemer)<-words(alphabet=NUCLEOTIDES,length=5)
+  Fivemer<-sapply(names(Fivemer),function(Word)return(sum( c(triMutability[substring(Word,3,5),1],triMutability[substring(Word,2,4),2],triMutability[substring(Word,1,3),3]),na.rm=TRUE)))
+  Fivemer<-Fivemer/sum(Fivemer)
+  return(Fivemer)
+}
+
+collapseFivemerToTri<-function(Fivemer,Weights=MutabilityWeights,position=1,NUC="A"){
+  Indices<-substring(names(Fivemer),3,3)==NUC
+  Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
+  tapply(which(Indices),Factors,function(i)weighted.mean(Fivemer[i],Weights[i],na.rm=TRUE))
+}
+
+
+
+CountFivemerToTri<-function(Fivemer,Weights=MutabilityWeights,position=1,NUC="A"){
+  Indices<-substring(names(Fivemer),3,3)==NUC
+  Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
+  tapply(which(Indices),Factors,function(i)sum(Weights[i],na.rm=TRUE))
+}
+
+#Uses the real counts of the mutated fivemers
+CountFivemerToTri2<-function(Fivemer,Counts=MutabilityCounts,position=1,NUC="A"){
+  Indices<-substring(names(Fivemer),3,3)==NUC
+  Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
+  tapply(which(Indices),Factors,function(i)sum(Counts[i],na.rm=TRUE))
+}
+
+bootstrap<-function(x=c(33,12,21),M=10000,alpha=0.05){
+N<-sum(x)
+if(N){
+p<-x/N
+k<-length(x)-1
+tmp<-rmultinom(M, size = N, prob=p)
+tmp_p<-apply(tmp,2,function(y)y/N)
+(apply(tmp_p,1,function(y)quantile(y,c(alpha/2/k,1-alpha/2/k))))
+}
+else return(matrix(0,2,length(x)))
+}
+
+
+
+
+bootstrap2<-function(x=c(33,12,21),n=10,M=10000,alpha=0.05){
+
+N<-sum(x)
+k<-length(x)
+y<-rep(1:k,x)
+tmp<-sapply(1:M,function(i)sample(y,n))
+if(n>1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[,j]==i)))/n
+if(n==1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[j]==i)))/n
+(apply(tmp_p,1,function(z)quantile(z,c(alpha/2/(k-1),1-alpha/2/(k-1)))))
+}
+
+
+
+p_value<-function(x=c(33,12,21),M=100000,x_obs=c(2,5,3)){
+n=sum(x_obs)
+N<-sum(x)
+k<-length(x)
+y<-rep(1:k,x)
+tmp<-sapply(1:M,function(i)sample(y,n))
+if(n>1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[,j]==i)))
+if(n==1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[j]==i)))
+tmp<-rbind(sapply(1:3,function(i)sum(tmp_p[i,]>=x_obs[i])/M),
+sapply(1:3,function(i)sum(tmp_p[i,]<=x_obs[i])/M))
+sapply(1:3,function(i){if(tmp[1,i]>=tmp[2,i])return(-tmp[2,i])else return(tmp[1,i])})
+}
+
+#"D:\\Sequences\\IMGT Germlines\\Human_SNPless_IGHJ.FASTA"
+# Remove SNPs from IMGT germline segment alleles
+generateUnambiguousRepertoire <- function(repertoireInFile,repertoireOutFile){
+  repertoireIn <- read.fasta(repertoireInFile, seqtype="DNA",as.string=T,set.attributes=F,forceDNAtolower=F)
+  alleleNames <- sapply(names(repertoireIn),function(x)strsplit(x,"|",fixed=TRUE)[[1]][2])
+  SNPs <- tapply(repertoireIn,sapply(alleleNames,function(x)strsplit(x,"*",fixed=TRUE)[[1]][1]),function(x){
+    Indices<-NULL
+    for(i in 1:length(x)){
+      firstSeq = s2c(x[[1]])
+      iSeq = s2c(x[[i]])
+      Indices<-c(Indices,which(firstSeq[1:320]!=iSeq[1:320] & firstSeq[1:320]!="." & iSeq[1:320]!="."  ))
+    }
+    return(sort(unique(Indices)))
+  })
+ repertoireOut <- repertoireIn
+ repertoireOut <- lapply(names(repertoireOut), function(repertoireName){
+                                        alleleName <- strsplit(repertoireName,"|",fixed=TRUE)[[1]][2]
+                                        geneSegmentName <- strsplit(alleleName,"*",fixed=TRUE)[[1]][1]
+                                        alleleSeq <- s2c(repertoireOut[[repertoireName]])
+                                        alleleSeq[as.numeric(unlist(SNPs[geneSegmentName]))] <- "N"
+                                        alleleSeq <- c2s(alleleSeq)
+                                        repertoireOut[[repertoireName]] <- alleleSeq
+                                      })
+  names(repertoireOut) <- names(repertoireIn)
+  write.fasta(repertoireOut,names(repertoireOut),file.out=repertoireOutFile)                                               
+                                      
+}
+
+
+
+
+
+
+############
+groupBayes2 = function(indexes, param_resultMat){
+  
+  BayesGDist_Focused_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2]+x[4])}))
+  BayesGDist_Focused_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(3,2,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[3]+x[2]+x[4])}))
+  #BayesGDist_Local_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2])}))
+  #BayesGDist_Local_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[3]+x[4])}))
+  #BayesGDist_Global_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2,3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2]+x[3]+x[4])}))
+  #BayesGDist_Global_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(1,2,3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[1]+x[2]+x[3]+x[4])}))
+  return ( list("BayesGDist_Focused_CDR"=BayesGDist_Focused_CDR,
+                "BayesGDist_Focused_FWR"=BayesGDist_Focused_FWR) )
+                #"BayesGDist_Local_CDR"=BayesGDist_Local_CDR,
+                #"BayesGDist_Local_FWR" = BayesGDist_Local_FWR))
+#                "BayesGDist_Global_CDR" = BayesGDist_Global_CDR,
+#                "BayesGDist_Global_FWR" = BayesGDist_Global_FWR) )
+
+
+}
+
+
+calculate_bayesG <- function( x=array(), N=array(), p=array(), max_sigma=20, length_sigma=4001){
+  G <- max(length(x),length(N),length(p))
+  x=array(x,dim=G)
+  N=array(N,dim=G)
+  p=array(p,dim=G)
+
+  indexOfZero = N>0 & p>0
+  N = N[indexOfZero]
+  x = x[indexOfZero]
+  p = p[indexOfZero]  
+  G <- length(x)
+  
+  if(G){
+    
+    cons<-array( dim=c(length_sigma,G) )
+    if(G==1) {
+    return(calculate_bayes(x=x[G],N=N[G],p=p[G],max_sigma=max_sigma,length_sigma=length_sigma))
+    }
+    else {
+      for(g in 1:G) cons[,g] <- calculate_bayes(x=x[g],N=N[g],p=p[g],max_sigma=max_sigma,length_sigma=length_sigma)
+      listMatG <- convolutionPowersOfTwoByTwos(cons,length_sigma=length_sigma)
+      y<-calculate_bayesGHelper(listMatG,length_sigma=length_sigma)
+      return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
+    }
+  }else{
+    return(NA)
+  }
+}
+
+
+calculate_bayesGHelper <- function( listMatG,length_sigma=4001 ){
+  matG <- listMatG[[1]]  
+  groups <- listMatG[[2]]
+  i = 1  
+  resConv <- matG[,i]
+  denom <- 2^groups[i]
+  if(length(groups)>1){
+    while( i<length(groups) ){
+      i = i + 1
+      resConv <- weighted_conv(resConv, matG[,i], w= {{2^groups[i]}/denom} ,length_sigma=length_sigma)
+      #cat({{2^groups[i]}/denom},"\n")
+      denom <- denom + 2^groups[i]
+    }
+  }
+  return(resConv)  
+}
+
+weighted_conv<-function(x,y,w=1,m=100,length_sigma=4001){
+lx<-length(x)
+ly<-length(y)
+if({lx<m}| {{lx*w}<m}| {{ly}<m}| {{ly*w}<m}){
+if(w<1){
+y1<-approx(1:ly,y,seq(1,ly,length.out=m))$y
+x1<-approx(1:lx,x,seq(1,lx,length.out=m/w))$y
+lx<-length(x1)
+ly<-length(y1)
+}
+else {
+y1<-approx(1:ly,y,seq(1,ly,length.out=m*w))$y
+x1<-approx(1:lx,x,seq(1,lx,length.out=m))$y
+lx<-length(x1)
+ly<-length(y1)
+}
+}
+else{
+x1<-x
+y1<-approx(1:ly,y,seq(1,ly,length.out=floor(lx*w)))$y
+ly<-length(y1)
+}
+tmp<-approx(x=1:(lx+ly-1),y=convolve(x1,rev(y1),type="open"),xout=seq(1,lx+ly-1,length.out=length_sigma))$y
+tmp[tmp<=0] = 0 
+return(tmp/sum(tmp))
+}
+
+########################
+
+
+
+
+mutabilityMatrixONE<-rep(0,4)
+names(mutabilityMatrixONE)<-NUCLEOTIDES
+
+  # triMutability Background Count
+  buildMutabilityModelONE <- function( inputMatrixIndex, model=0 , multipleMutation=0, seqWithStops=0, stopMutations=0){
+    
+    #rowOrigMatInput = matInput[inputMatrixIndex,]    
+    seqGL =  gsub("-", "", matInput[inputMatrixIndex,2])
+    seqInput = gsub("-", "", matInput[inputMatrixIndex,1])    
+    matInput[inputMatrixIndex,] <<- c(seqInput,seqGL)
+    seqLength = nchar(seqGL)      
+    mutationCount <- analyzeMutations(inputMatrixIndex, model, multipleMutation, seqWithStops)
+    BackgroundMatrix = mutabilityMatrixONE
+    MutationMatrix = mutabilityMatrixONE    
+    MutationCountMatrix = mutabilityMatrixONE    
+    if(!is.na(mutationCount)){
+      if((stopMutations==0 & model==0) | (stopMutations==1 & (sum(mutationCount=="Stop")<length(mutationCount))) | (model==1 & (sum(mutationCount=="S")>0)) ){ 
+                  
+#         ONEmerStartPos = 1:(seqLength)
+#         ONEmerLength <- length(ONEmerStartPos)
+        ONEmerGL <- s2c(seqGL)
+        ONEmerSeq <- s2c(seqInput)
+    
+        #Background
+        for(ONEmerIndex in 1:seqLength){
+          ONEmer = ONEmerGL[ONEmerIndex]
+          if(ONEmer!="N"){
+            ONEmerCodonPos = getCodonPos(ONEmerIndex)
+            ONEmerReadingFrameCodon = c2s(ONEmerGL[ONEmerCodonPos]) 
+            ONEmerReadingFrameCodonInputSeq = c2s(ONEmerSeq[ONEmerCodonPos] )         
+            
+            # All mutations model
+            #if(!any(grep("N",ONEmerReadingFrameCodon))){
+              if(model==0){
+                if(stopMutations==0){
+                  if(!any(grep("N",ONEmerReadingFrameCodonInputSeq)))
+                    BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + 1)              
+                }else{
+                  if( !any(grep("N",ONEmerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)!="*"){
+                    positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)#positionsWithinCodon[(ONEmerCodonPos[1]%%3)+1]
+                    BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + probNonStopMutations[ONEmerReadingFrameCodon,positionWithinCodon])
+                  }
+                }
+              }else{ # Only silent mutations
+                if( !any(grep("N",ONEmerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)!="*" & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)==translateCodonToAminoAcid(ONEmerReadingFrameCodon) ){
+                  positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)
+                  BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + probSMutations[ONEmerReadingFrameCodon,positionWithinCodon])
+                }
+              }
+            }
+          }
+        }
+        
+        #Mutations
+        if(stopMutations==1) mutationCount = mutationCount[mutationCount!="Stop"]
+        if(model==1) mutationCount = mutationCount[mutationCount=="S"]  
+        mutationPositions = as.numeric(names(mutationCount))
+        mutationCount = mutationCount[mutationPositions>2 & mutationPositions<(seqLength-1)]
+        mutationPositions =  mutationPositions[mutationPositions>2 & mutationPositions<(seqLength-1)]
+        countMutations = 0 
+        for(mutationPosition in mutationPositions){
+          ONEmerIndex = mutationPosition
+          ONEmer = ONEmerSeq[ONEmerIndex]
+          GLONEmer = ONEmerGL[ONEmerIndex]
+          ONEmerCodonPos = getCodonPos(ONEmerIndex)
+          ONEmerReadingFrameCodon = c2s(ONEmerSeq[ONEmerCodonPos])  
+          ONEmerReadingFrameCodonGL =c2s(ONEmerGL[ONEmerCodonPos])  
+          if(!any(grep("N",ONEmer)) & !any(grep("N",GLONEmer))){
+            if(model==0){
+                countMutations = countMutations + 1              
+                MutationMatrix[GLONEmer] <- (MutationMatrix[GLONEmer] + 1)
+                MutationCountMatrix[GLONEmer] <- (MutationCountMatrix[GLONEmer] + 1)             
+            }else{
+              if( translateCodonToAminoAcid(ONEmerReadingFrameCodonGL)!="*" ){
+                  countMutations = countMutations + 1
+                  positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)
+                  glNuc =  substr(ONEmerReadingFrameCodonGL,positionWithinCodon,positionWithinCodon)
+                  inputNuc =  substr(ONEmerReadingFrameCodon,positionWithinCodon,positionWithinCodon)
+                  MutationMatrix[GLONEmer] <- (MutationMatrix[GLONEmer] + substitution[glNuc,inputNuc])
+                  MutationCountMatrix[GLONEmer] <- (MutationCountMatrix[GLONEmer] + 1)                                    
+              }                
+            }                  
+          }              
+        }
+        
+        seqMutability = MutationMatrix/BackgroundMatrix
+        seqMutability = seqMutability/sum(seqMutability,na.rm=TRUE)
+        #cat(inputMatrixIndex,"\t",countMutations,"\n")
+        return(list("seqMutability"  = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix, "BackgroundMatrix"=BackgroundMatrix))      
+#         tmp<-list("seqMutability"  = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix)
+      }        
+    }
+  
+################
+# $Id: trim.R 989 2006-10-29 15:28:26Z ggorjan $
+
+trim <- function(s, recode.factor=TRUE, ...)
+  UseMethod("trim", s)
+
+trim.default <- function(s, recode.factor=TRUE, ...)
+  s
+
+trim.character <- function(s, recode.factor=TRUE, ...)
+{
+  s <- sub(pattern="^ +", replacement="", x=s)
+  s <- sub(pattern=" +$", replacement="", x=s)
+  s
+}
+
+trim.factor <- function(s, recode.factor=TRUE, ...)
+{
+  levels(s) <- trim(levels(s))
+  if(recode.factor) {
+    dots <- list(x=s, ...)
+    if(is.null(dots$sort)) dots$sort <- sort
+    s <- do.call(what=reorder.factor, args=dots)
+  }
+  s
+}
+
+trim.list <- function(s, recode.factor=TRUE, ...)
+  lapply(s, trim, recode.factor=recode.factor, ...)
+
+trim.data.frame <- function(s, recode.factor=TRUE, ...)
+{
+  s[] <- trim.list(s, recode.factor=recode.factor, ...)
+  s
+}
+#######################################
+# Compute the expected for each sequence-germline pair by codon 
+getExpectedIndividualByCodon <- function(matInput){    
+if( any(grep("multicore",search())) ){  
+  facGL <- factor(matInput[,2])
+  facLevels = levels(facGL)
+  LisGLs_MutabilityU = mclapply(1:length(facLevels),  function(x){
+    computeMutabilities(facLevels[x])
+  })
+  facIndex = match(facGL,facLevels)
+  
+  LisGLs_Mutability = mclapply(1:nrow(matInput),  function(x){
+    cInput = rep(NA,nchar(matInput[x,1]))
+    cInput[s2c(matInput[x,1])!="N"] = 1
+    LisGLs_MutabilityU[[facIndex[x]]] * cInput                                                   
+  })
+  
+  LisGLs_Targeting =  mclapply(1:dim(matInput)[1],  function(x){
+    computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
+  })
+  
+  LisGLs_MutationTypes  = mclapply(1:length(matInput[,2]),function(x){
+    #print(x)
+    computeMutationTypes(matInput[x,2])
+  })
+  
+  LisGLs_R_Exp = mclapply(1:nrow(matInput),  function(x){
+    Exp_R <-  rollapply(as.zoo(1:readEnd),width=3,by=3,
+                        function(codonNucs){                                                      
+                          RPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="R") 
+                          sum( LisGLs_Targeting[[x]][,codonNucs][RPos], na.rm=T ) 
+                        }
+    )                                                   
+  })
+  
+  LisGLs_S_Exp = mclapply(1:nrow(matInput),  function(x){
+    Exp_S <-  rollapply(as.zoo(1:readEnd),width=3,by=3,
+                        function(codonNucs){                                                      
+                          SPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="S")   
+                          sum( LisGLs_Targeting[[x]][,codonNucs][SPos], na.rm=T )
+                        }
+    )                                                 
+  })                                                
+  
+  Exp_R = matrix(unlist(LisGLs_R_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)  
+  Exp_S = matrix(unlist(LisGLs_S_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)  
+  return( list( "Expected_R"=Exp_R, "Expected_S"=Exp_S) )
+  }else{
+    facGL <- factor(matInput[,2])
+    facLevels = levels(facGL)
+    LisGLs_MutabilityU = lapply(1:length(facLevels),  function(x){
+      computeMutabilities(facLevels[x])
+    })
+    facIndex = match(facGL,facLevels)
+    
+    LisGLs_Mutability = lapply(1:nrow(matInput),  function(x){
+      cInput = rep(NA,nchar(matInput[x,1]))
+      cInput[s2c(matInput[x,1])!="N"] = 1
+      LisGLs_MutabilityU[[facIndex[x]]] * cInput                                                   
+    })
+    
+    LisGLs_Targeting =  lapply(1:dim(matInput)[1],  function(x){
+      computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
+    })
+    
+    LisGLs_MutationTypes  = lapply(1:length(matInput[,2]),function(x){
+      #print(x)
+      computeMutationTypes(matInput[x,2])
+    })
+    
+    LisGLs_R_Exp = lapply(1:nrow(matInput),  function(x){
+      Exp_R <-  rollapply(as.zoo(1:readEnd),width=3,by=3,
+                          function(codonNucs){                                                      
+                            RPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="R") 
+                            sum( LisGLs_Targeting[[x]][,codonNucs][RPos], na.rm=T ) 
+                          }
+      )                                                   
+    })
+    
+    LisGLs_S_Exp = lapply(1:nrow(matInput),  function(x){
+      Exp_S <-  rollapply(as.zoo(1:readEnd),width=3,by=3,
+                          function(codonNucs){                                                      
+                            SPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="S")   
+                            sum( LisGLs_Targeting[[x]][,codonNucs][SPos], na.rm=T )
+                          }
+      )                                                 
+    })                                                
+    
+    Exp_R = matrix(unlist(LisGLs_R_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)  
+    Exp_S = matrix(unlist(LisGLs_S_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)  
+    return( list( "Expected_R"=Exp_R, "Expected_S"=Exp_S) )    
+  }
+}
+
+# getObservedMutationsByCodon <- function(listMutations){
+#   numbSeqs <- length(listMutations) 
+#   obsMu_R <- matrix(0,nrow=numbSeqs,ncol=readEnd/3,dimnames=list(c(1:numbSeqs),c(1:(readEnd/3))))
+#   obsMu_S <- obsMu_R
+#   temp <- mclapply(1:length(listMutations), function(i){
+#     arrMutations = listMutations[[i]]
+#     RPos = as.numeric(names(arrMutations)[arrMutations=="R"])
+#     RPos <- sapply(RPos,getCodonNumb)                                                                    
+#     if(any(RPos)){
+#       tabR <- table(RPos)
+#       obsMu_R[i,as.numeric(names(tabR))] <<- tabR
+#     }                                    
+#     
+#     SPos = as.numeric(names(arrMutations)[arrMutations=="S"])
+#     SPos <- sapply(SPos,getCodonNumb)
+#     if(any(SPos)){
+#       tabS <- table(SPos)
+#       obsMu_S[i,names(tabS)] <<- tabS
+#     }                                          
+#   }
+#   )
+#   return( list( "Observed_R"=obsMu_R, "Observed_S"=obsMu_S) ) 
+# }
+
+getObservedMutationsByCodon <- function(listMutations){
+  numbSeqs <- length(listMutations) 
+  obsMu_R <- matrix(0,nrow=numbSeqs,ncol=readEnd/3,dimnames=list(c(1:numbSeqs),c(1:(readEnd/3))))
+  obsMu_S <- obsMu_R
+  temp <- lapply(1:length(listMutations), function(i){
+    arrMutations = listMutations[[i]]
+    RPos = as.numeric(names(arrMutations)[arrMutations=="R"])
+    RPos <- sapply(RPos,getCodonNumb)                                                                    
+    if(any(RPos)){
+      tabR <- table(RPos)
+      obsMu_R[i,as.numeric(names(tabR))] <<- tabR
+    }                                    
+    
+    SPos = as.numeric(names(arrMutations)[arrMutations=="S"])
+    SPos <- sapply(SPos,getCodonNumb)
+    if(any(SPos)){
+      tabS <- table(SPos)
+      obsMu_S[i,names(tabS)] <<- tabS
+    }                                          
+  }
+  )
+  return( list( "Observed_R"=obsMu_R, "Observed_S"=obsMu_S) ) 
+}
+