# Lab 9 for the University of Tulsa's CS-6643 Bioinformatics Course
# Pairwise Sequence Alignment with Dynamic Programming
# Professor: Dr. McKinney, Fall 2022
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


#### Part A: Specifying the Input
## Score Rules and Seqs
x_str <- "ATAC"   # side sequence
y_str <- "GTGTAC" # top sequence
match_score <- 3
mismatch_score <- -1
gap_penalty <- -4

## Substitution Matrix
dna.letters<-c("A","C","G","T")
num.letters <- length(dna.letters)
S<-data.frame(matrix(0,nrow=num.letters,ncol=num.letters))  # data frame
rownames(S)<-dna.letters; colnames(S)<-dna.letters
for (i in 1:4){
  for (j in 1:4){
    if(dna.letters[i]==dna.letters[j]){ 
      S[i,j]<- match_score
    }
    else{ 
      S[i,j]<- mismatch_score
    }
  }
}

#### Part B: Alignment Score Matrix (F) and Traceback Matrix (T)
x <- unlist(strsplit(x_str, ""))
y <- unlist(strsplit(y_str, ""))
x.len <- length(x) 
y.len <- length(y) 

Fmat<-matrix(0,nrow=x.len+1,ncol=y.len+1) 
Tmat<-Fmat # 0's to start

rownames(Fmat)<-c("-",x); colnames(Fmat)<-c("-",y)
rownames(Tmat)<-c("-",x); colnames(Tmat)<-c("-",y)

# create first row and column
Fmat[,1]<- seq(from=0,len=x.len+1,by=-abs(gap_penalty))  
Fmat[1,]<- seq(from=0,len=y.len+1,by=-abs(gap_penalty))
Tmat[,1]<- rep(2,x.len+1)  # 2 means align with a gap in the upper seq
Tmat[1,]<- rep(3,y.len+1)  # 3 means align with a gap in the side seq


#### Part C: Building Fmat and Tmat
## Building Fmat and Tmat
for (i in 2:nrow(Fmat)){
  for (j in 2:ncol(Fmat)){    # use F recursive rules 
    test_three_cases <- c(Fmat[i-1, j-1] + S[rownames(Fmat)[i], colnames(Fmat)[j]],   # 1 mis/match
                          Fmat[i-1, j] + gap_penalty,    # 2 up-gap
                          Fmat[i, j-1] + gap_penalty)  # 3 left-gap
    Fmat[i,j]=max(test_three_cases)
    Tmat[i,j]=which.max(test_three_cases)
  }
}
final_score <- Fmat[nrow(Fmat),ncol(Fmat)]

## Aligning from Tmat
n <- nrow(Tmat) 
m <- ncol(Tmat)
seq_align <- character()
while( (n+m)!=2 ){
  if (Tmat[n,m]==1){
    curr_align_col <- rbind(x[n-1],y[m-1]) 
    seq_align <- cbind(curr_align_col,seq_align)
    n <- n-1
    m <- m-1
  }else if(Tmat[n,m]==2){ 
    curr_align_col <- rbind(x[n-1],"-")
    seq_align <- cbind(curr_align_col,seq_align)
    n <- n-1
  }else{
    curr_align_col <- rbind("-",y[m-1])
    seq_align <- cbind(curr_align_col,seq_align)
    m <- m-1
  }           
} # end while
seq_align


#### Part D: Convert to functions
make.alignment.matrices <- function(x_str, y_str, match_score, mismatch_score, 
                                    gap_penalty){
  ## Substitution Matrix
  dna.letters<-c("A","C","G","T")
  num.letters <- length(dna.letters)
  S<-data.frame(matrix(0,nrow=num.letters,ncol=num.letters))  # data frame
  rownames(S)<-dna.letters; colnames(S)<-dna.letters
  for (i in 1:4){
    for (j in 1:4){
      if(dna.letters[i]==dna.letters[j]){ 
        S[i,j]<- match_score
      }
      else{ 
        S[i,j]<- mismatch_score
      }
    }
  }
  
  ## F Matrix and T Matrix
  x <- unlist(strsplit(x_str, ""))
  y <- unlist(strsplit(y_str, ""))
  x.len <- length(x) 
  y.len <- length(y) 
  
  Fmat<-matrix(0,nrow=x.len+1,ncol=y.len+1) 
  Tmat<-Fmat # 0's to start
  
  rownames(Fmat)<-c("-",x); colnames(Fmat)<-c("-",y)
  rownames(Tmat)<-c("-",x); colnames(Tmat)<-c("-",y)
  
  # create first row and column
  Fmat[,1]<- seq(from=0,len=x.len+1,by=-abs(gap_penalty))  
  Fmat[1,]<- seq(from=0,len=y.len+1,by=-abs(gap_penalty))
  Tmat[,1]<- rep(2,x.len+1)  # 2 means align with a gap in the upper seq
  Tmat[1,]<- rep(3,y.len+1)  # 3 means align with a gap in the side seq
  
  ## Building Fmat and Tmat
  for (i in 2:nrow(Fmat)){
    for (j in 2:ncol(Fmat)){    # use F recursive rules 
      test_three_cases <- c(Fmat[i-1, j-1] + S[rownames(Fmat)[i], colnames(Fmat)[j]],   # 1 mis/match
                            Fmat[i-1, j] + gap_penalty,    # 2 up-gap
                            Fmat[i, j-1] + gap_penalty)  # 3 left-gap
      Fmat[i,j]=max(test_three_cases)
      Tmat[i,j]=which.max(test_three_cases)
    }
  }
  final_score <- Fmat[nrow(Fmat),ncol(Fmat)]
  return(list(Fmat=Fmat, Tmat=Tmat, score_out=final_score))
}

# load new input
x_str2 <- "GATTA"  # side sequence
y_str2 <- "GAATTC" # top sequence
match_score <- 2
mismatch_score <- -1
gap_penalty <- -2

align.list2 <- make.alignment.matrices(x_str2, y_str2, match_score,
                                       mismatch_score, gap_penalty)
align.list2$Fmat
align.list2$Tmat
align.list2$score_out  
  
if (!require("gplots")) install.packages("gplots")
library(gplots)
Fmat2 <- align.list2$Fmat
col = c("black","blue","red","yellow","green")
breaks = seq(min(Fmat2),max(Fmat2),len=length(col)+1)

heatmap.2(Fmat2[-1,-1], dendrogram='none', density.info="none", 
          Rowv=FALSE, Colv=FALSE, trace='none', 
          breaks = breaks, col = col,
          sepwidth=c(0.01,0.01),
          sepcolor="black",
          colsep=1:ncol(Fmat2),
          rowsep=1:nrow(Fmat2))

#### Part E: Traceback Matrix
show.alignment <- function(x_str,y_str,Tmat){ 
  ################ create the alignment
  # input Tmat and the two sequences: x side seq and y is top seq
  # make character vectors out of the strings
  x<-unlist(strsplit(x_str,""))
  y<-unlist(strsplit(y_str,""))
  
  n<-nrow(Tmat) # start at bottom right of Tmat
  m<-ncol(Tmat)
  alignment<-character()
  while( (n+m)!=2 ){
    if (Tmat[n,m]==1){
      # subtract 1 from x and y indices because they are
      # one row/col smaller than Tmat
      curr_align_col <- rbind(x[n-1],y[m-1]) 
      alignment <- cbind(curr_align_col,alignment)
      n=n-1; m=m-1; # move back diagonally
    }else if(Tmat[n,m]==2){ 
      curr_align_col <- rbind(x[n-1],"-") # put gap in top seq
      alignment <- cbind(curr_align_col,alignment)
      n=n-1 # move up
    }else{
      curr_align_col <- rbind("-",y[m-1]) # put gap in side seq
      alignment <- cbind(curr_align_col,alignment)
      m=m-1 # move left
    }           
  } # end while
  return(alignment)
} # end function

alignment2 <- show.alignment(x_str2,y_str2,align.list2$Tmat)
alignment2
write.table(alignment2,row.names=F,col.names=F,quote=F)

## Input 3
x_str3 <- "ATCGT"  # side sequence
y_str3 <- "TGGTG" # top sequence
match_score <- 1
mismatch_score <- -2
gap_penalty <- -1 

align.list3 <- make.alignment.matrices(x_str3, y_str3, match_score,
                                       mismatch_score, gap_penalty)
align.list3$Fmat
Fmat3 <- align.list3$Fmat
align.list3$Tmat
align.list3$score_out 

col = c("black","blue","red","yellow","green")
breaks = seq(min(Fmat3),max(Fmat3),len=length(col)+1)

heatmap.2(Fmat3[-1,-1], dendrogram='none', density.info="none", 
          Rowv=FALSE, Colv=FALSE, trace='none', 
          breaks = breaks, col = col,
          sepwidth=c(0.01,0.01),
          sepcolor="black",
          colsep=1:ncol(Fmat3),
          rowsep=1:nrow(Fmat3))

alignment3 <- show.alignment(x_str3,y_str3,align.list3$Tmat)
alignment3
write.table(alignment3,row.names=F,col.names=F,quote=F)