view baseline/Baseline_Main.r @ 97:fbc6307dd83b draft

planemo upload commit 91fd26f19cb20fc2ee4d87d7edbc9ce6099752b0

#########################################################################################
# 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)
#########################################################################################

op <- options();
options(showWarnCalls=FALSE, showErrorCalls=FALSE, warn=-1)
library('seqinr')
if( F & Sys.info()[1]=="Linux"){
  library("multicore")
}

# Load functions and initialize global variables
source("Baseline_Functions.r")

# Initialize parameters with user provided arguments
  arg <- commandArgs(TRUE)                       
  #arg = c(2,1,5,5,0,1,"1:26:38:55:65:104:116", "test.fasta","","sample")
  #arg = c(1,1,5,5,0,1,"1:38:55:65:104:116:200", "test.fasta","","sample")
  #arg = c(1,1,5,5,1,1,"1:26:38:55:65:104:116", "/home/mu37/Wu/Wu_Cloned_gapped_sequences_D-masked.fasta","/home/mu37/Wu/","Wu")
  testID <- as.numeric(arg[1])                    # 1 = Focused, 2 = Local
  species <- as.numeric(arg[2])                   # 1 = Human. 2 = Mouse
  substitutionModel <- as.numeric(arg[3])         # 0 = Uniform substitution, 1 = Smith DS et al. 1996, 5 = FiveS
  mutabilityModel <- as.numeric(arg[4])           # 0 = Uniform mutablity, 1 = Tri-nucleotide (Shapiro GS et al. 2002)  , 5 = FiveS
  clonal <- as.numeric(arg[5])                    # 0 = Independent sequences, 1 = Clonally related, 2 = Clonally related & only non-terminal mutations
  fixIndels <- as.numeric(arg[6])                 # 0 = Do nothing, 1 = Try and fix Indels
  region <- as.numeric(strsplit(arg[7],":")[[1]]) # StartPos:LastNucleotideF1:C1:F2:C2:F3:C3
  inputFilePath <- arg[8]                         # Full path to input file
  outputPath <- arg[9]                            # Full path to location of output files
  outputID <- arg[10]                             # ID for session output  
  

  if(testID==5){
    traitChangeModel <- 1
    if( !is.na(any(arg[11])) ) traitChangeModel <- as.numeric(arg[11])    # 1 <- Chothia 1998
    initializeTraitChange(traitChangeModel)    
  }
  
# Initialize other parameters/variables
    
  # Initialzie the codon table ( definitions of R/S )
  computeCodonTable(testID) 

  # Initialize   
  # Test Name
  testName<-"Focused"
  if(testID==2) testName<-"Local"
  if(testID==3) testName<-"Imbalanced"    
  if(testID==4) testName<-"ImbalancedSilent"    
    
  # Indel placeholders initialization
  indelPos <- NULL
  delPos <- NULL
  insPos <- NULL

  # Initialize in Tranistion & Mutability matrixes
  substitution <- initializeSubstitutionMatrix(substitutionModel,species)
  mutability <- initializeMutabilityMatrix(mutabilityModel,species)
  
  # FWR/CDR boundaries
  flagTrim <- F
  if( is.na(region[7])){
    flagTrim <- T
    region[7]<-region[6]
  }
  readStart = min(region,na.rm=T)
  readEnd = max(region,na.rm=T)
  if(readStart>1){
    region = region - (readStart - 1)
  }
  region_Nuc = c( (region[1]*3-2) , (region[2:7]*3) )
  region_Cod = region
  
  readStart = (readStart*3)-2
  readEnd = (readEnd*3)
    
    FWR_Nuc <- c( rep(TRUE,(region_Nuc[2])),
                  rep(FALSE,(region_Nuc[3]-region_Nuc[2])),
                  rep(TRUE,(region_Nuc[4]-region_Nuc[3])),
                  rep(FALSE,(region_Nuc[5]-region_Nuc[4])),
                  rep(TRUE,(region_Nuc[6]-region_Nuc[5])),
                  rep(FALSE,(region_Nuc[7]-region_Nuc[6]))
                )
    CDR_Nuc <- (1-FWR_Nuc)
    CDR_Nuc <- as.logical(CDR_Nuc)
    FWR_Nuc_Mat <- matrix( rep(FWR_Nuc,4), ncol=length(FWR_Nuc), nrow=4, byrow=T)
    CDR_Nuc_Mat <- matrix( rep(CDR_Nuc,4), ncol=length(CDR_Nuc), nrow=4, byrow=T)
    
    FWR_Codon <- c( rep(TRUE,(region[2])),
                  rep(FALSE,(region[3]-region[2])),
                  rep(TRUE,(region[4]-region[3])),
                  rep(FALSE,(region[5]-region[4])),
                  rep(TRUE,(region[6]-region[5])),
                  rep(FALSE,(region[7]-region[6]))
                )
    CDR_Codon <- (1-FWR_Codon)
    CDR_Codon <- as.logical(CDR_Codon)


# Read input FASTA file
  tryCatch(
    inputFASTA <- baseline.read.fasta(inputFilePath, seqtype="DNA",as.string=T,set.attributes=F,forceDNAtolower=F)
    , error = function(ex){
      cat("Error|Error reading input. Please enter or upload a valid FASTA file.\n")
      q()
    }
  )
  
  if (length(inputFASTA)==1) {
    cat("Error|Error reading input. Please enter or upload a valid FASTA file.\n")
    q()
  }

  # Process sequence IDs/names
  names(inputFASTA) <- sapply(names(inputFASTA),function(x){trim(x)})
  
  # Convert non nucleotide characters to N
  inputFASTA[length(inputFASTA)] = gsub("\t","",inputFASTA[length(inputFASTA)])
  inputFASTA <- lapply(inputFASTA,replaceNonFASTAChars)

  # Process the FASTA file and conver to Matrix[inputSequence, germlineSequence]
  processedInput <- processInputAdvanced(inputFASTA)
  matInput <- processedInput[[1]]
  germlines <- processedInput[[2]]
  lenGermlines = length(unique(germlines))
  groups <- processedInput[[3]]
  lenGroups = length(unique(groups))
  rm(processedInput)
  rm(inputFASTA)

#   # remove clones with less than 2 seqeunces
#   tableGL <- table(germlines)
#   singletons <- which(tableGL<8)
#   rowsToRemove <- match(singletons,germlines)
#   if(any(rowsToRemove)){    
#     matInput <- matInput[-rowsToRemove,]
#     germlines <- germlines[-rowsToRemove]    
#     groups <- groups[-rowsToRemove]
#   }
# 
#   # remove unproductive seqs
#   nonFuctionalSeqs <- sapply(rownames(matInput),function(x){any(grep("unproductive",x))})
#   if(any(nonFuctionalSeqs)){
#     if(sum(nonFuctionalSeqs)==length(germlines)){
#       write.table("Unproductive",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
#       q()      
#     }
#     matInput <- matInput[-which(nonFuctionalSeqs),]
#     germlines <- germlines[-which(nonFuctionalSeqs)]
#     germlines[1:length(germlines)] <- 1:length(germlines)
#     groups <- groups[-which(nonFuctionalSeqs)]
#   }
# 
#   if(class(matInput)=="character"){
#     write.table("All unproductive seqs",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
#     q()    
#   }
#   
#   if(nrow(matInput)<10 | is.null(nrow(matInput))){
#     write.table(paste(nrow(matInput), "seqs only",sep=""),file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
#     q()
#   }

# replace leading & trailing "-" with "N:
  matInput <- t(apply(matInput,1,replaceLeadingTrailingDashes,readEnd))
    
  # Trim (nucleotide) input sequences to the last codon
  #matInput[,1] <- apply(matrix(matInput[,1]),1,trimToLastCodon) 

#   # Check for Indels
#   if(fixIndels){
#     delPos <- fixDeletions(matInput)
#     insPos <- fixInsertions(matInput)
#   }else{
#     # Check for indels
#     indelPos <- checkForInDels(matInput)
#     indelPos <- apply(cbind(indelPos[[1]],indelPos[[2]]),1,function(x){(x[1]==T & x[2]==T)})
#   }
  
  # If indels are present, remove mutations in the seqeunce & throw warning at end
  #matInput[indelPos,] <- apply(matrix(matInput[indelPos,],nrow=sum(indelPos),ncol=2),1,function(x){x[1]=x[2]; return(x) })
  
  colnames(matInput)=c("Input","Germline")

  # If seqeunces are clonal, create effective sequence for each clone & modify germline/group definitions
  germlinesOriginal = NULL
  if(clonal){
    germlinesOriginal <- germlines
    collapseCloneResults <- tapply(1:nrow(matInput),germlines,function(i){
                                                                collapseClone(matInput[i,1],matInput[i[1],2],readEnd,nonTerminalOnly=(clonal-1))
                                                              })
    matInput = t(sapply(collapseCloneResults,function(x){return(x[[1]])}))
    names_groups = tapply(groups,germlines,function(x){names(x[1])})  
    groups = tapply(groups,germlines,function(x){array(x[1],dimnames=names(x[1]))})  
    names(groups) = names_groups
  
    names_germlines =  tapply(germlines,germlines,function(x){names(x[1])})  
    germlines = tapply(   germlines,germlines,function(x){array(x[1],dimnames=names(x[1]))}   )
    names(germlines) = names_germlines
    matInputErrors = sapply(collapseCloneResults,function(x){return(x[[2]])})  
  }


# Selection Analysis

  
#  if (length(germlines)>sequenceLimit) {
#    # Code to parallelize processing goes here
#    stop( paste("Error: Cannot process more than ", Upper_limit," sequences",sep="") )
#  }

#  if (length(germlines)<sequenceLimit) {}
  
    # Compute expected mutation frequencies
    matExpected <- getExpectedIndividual(matInput)
    
    # Count observed number of mutations in the different regions
    mutations <- lapply( 1:nrow(matInput),  function(i){
                                              #cat(i,"\n")
                                              seqI = s2c(matInput[i,1])
                                              seqG = s2c(matInput[i,2])
                                              matIGL = matrix(c(seqI,seqG),ncol=length(seqI),nrow=2,byrow=T)    
                                              retVal <- NA
                                              tryCatch(
                                                retVal <- analyzeMutations2NucUri(matIGL)
                                                , error = function(ex){
                                                  retVal <- NA
                                                }
                                              )                                              
                                              
                                              
                                              return( retVal )
                                            })

    matObserved <- t(sapply( mutations, processNucMutations2 ))
    numberOfSeqsWithMutations <- numberOfSeqsWithMutations(matObserved, testID)

    #if(sum(numberOfSeqsWithMutations)==0){
    #  write.table("No mutated sequences",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
    #  q()      
    #}
    
    matMutationInfo <- cbind(matObserved,matExpected)
    rm(matObserved,matExpected)
    
     
    #Bayesian  PDFs
    bayes_pdf = computeBayesianScore(matMutationInfo, test=testName, max_sigma=20,length_sigma=4001)
    bayesPDF_cdr = bayes_pdf[[1]]
    bayesPDF_fwr = bayes_pdf[[2]]    
    rm(bayes_pdf)

    bayesPDF_germlines_cdr = tapply(bayesPDF_cdr,germlines,function(x) groupPosteriors(x,length_sigma=4001))
    bayesPDF_germlines_fwr = tapply(bayesPDF_fwr,germlines,function(x) groupPosteriors(x,length_sigma=4001))
    
    bayesPDF_groups_cdr = tapply(bayesPDF_cdr,groups,function(x) groupPosteriors(x,length_sigma=4001))
    bayesPDF_groups_fwr = tapply(bayesPDF_fwr,groups,function(x) groupPosteriors(x,length_sigma=4001))
    
    if(lenGroups>1){
      groups <- c(groups,lenGroups+1)
      names(groups)[length(groups)] = "All sequences combined"
      bayesPDF_groups_cdr[[lenGroups+1]] =   groupPosteriors(bayesPDF_groups_cdr,length_sigma=4001)
      bayesPDF_groups_fwr[[lenGroups+1]] =   groupPosteriors(bayesPDF_groups_fwr,length_sigma=4001)
    }
    
    #Bayesian  Outputs
    bayes_cdr =  t(sapply(bayesPDF_cdr,calcBayesOutputInfo))
    bayes_fwr =  t(sapply(bayesPDF_fwr,calcBayesOutputInfo))
    bayes_germlines_cdr =  t(sapply(bayesPDF_germlines_cdr,calcBayesOutputInfo))
    bayes_germlines_fwr =  t(sapply(bayesPDF_germlines_fwr,calcBayesOutputInfo))
    bayes_groups_cdr =  t(sapply(bayesPDF_groups_cdr,calcBayesOutputInfo))
    bayes_groups_fwr =  t(sapply(bayesPDF_groups_fwr,calcBayesOutputInfo))
    
    #P-values
    simgaP_cdr = sapply(bayesPDF_cdr,computeSigmaP)
    simgaP_fwr = sapply(bayesPDF_fwr,computeSigmaP)
    
    simgaP_germlines_cdr = sapply(bayesPDF_germlines_cdr,computeSigmaP)
    simgaP_germlines_fwr = sapply(bayesPDF_germlines_fwr,computeSigmaP)
    
    simgaP_groups_cdr = sapply(bayesPDF_groups_cdr,computeSigmaP)
    simgaP_groups_fwr = sapply(bayesPDF_groups_fwr,computeSigmaP)
    
    
    #Format output
    
    # Round expected mutation frequencies to 3 decimal places
    matMutationInfo[germlinesOriginal[indelPos],] = NA
    if(nrow(matMutationInfo)==1){
      matMutationInfo[5:8] = round(matMutationInfo[,5:8]/sum(matMutationInfo[,5:8],na.rm=T),3)
    }else{
      matMutationInfo[,5:8] = t(round(apply(matMutationInfo[,5:8],1,function(x){ return(x/sum(x,na.rm=T)) }),3))
    }
    
    listPDFs = list()
    nRows = length(unique(groups)) + length(unique(germlines)) + length(groups)
    
    matOutput = matrix(NA,ncol=18,nrow=nRows)
    rowNumb = 1
    for(G in unique(groups)){
      #print(G)
      matOutput[rowNumb,c(1,2,11:18)] = c("Group",names(groups)[groups==G][1],bayes_groups_cdr[G,],bayes_groups_fwr[G,],simgaP_groups_cdr[G],simgaP_groups_fwr[G])
      listPDFs[[rowNumb]] = list("CDR"=bayesPDF_groups_cdr[[G]],"FWR"=bayesPDF_groups_fwr[[G]])
      names(listPDFs)[rowNumb] = names(groups[groups==paste(G)])[1]
      #if(names(groups)[which(groups==G)[1]]!="All sequences combined"){
      gs = unique(germlines[groups==G])
      rowNumb = rowNumb+1
      if( !is.na(gs) ){
        for( g in gs ){
          matOutput[rowNumb,c(1,2,11:18)] = c("Germline",names(germlines)[germlines==g][1],bayes_germlines_cdr[g,],bayes_germlines_fwr[g,],simgaP_germlines_cdr[g],simgaP_germlines_fwr[g])
          listPDFs[[rowNumb]] = list("CDR"=bayesPDF_germlines_cdr[[g]],"FWR"=bayesPDF_germlines_fwr[[g]])
          names(listPDFs)[rowNumb] = names(germlines[germlines==paste(g)])[1]
          rowNumb = rowNumb+1
          indexesOfInterest = which(germlines==g)
          numbSeqsOfInterest =  length(indexesOfInterest)
          rowNumb = seq(rowNumb,rowNumb+(numbSeqsOfInterest-1))
          matOutput[rowNumb,] = matrix(   c(  rep("Sequence",numbSeqsOfInterest),
                                              rownames(matInput)[indexesOfInterest],
                                              c(matMutationInfo[indexesOfInterest,1:4]),
                                              c(matMutationInfo[indexesOfInterest,5:8]),
                                              c(bayes_cdr[indexesOfInterest,]),
                                              c(bayes_fwr[indexesOfInterest,]),
                                              c(simgaP_cdr[indexesOfInterest]),
                                              c(simgaP_fwr[indexesOfInterest])                                              
          ), ncol=18, nrow=numbSeqsOfInterest,byrow=F)
          increment=0
          for( ioi in indexesOfInterest){
            listPDFs[[min(rowNumb)+increment]] =  list("CDR"=bayesPDF_cdr[[ioi]] , "FWR"=bayesPDF_fwr[[ioi]])
            names(listPDFs)[min(rowNumb)+increment] = rownames(matInput)[ioi]
            increment = increment + 1
          }
          rowNumb=max(rowNumb)+1

        }
      }
    }
    colsToFormat = 11:18
    matOutput[,colsToFormat] = formatC(  matrix(as.numeric(matOutput[,colsToFormat]), nrow=nrow(matOutput), ncol=length(colsToFormat)) ,  digits=3)
    matOutput[matOutput== " NaN"] = NA
    
    
    
    colnames(matOutput) = c("Type", "ID", "Observed_CDR_R", "Observed_CDR_S", "Observed_FWR_R", "Observed_FWR_S",
                            "Expected_CDR_R", "Expected_CDR_S", "Expected_FWR_R", "Expected_FWR_S",
                            paste( rep(testName,6), rep(c("Sigma","CIlower","CIupper"),2),rep(c("CDR","FWR"),each=3), sep="_"),
                            paste( rep(testName,2), rep("P",2),c("CDR","FWR"), sep="_")
    )
    fileName = paste(outputPath,outputID,".txt",sep="")
    write.table(matOutput,file=fileName,quote=F,sep="\t",row.names=T,col.names=NA)
    fileName = paste(outputPath,outputID,".RData",sep="")
    save(listPDFs,file=fileName)

indelWarning = FALSE
if(sum(indelPos)>0){
  indelWarning = "<P>Warning: The following sequences have either gaps and/or deletions, and have been ommited from the analysis.";
  indelWarning = paste( indelWarning , "<UL>", sep="" )
  for(indels in names(indelPos)[indelPos]){
    indelWarning = paste( indelWarning , "<LI>", indels, "</LI>", sep="" )
  }
  indelWarning = paste( indelWarning , "</UL></P>", sep="" )
}

cloneWarning = FALSE
if(clonal==1){
  if(sum(matInputErrors)>0){
    cloneWarning = "<P>Warning: The following clones have sequences of unequal length.";
    cloneWarning = paste( cloneWarning , "<UL>", sep="" )
    for(clone in names(matInputErrors)[matInputErrors]){
      cloneWarning = paste( cloneWarning , "<LI>", names(germlines)[as.numeric(clone)], "</LI>", sep="" )
    }
    cloneWarning = paste( cloneWarning , "</UL></P>", sep="" )
  }
}
cat(paste("Success",outputID,indelWarning,cloneWarning,sep="|"))