# HG changeset patch # User davidvanzessen # Date 1405428229 14400 # Node ID d685e7ba0ed4673f31d577158c06ffbc0b8d3217 Uploaded diff -r 000000000000 -r d685e7ba0ed4 Baseline_Functions.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Baseline_Functions.r Tue Jul 15 08:43:49 2014 -0400 @@ -0,0 +1,2286 @@ +######################################################################################### +# 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 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({lx1){ + while( i1 & 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")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( i0)) ){ + +# 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) ) +} + diff -r 000000000000 -r d685e7ba0ed4 Baseline_Main.r --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Baseline_Main.r Tue Jul 15 08:43:49 2014 -0400 @@ -0,0 +1,389 @@ +######################################################################################### +# 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( 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)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 = "

Warning: The following sequences have either gaps and/or deletions, and have been ommited from the analysis."; + indelWarning = paste( indelWarning , "

    ", sep="" ) + for(indels in names(indelPos)[indelPos]){ + indelWarning = paste( indelWarning , "
  • ", indels, "
    • ", sep="" ) + } + indelWarning = paste( indelWarning , "

", sep="" ) +} + +cloneWarning = FALSE +if(clonal==1){ + if(sum(matInputErrors)>0){ + cloneWarning = "

Warning: The following clones have sequences of unequal length."; + cloneWarning = paste( cloneWarning , "

    ", sep="" ) + for(clone in names(matInputErrors)[matInputErrors]){ + cloneWarning = paste( cloneWarning , "
  • ", names(germlines)[as.numeric(clone)], "
    • ", sep="" ) + } + cloneWarning = paste( cloneWarning , "

", sep="" ) + } +} +cat(paste("Success",outputID,indelWarning,cloneWarning,sep="|")) diff -r 000000000000 -r d685e7ba0ed4 FiveS_Mutability.RData Binary file FiveS_Mutability.RData has changed diff -r 000000000000 -r d685e7ba0ed4 FiveS_Substitution.RData Binary file FiveS_Substitution.RData has changed diff -r 000000000000 -r d685e7ba0ed4 baseline.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/baseline.xml Tue Jul 15 08:43:49 2014 -0400 @@ -0,0 +1,50 @@ + + Bayesian Estimation of Antigen-Driven Selection + + wrapper.sh $ss $species $substitution $mutability $clonal $fixindels "$boundaries" $in_file $out_file.files_path result $out_file $out_file2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Gur Yaari; Mohamed Uduman; Steven H. Kleinstein. Quantifying selection in high-throughput Immunoglobulin sequencing data sets. Nucleic Acids Res. 2012 May 27. + + Mohamed Uduman; Gur Yaari; Uri Hershberg; Mark J. Shlomchik; Steven H. Kleinstein. Detecting selection in immunoglobulin sequences. Nucleic Acids Res. 2011 Jul;39(Web Server issue):W499-504. + + + diff -r 000000000000 -r d685e7ba0ed4 wrapper.sh --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/wrapper.sh Tue Jul 15 08:43:49 2014 -0400 @@ -0,0 +1,33 @@ +#!/bin/bash +dir="$(cd "$(dirname "$0")" && pwd)" + +testID=$1 +species=$2 +substitutionModel=$3 +mutabilityModel=$4 +clonal=$5 +fixIndels=$6 +region=$7 +input=$8 +outDir=$9 +#mkdir $outDir +outID=${10} +outFile=${11} +outFile2=${12} + +cd $dir +Rscript --verbose $dir/Baseline_Main.r $testID $species $substitutionModel $mutabilityModel $clonal $fixIndels $region $input $dir/ $outID 2> /dev/null + +mv $dir/${outID}.txt $outFile +mv $dir/${outID}.RData $outFile2 + +#tail -n+2 $outDir/${outID}.txt > $outDir/tmp.txt +#echo "" >> $outFile +#echo "" >> $outFile +#while read 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 Focused_Sigma_CDR Focused_CIlower_CDR Focused_CIupper_CDR Focused_Sigma_FWR Focused_CIlower_FWR Focused_CIupper_FWR Focused_P_CDR Focused_P_FWR +#do +# echo "" >> $outFile +#done < $outDir/tmp.txt +#echo "
TypeIDObserved_CDR_RObserved_CDR_SObserved_FWR_RObserved_FWR_SExpected_CDR_RExpected_CDR_SExpected_FWR_RExpected_FWR_SFocused_Sigma_CDRFocused_CIlower_CDRFocused_CIupper_CDRFocused_Sigma_FWRFocused_CIlower_FWRFocused_CIupper_FWRFocused_P_CDRFocused_P_FWR
$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$Focused_Sigma_CDR$Focused_CIlower_CDR$Focused_CIupper_CDR$Focused_Sigma_FWR$Focused_CIlower_FWR$Focused_CIupper_FWR$Focused_P_CDR$Focused_P_FWR
" >> $outFile +#rm $outDir/tmp.txt +