SIMPAL - SIMple Plasmid ALignemnt

SIMPAL is a small script written in R with the purpose for simple N->1 gene alignment and mutation analysis. The idea is to have a handy tool to analyze some tens to thousands of equivalent data sets containing just one or a few mutation as it is the situation within a gene library derived from random mutagenesis protocols. The correct ORF is determined by a unique identifier (Fig. 1) , a short recognition sequence of 12-15 nucleotides upstream of the target gene which also serves as the starting point for further alignment. Before you can start you need to tell the script just three information. First of course the wild-type sequence within a certain plasmid system (intended) or even as stand-alone gene both written as plain text into a single file. Secondly the unique identifier also given as a plain text within a another file. The third information is the directory where the script could find your sequencing library in FASTA-format. Assuming everything works well, the output of the script contains some general statistics about the absolute number and distribution of single and multiple mutations. Moreover, a so called level plot and its corresponding Nx20 matrix are generated representing a mutation frequency diagram with a amino acid to amino acid assignment. The results from a mutation analysis of up to several thousands of data sets is thus concentrated within a single plot (Fig. 2).

SIMPAL working principle — Fig1. Principle of sequence alignment and analysis
using SIMPAL (Graphic made with Inkscape).

Example of mutation level plot — Fig2. Example of mutation analysis along 96 data sets (Output from R script).

SIMPAL source code

names<-c(list.files(path = "C:/.../SIMPAL/Samples/"))
data=0
for (i in 1:length(names)){
   file<-read.table(paste("C:/.../SIMPAL/Samples/",names[i], sep=""),colClasses = "character",skip="1",nrows=1,header=FALSE)
   s1<-as.character(file)
   s2<-c(unlist(strsplit(s1,"")))
   s3<-toupper(s2) #Großbuchstaben 
   data<-list(data,s3)
 }

1    findwords <- function(tf) {
2       # read in the words from the file, into a vector of mode character
3       txt <- scan(tf,"")
4       wl <- list()
5       for (i in 1:length(txt)) {
6          wrd <- txt[i]  # ith word in input file
7          wl[[wrd]] <- c(wl[[wrd]],i)
8       }
9       return(wl)
10    }


#MasterFASTA
mf<-read.table("C:/.../SIMPAL/MasterFASTA.txt",colClasses = "character",nrows=1,header=FALSE) #Einlesen, nur eine Zeile
mf2<-as.character(mf)
fasta<-c(unlist(strsplit(mf2,""))) #in Einzelbuchstaben zerlegen und als Einzelelemte mit "unlist" zuordnen

#MasterProtein
mp<-read.table("C:/AG_Möglich/Directed_Evolution/Simpal/MasterPROT.txt",colClasses = "character",nrows=1,header=FALSE)
mp2<-as.character(mp)
prot<-c(unlist(strsplit(mp2,"")))

#Code-Sun
cs<-read.table("C:/.../SIMPAL/CodeSun.txt",header=FALSE)

#Adress Code-Sun mxn:=64x4
x<-as.character(cs[50,2]) #One letter of postion 50 
x<-as.character(cs[50,3]) #Three letter of portion 50
x<-as.character(cs[50,1]) #Codon of portion 50

#Samples
s2<-read.table("C:/.../SIMPAL/Samples/2.fas",colClasses = "character",skip="1",nrows=1,header=FALSE)
s22<-as.character(s2)
samp22<-c(unlist(strsplit(s22,"")))
samp2=toupper(samp22) 

s1<-read.table("C:/.../SIMPAL/Samples/1.fas",colClasses = "character",skip="1",nrows=1,header=FALSE)
s12<-as.character(s1)
samp11<-c(unlist(strsplit(s12,"")))
samp1=toupper(samp11)

#Compare codons
cod.compare<-function(cs,f){
a=0
  for(i in 1:dim(cs)[1]){
   if(f==as.character(cs[i,1])){a=as.character(cs[i,3])
  }}
  return(a)
}
cod.compare(cs,x[3,1])

#Set identifier
n="GTGGTGAATGTG"
i<-as.character(n)
id<-c(unlist(strsplit(i,"")))

#Find identifier
find<-function(samp2,id){
pos=0
  for(i in 1:length(samp2)){
    if(samp2[i]==id[1] & samp2[i+1]==id[2] & samp2[i+2]==id[3] & samp2[i+3]==id[4] & samp2[i+4]==id[5] & samp2[i+5]==id[6] & samp2[i+6]==id[7]){
      pos=i
      }}
    if(pos!=0){return(pos)}
     else{return("No hit")}
}
find(samp2,id)   
n<-find(samp2,id) 
samp2<-samp2[-1:-n+1] 
samp2<-samp2[1:length(fasta)]


  #Problems with TRUE/FALSE --> need to start from nucleotide 2
  find<-function(samp1,id){
  pos=0
    for(i in 1:length(samp1)){
      if(samp1[i]==id[1] & samp1[i+1]==id[2] & samp1[i+2]==id[3] & samp1[i+3]==id[4] & samp1[i+4]==id[5] & samp1[i+5]==id[6] & samp1[i+7]==id[7]){
        pos=i
        }}
      if(pos!=0){return(pos)}
       else{return("No hit")}
  }
  n<-find(samp1,id) 
  samp1<-samp1[-1:-n+1]
  samp1<-samp1[1:length(fasta)]


#Compare Nucleotide by Nucleotide 

nbn.compare<-function(samp1,fasta){
nucF=0
nucS=0
pos=0
mat=0
 for(i in 1:length(fasta)){
   if(fasta[i]!=samp1[i]){
     nucF<-rbind(nucF,fasta[i])
     nucS<-rbind(nucS,samp1[i])
     pos<-rbind(pos,i)
   }
  }
 mat<-cbind(nucF,nucS,pos)
 return(mat)
}
nbn.compare(samp1,fasta)


#Join Strings 
f<-c(e1,e2,e3,...)
g=paste(f, collapse="")

g=paste(f<-c(a,b,c), collapse="")


#Erweiterung nbn.compare (x%%y!=0 modulo not zero) 


nbn.compare<-function(samp2,fasta){
nucF=0
nucS=0
mat=0
pos=0
h=0
a=0
 for(i in 4:length(fasta)){
   if(fasta[i]!=samp2[i]){
    if(i%%3==0){
      a<-paste(h<-c(samp2[i-2],samp2[i-1],samp2[i]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3!=0){
      a<-paste(h<-c(samp2[i-1],samp2[i],samp2[i+1]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3==0){
      a<-paste(h<-c(samp2[i],samp2[i+1],samp2[i+2]), collapse="")
    }  
  nucS<-rbind(nucS,a)
  print(nucS)
  print(length(nucS))
  nucF<-rbind(nucF,fasta[i])
  print(length(nucF))
  pos<-rbind(pos,i) 
  print(length(pos))
  }
 }
 mat<-cbind(nucS,pos,nucF)
 return(mat)
}
nbn.compare(samp2,fasta)


#Erweiterung 2 nbn.compare (+cod.compare) 


nbn.compare<-function(samp2,fasta,cs){
nucF=0
nucS=0
mat=0
posAA=0
posN=0
h=0
a=0
b=0
c=0
d=0
e=0
AAS=0
AAF=0
codF=0
 for(i in 4:length(fasta)){
   if(fasta[i]!=samp2[i]){
    if(i%%3==0){
      a<-paste(h<-c(samp2[i-2],samp2[i-1],samp2[i]), collapse="")
      c=(i/3)-3 #Realpostion abziehen
      d<-paste(h<-c(fasta[i-2],fasta[i-1],fasta[i]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3!=0){
      a<-paste(h<-c(samp2[i-1],samp2[i],samp2[i+1]), collapse="")
      c=((i+1)/3)-3
      d<-paste(h<-c(fasta[i-1],fasta[i],fasta[i+1]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3==0){
      a<-paste(h<-c(samp2[i],samp2[i+1],samp2[i+2]), collapse="")
      c=((i+2)/3)-3 
      d<-paste(h<-c(fasta[i],fasta[i+1],fasta[i+2]), collapse="")
    }  
  nucS<-rbind(nucS,a)
  nucF<-rbind(nucF,d)
  posAA<-rbind(posAA,c) 
  posN<-rbind(posN,i-9)
  for(i in 1:dim(cs)[1]){
   if(a==as.character(cs[i,1])){
    b=as.character(cs[i,2])
    AAS<-rbind(AAS,b)
   }
   if(d==as.character(cs[i,1])){
    e=as.character(cs[i,2])
    AAF<-rbind(AAF,e)
   }
  }
  }
 }
 mat<-cbind(posN,nucF,AAF,posAA,AAS,nucS)
 return(mat)
}
nbn.compare(samp2,fasta,cs)


#Automatization 

mf<-read.table("C:/.../SIMPAL/MasterFASTA.txt",colClasses = "character",nrows=1,header=FALSE) #Einlesen, nur eine Zeile
mf2<-as.character(mf)
fasta<-c(unlist(strsplit(mf2,""))) #in Einzelbuchstaben zerlegen und als Einzelelemte mit "unlist" zuordnen


mp<-read.table("C:/.../SIMPAL/MasterPROT.txt",colClasses = "character",nrows=1,header=FALSE)
mp2<-as.character(mp)
prot<-c(unlist(strsplit(mp2,"")))


cs<-read.table("C:/.../SIMPAL/CodeSun.txt",header=FALSE)


s2<-read.table("C:/.../SIMPAL/Samples/2.fas",colClasses = "character",skip="1",nrows=1,header=FALSE)
s22<-as.character(s2)
samp22<-c(unlist(strsplit(s22,"")))
samp2=toupper(samp22) 

s1<-read.table("C:/.../SIMPAL/Samples/1.fas",colClasses = "character",skip="1",nrows=1,header=FALSE)
s12<-as.character(s1)
samp11<-c(unlist(strsplit(s12,"")))
samp1=toupper(samp11)



n="GTGGTGAATGTG"
i<-as.character(n)
id<-c(unlist(strsplit(i,"")))


find<-function(samp2,id){
pos=0
  for(i in 1:length(samp2)){
    if(samp2[i]==id[1] & samp2[i+1]==id[2] & samp2[i+2]==id[3] & samp2[i+3]==id[4] & samp2[i+4]==id[5] & samp2[i+5]==id[6] & samp2[i+6]==id[7]){
      pos=i
      }}
    if(pos!=0){return(pos)}
     else{return("No hit")}
}
find(samp2,id)   
n<-find(samp2,id) 
samp2<-samp2[-1:-n+1] 
samp2<-samp2[1:length(fasta)]


nbn.compare<-function(samp2,fasta,cs){
nucF=0
nucS=0
mat=0
posAA=0
posN=0
h=0
a=0
b=0
c=0
d=0
e=0
AAS=0
AAF=0
codF=0
 for(i in 4:length(fasta)){
   if(fasta[i]!=samp2[i]){
    if(i%%3==0){
      a<-paste(h<-c(samp2[i-2],samp2[i-1],samp2[i]), collapse="")
      c=(i/3)-3 #Realpostion abziehen
      d<-paste(h<-c(fasta[i-2],fasta[i-1],fasta[i]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3!=0){
      a<-paste(h<-c(samp2[i-1],samp2[i],samp2[i+1]), collapse="")
      c=((i+1)/3)-3
      d<-paste(h<-c(fasta[i-1],fasta[i],fasta[i+1]), collapse="")
    }
    if(i%%3!=0 & (i-1)%%3==0){
      a<-paste(h<-c(samp2[i],samp2[i+1],samp2[i+2]), collapse="")
      c=((i+2)/3)-3 
      d<-paste(h<-c(fasta[i],fasta[i+1],fasta[i+2]), collapse="")
    }  
   nucS<-rbind(nucS,a)
   nucF<-rbind(nucF,d)
   posAA<-rbind(posAA,c) 
   posN<-rbind(posN,i-9)
  for(i in 1:dim(cs)[1]){
   if(a==as.character(cs[i,1])){
    b=as.character(cs[i,2])
    AAS<-rbind(AAS,b)
   }
   if(d==as.character(cs[i,1])){
    e=as.character(cs[i,2])
    AAF<-rbind(AAF,e)
   }
  }
  }
 }
 mat<-cbind(posN,nucF,AAF,posAA,AAS,nucS)
 return(mat)
}
nbn.compare(samp2,fasta,cs)