.packageName <- "aaMI"
aaMI <- function(file){
  AA <- c(LETTERS[-c(2,10,15,21,24,26)],"-")
  AA.F <- matrix(0,nrow=21,ncol=ncol(file),dimnames=list(AA,seq(1:ncol(file))))
  H <- vector(length=ncol(file))
  MIt <- matrix(0,nrow=ncol(file),ncol=ncol(file))

  for(i in 1:ncol(file))
    AA.F[names(summary(file[[i]])),i] <- (summary(file[[i]])/length(file[[i]]))

  for(i in 1:ncol(file))
   H[i] <- -1*sum(AA.F[AA.F[,i] != 0,i]*log2(AA.F[AA.F[,i] != 0,i]))

  for(i in 1:ncol(file)){
    for(j in i:ncol(file)){
      FrqT <- table(file[,i],file[,j])/sum(table(file[,i],file[,j]))
      MIt[i,j] <- H[i] + H[j] + sum(FrqT[FrqT != 0] * log2(FrqT[FrqT != 0]))
    }
  }
  return(MIt)
}
aaMIn <- function(file){
  AA <- c(LETTERS[-c(2,10,15,21,24,26)],"-")
  AA.F <- matrix(0,nrow=21,ncol=ncol(file),dimnames=list(AA,seq(1:ncol(file))))
  H <- vector(length=ncol(file))
  MIt <- matrix(0,nrow=ncol(file),ncol=ncol(file))

  for(i in 1:ncol(file))
    AA.F[names(summary(file[[i]])),i] <- (summary(file[[i]])/length(file[[i]]))

  for(i in 1:ncol(file))
   H[i] <- -1*sum(AA.F[AA.F[,i] != 0,i]*log2(AA.F[AA.F[,i] != 0,i]))

  for(i in 1:ncol(file)){
    for(j in i:ncol(file)){
      FrqT <- table(file[,i],file[,j])/sum(table(file[,i],file[,j]))
     JointH <- sum(FrqT[FrqT != 0] * log2(FrqT[FrqT != 0]))
      MIt[i,j] <- (H[i]+H[j] + JointH)/(-1*JointH)
    }
  }
  return(MIt)
}
read.CX <- function(file){
  CX1 <- read.table(file,header=F,fill=T)
  CX1 <- CX1[2:nrow(CX1),1:2] #remove "ClustalX (n) multiple sequence alignment"
  CX1[,1] <- as.character(CX1[,1])
  CX1[,2] <- as.character(CX1[,2])
  SeqID <- unique(grep("[[:alnum:]]",CX1[,1],value=T))
  CX2 <- matrix("",nrow=length(SeqID),ncol=2)
  CX2[,1] <- SeqID
  for(i in 1:length(SeqID)){
    rep <- which(CX1[,1] == CX2[i,1])
    for(j in 1:length(rep))
      CX2[i,2] <- paste(CX2[i,2],CX1[rep[j],2],sep="")
  }
  CX2[,2] <- toupper(CX2[,2]) # Make sure all aa symbols are upper case.
  align1 <- matrix(NA,nrow=nrow(CX2),ncol=nchar(CX2[1,2]),
    dimnames=list(CX2[,1],seq(1:nchar(CX2[1,2]))))
  for(i in 1:nrow(align1))
    align1[i,] <- unlist(strsplit(CX2[i,2],split=""))
  align1 <- as.data.frame(align1)
  return(align1)
}
read.FASTA <- function(file){
  FA1 <- read.table(file,header=F,fill=T,as.is=T)
  SeqID <- grep(">",FA1[[1]]) # extract row numbers of sequence names
  skip <- nrow(FA1)/length(SeqID)
  FA2 <- matrix("",nrow=length(SeqID),ncol=2)
  for(i in 0:nrow(FA1)-1){
    F2row <- trunc(i/skip)+1
    if((i+1) %in% SeqID) {FA2[F2row,1] <- FA1[i+1,]
     }else{
       FA2[F2row,2] <- paste(FA2[F2row,2],FA1[i+1,],sep="")
     } # combine multiple rows of sequence into one row 
  }
  FA2[,2] <- toupper(FA2[,2]) # Make sure all aa symbols are upper case.
  align1 <- matrix(NA,nrow=nrow(FA2),ncol=nchar(as.character(FA2[1,2])),
                   dimnames=list(FA2[,1],seq(1:nchar(as.character(FA2[1,2])))))
  for(i in 1:nrow(align1))
    align1[i,] <- unlist(strsplit(as.character(FA2[i,2]),split=""))
  align1 <- as.data.frame(align1)
  row.names(align1) <- gsub(">","",row.names(align1))
  return(align1)
}
read.GDoc <- function(file){
  GD1 <- read.table(file,header=F,fill=T)
  GD1 <- GD1[(which(GD1[,1] == "1")):nrow(GD1),1:6] # remove initial gibberish
  for(i in 1:ncol(GD1))
    GD1[,i] <- as.character(GD1[,i])
  SeqID <- unique(grep("[[:alpha:]]",GD1[,1],value=T)) # extract sequence names
  GD2 <- matrix("",nrow=length(SeqID),ncol=2)
  GD2[,1] <- SeqID
  for(i in 1:length(SeqID)){
    rep <- which(GD1[,1] == GD2[i,1])
    for(j in 1:length(rep)){
      for(k in 2:ncol(GD1))
        GD2[i,2] <- paste(GD2[i,2],GD1[rep[j],k],sep="")
    }
  }
  GD2[,2] <- toupper(GD2[,2]) # Make sure all aa symbols are upper case.
  GD2[,2] <- gsub("\\.","-",GD2[,2]) # change style of gap character
  align1 <- matrix(NA,nrow=nrow(GD2),ncol=nchar(GD2[1,2]),
    dimnames=list(GD2[,1],seq(1:nchar(GD2[1,2]))))
  for(i in 1:nrow(align1))
    align1[i,] <- unlist(strsplit(GD2[i,2],split=""))
  align1 <- as.data.frame(align1)
  return(align1)
}
