Mercurial > repos > davidvanzessen > mutation_analysis

comparison baseline/Baseline_Main.r @ 0:8a5a2abbb870 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 29 Aug 2016 05:36:10 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8a5a2abbb870
1 #########################################################################################
2 # License Agreement
3 #
4 # THIS WORK IS PROVIDED UNDER THE TERMS OF THIS CREATIVE COMMONS PUBLIC LICENSE
5 # ("CCPL" OR "LICENSE"). THE WORK IS PROTECTED BY COPYRIGHT AND/OR OTHER
6 # APPLICABLE LAW. ANY USE OF THE WORK OTHER THAN AS AUTHORIZED UNDER THIS LICENSE
7 # OR COPYRIGHT LAW IS PROHIBITED.
8 #
9 # BY EXERCISING ANY RIGHTS TO THE WORK PROVIDED HERE, YOU ACCEPT AND AGREE TO BE
10 # BOUND BY THE TERMS OF THIS LICENSE. TO THE EXTENT THIS LICENSE MAY BE CONSIDERED
11 # TO BE A CONTRACT, THE LICENSOR GRANTS YOU THE RIGHTS CONTAINED HERE IN
12 # CONSIDERATION OF YOUR ACCEPTANCE OF SUCH TERMS AND CONDITIONS.
13 #
14 # BASELIne: Bayesian Estimation of Antigen-Driven Selection in Immunoglobulin Sequences
15 # Coded by: Mohamed Uduman & Gur Yaari
16 # Copyright 2012 Kleinstein Lab
17 # Version: 1.3 (01/23/2014)
18 #########################################################################################
19
20 op <- options();
21 options(showWarnCalls=FALSE, showErrorCalls=FALSE, warn=-1)
22 library('seqinr')
23 if( F & Sys.info()[1]=="Linux"){
24 library("multicore")
25 }
26
27 # Load functions and initialize global variables
28 source("Baseline_Functions.r")
29
30 # Initialize parameters with user provided arguments
31 arg <- commandArgs(TRUE)
32 #arg = c(2,1,5,5,0,1,"1:26:38:55:65:104:116", "test.fasta","","sample")
33 #arg = c(1,1,5,5,0,1,"1:38:55:65:104:116:200", "test.fasta","","sample")
34 #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")
35 testID <- as.numeric(arg[1]) # 1 = Focused, 2 = Local
36 species <- as.numeric(arg[2]) # 1 = Human. 2 = Mouse
37 substitutionModel <- as.numeric(arg[3]) # 0 = Uniform substitution, 1 = Smith DS et al. 1996, 5 = FiveS
38 mutabilityModel <- as.numeric(arg[4]) # 0 = Uniform mutablity, 1 = Tri-nucleotide (Shapiro GS et al. 2002) , 5 = FiveS
39 clonal <- as.numeric(arg[5]) # 0 = Independent sequences, 1 = Clonally related, 2 = Clonally related & only non-terminal mutations
40 fixIndels <- as.numeric(arg[6]) # 0 = Do nothing, 1 = Try and fix Indels
41 region <- as.numeric(strsplit(arg[7],":")[[1]]) # StartPos:LastNucleotideF1:C1:F2:C2:F3:C3
42 inputFilePath <- arg[8] # Full path to input file
43 outputPath <- arg[9] # Full path to location of output files
44 outputID <- arg[10] # ID for session output
45
46
47 if(testID==5){
48 traitChangeModel <- 1
49 if( !is.na(any(arg[11])) ) traitChangeModel <- as.numeric(arg[11]) # 1 <- Chothia 1998
50 initializeTraitChange(traitChangeModel)
51 }
52
53 # Initialize other parameters/variables
54
55 # Initialzie the codon table ( definitions of R/S )
56 computeCodonTable(testID)
57
58 # Initialize
59 # Test Name
60 testName<-"Focused"
61 if(testID==2) testName<-"Local"
62 if(testID==3) testName<-"Imbalanced"
63 if(testID==4) testName<-"ImbalancedSilent"
64
65 # Indel placeholders initialization
66 indelPos <- NULL
67 delPos <- NULL
68 insPos <- NULL
69
70 # Initialize in Tranistion & Mutability matrixes
71 substitution <- initializeSubstitutionMatrix(substitutionModel,species)
72 mutability <- initializeMutabilityMatrix(mutabilityModel,species)
73
74 # FWR/CDR boundaries
75 flagTrim <- F
76 if( is.na(region[7])){
77 flagTrim <- T
78 region[7]<-region[6]
79 }
80 readStart = min(region,na.rm=T)
81 readEnd = max(region,na.rm=T)
82 if(readStart>1){
83 region = region - (readStart - 1)
84 }
85 region_Nuc = c( (region[1]*3-2) , (region[2:7]*3) )
86 region_Cod = region
87
88 readStart = (readStart*3)-2
89 readEnd = (readEnd*3)
90
91 FWR_Nuc <- c( rep(TRUE,(region_Nuc[2])),
92 rep(FALSE,(region_Nuc[3]-region_Nuc[2])),
93 rep(TRUE,(region_Nuc[4]-region_Nuc[3])),
94 rep(FALSE,(region_Nuc[5]-region_Nuc[4])),
95 rep(TRUE,(region_Nuc[6]-region_Nuc[5])),
96 rep(FALSE,(region_Nuc[7]-region_Nuc[6]))
97 )
98 CDR_Nuc <- (1-FWR_Nuc)
99 CDR_Nuc <- as.logical(CDR_Nuc)
100 FWR_Nuc_Mat <- matrix( rep(FWR_Nuc,4), ncol=length(FWR_Nuc), nrow=4, byrow=T)
101 CDR_Nuc_Mat <- matrix( rep(CDR_Nuc,4), ncol=length(CDR_Nuc), nrow=4, byrow=T)
102
103 FWR_Codon <- c( rep(TRUE,(region[2])),
104 rep(FALSE,(region[3]-region[2])),
105 rep(TRUE,(region[4]-region[3])),
106 rep(FALSE,(region[5]-region[4])),
107 rep(TRUE,(region[6]-region[5])),
108 rep(FALSE,(region[7]-region[6]))
109 )
110 CDR_Codon <- (1-FWR_Codon)
111 CDR_Codon <- as.logical(CDR_Codon)
112
113
114 # Read input FASTA file
115 tryCatch(
116 inputFASTA <- baseline.read.fasta(inputFilePath, seqtype="DNA",as.string=T,set.attributes=F,forceDNAtolower=F)
117 , error = function(ex){
118 cat("Error|Error reading input. Please enter or upload a valid FASTA file.\n")
119 q()
120 }
121 )
122
123 if (length(inputFASTA)==1) {
124 cat("Error|Error reading input. Please enter or upload a valid FASTA file.\n")
125 q()
126 }
127
128 # Process sequence IDs/names
129 names(inputFASTA) <- sapply(names(inputFASTA),function(x){trim(x)})
130
131 # Convert non nucleotide characters to N
132 inputFASTA[length(inputFASTA)] = gsub("\t","",inputFASTA[length(inputFASTA)])
133 inputFASTA <- lapply(inputFASTA,replaceNonFASTAChars)
134
135 # Process the FASTA file and conver to Matrix[inputSequence, germlineSequence]
136 processedInput <- processInputAdvanced(inputFASTA)
137 matInput <- processedInput[[1]]
138 germlines <- processedInput[[2]]
139 lenGermlines = length(unique(germlines))
140 groups <- processedInput[[3]]
141 lenGroups = length(unique(groups))
142 rm(processedInput)
143 rm(inputFASTA)
144
145 # # remove clones with less than 2 seqeunces
146 # tableGL <- table(germlines)
147 # singletons <- which(tableGL<8)
148 # rowsToRemove <- match(singletons,germlines)
149 # if(any(rowsToRemove)){
150 # matInput <- matInput[-rowsToRemove,]
151 # germlines <- germlines[-rowsToRemove]
152 # groups <- groups[-rowsToRemove]
153 # }
154 #
155 # # remove unproductive seqs
156 # nonFuctionalSeqs <- sapply(rownames(matInput),function(x){any(grep("unproductive",x))})
157 # if(any(nonFuctionalSeqs)){
158 # if(sum(nonFuctionalSeqs)==length(germlines)){
159 # write.table("Unproductive",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
160 # q()
161 # }
162 # matInput <- matInput[-which(nonFuctionalSeqs),]
163 # germlines <- germlines[-which(nonFuctionalSeqs)]
164 # germlines[1:length(germlines)] <- 1:length(germlines)
165 # groups <- groups[-which(nonFuctionalSeqs)]
166 # }
167 #
168 # if(class(matInput)=="character"){
169 # write.table("All unproductive seqs",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
170 # q()
171 # }
172 #
173 # if(nrow(matInput)<10 | is.null(nrow(matInput))){
174 # write.table(paste(nrow(matInput), "seqs only",sep=""),file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
175 # q()
176 # }
177
178 # replace leading & trailing "-" with "N:
179 matInput <- t(apply(matInput,1,replaceLeadingTrailingDashes,readEnd))
180
181 # Trim (nucleotide) input sequences to the last codon
182 #matInput[,1] <- apply(matrix(matInput[,1]),1,trimToLastCodon)
183
184 # # Check for Indels
185 # if(fixIndels){
186 # delPos <- fixDeletions(matInput)
187 # insPos <- fixInsertions(matInput)
188 # }else{
189 # # Check for indels
190 # indelPos <- checkForInDels(matInput)
191 # indelPos <- apply(cbind(indelPos[[1]],indelPos[[2]]),1,function(x){(x[1]==T & x[2]==T)})
192 # }
193
194 # If indels are present, remove mutations in the seqeunce & throw warning at end
195 #matInput[indelPos,] <- apply(matrix(matInput[indelPos,],nrow=sum(indelPos),ncol=2),1,function(x){x[1]=x[2]; return(x) })
196
197 colnames(matInput)=c("Input","Germline")
198
199 # If seqeunces are clonal, create effective sequence for each clone & modify germline/group definitions
200 germlinesOriginal = NULL
201 if(clonal){
202 germlinesOriginal <- germlines
203 collapseCloneResults <- tapply(1:nrow(matInput),germlines,function(i){
204 collapseClone(matInput[i,1],matInput[i[1],2],readEnd,nonTerminalOnly=(clonal-1))
205 })
206 matInput = t(sapply(collapseCloneResults,function(x){return(x[[1]])}))
207 names_groups = tapply(groups,germlines,function(x){names(x[1])})
208 groups = tapply(groups,germlines,function(x){array(x[1],dimnames=names(x[1]))})
209 names(groups) = names_groups
210
211 names_germlines = tapply(germlines,germlines,function(x){names(x[1])})
212 germlines = tapply( germlines,germlines,function(x){array(x[1],dimnames=names(x[1]))} )
213 names(germlines) = names_germlines
214 matInputErrors = sapply(collapseCloneResults,function(x){return(x[[2]])})
215 }
216
217
218 # Selection Analysis
219
220
221 # if (length(germlines)>sequenceLimit) {
222 # # Code to parallelize processing goes here
223 # stop( paste("Error: Cannot process more than ", Upper_limit," sequences",sep="") )
224 # }
225
226 # if (length(germlines)<sequenceLimit) {}
227
228 # Compute expected mutation frequencies
229 matExpected <- getExpectedIndividual(matInput)
230
231 # Count observed number of mutations in the different regions
232 mutations <- lapply( 1:nrow(matInput), function(i){
233 #cat(i,"\n")
234 seqI = s2c(matInput[i,1])
235 seqG = s2c(matInput[i,2])
236 matIGL = matrix(c(seqI,seqG),ncol=length(seqI),nrow=2,byrow=T)
237 retVal <- NA
238 tryCatch(
239 retVal <- analyzeMutations2NucUri(matIGL)
240 , error = function(ex){
241 retVal <- NA
242 }
243 )
244
245
246 return( retVal )
247 })
248
249 matObserved <- t(sapply( mutations, processNucMutations2 ))
250 numberOfSeqsWithMutations <- numberOfSeqsWithMutations(matObserved, testID)
251
252 #if(sum(numberOfSeqsWithMutations)==0){
253 # write.table("No mutated sequences",file=paste(outputPath,outputID,".txt",sep=""),quote=F,sep="\t",row.names=F,col.names=T)
254 # q()
255 #}
256
257 matMutationInfo <- cbind(matObserved,matExpected)
258 rm(matObserved,matExpected)
259
260
261 #Bayesian PDFs
262 bayes_pdf = computeBayesianScore(matMutationInfo, test=testName, max_sigma=20,length_sigma=4001)
263 bayesPDF_cdr = bayes_pdf[[1]]
264 bayesPDF_fwr = bayes_pdf[[2]]
265 rm(bayes_pdf)
266
267 bayesPDF_germlines_cdr = tapply(bayesPDF_cdr,germlines,function(x) groupPosteriors(x,length_sigma=4001))
268 bayesPDF_germlines_fwr = tapply(bayesPDF_fwr,germlines,function(x) groupPosteriors(x,length_sigma=4001))
269
270 bayesPDF_groups_cdr = tapply(bayesPDF_cdr,groups,function(x) groupPosteriors(x,length_sigma=4001))
271 bayesPDF_groups_fwr = tapply(bayesPDF_fwr,groups,function(x) groupPosteriors(x,length_sigma=4001))
272
273 if(lenGroups>1){
274 groups <- c(groups,lenGroups+1)
275 names(groups)[length(groups)] = "All sequences combined"
276 bayesPDF_groups_cdr[[lenGroups+1]] = groupPosteriors(bayesPDF_groups_cdr,length_sigma=4001)
277 bayesPDF_groups_fwr[[lenGroups+1]] = groupPosteriors(bayesPDF_groups_fwr,length_sigma=4001)
278 }
279
280 #Bayesian Outputs
281 bayes_cdr = t(sapply(bayesPDF_cdr,calcBayesOutputInfo))
282 bayes_fwr = t(sapply(bayesPDF_fwr,calcBayesOutputInfo))
283 bayes_germlines_cdr = t(sapply(bayesPDF_germlines_cdr,calcBayesOutputInfo))
284 bayes_germlines_fwr = t(sapply(bayesPDF_germlines_fwr,calcBayesOutputInfo))
285 bayes_groups_cdr = t(sapply(bayesPDF_groups_cdr,calcBayesOutputInfo))
286 bayes_groups_fwr = t(sapply(bayesPDF_groups_fwr,calcBayesOutputInfo))
287
288 #P-values
289 simgaP_cdr = sapply(bayesPDF_cdr,computeSigmaP)
290 simgaP_fwr = sapply(bayesPDF_fwr,computeSigmaP)
291
292 simgaP_germlines_cdr = sapply(bayesPDF_germlines_cdr,computeSigmaP)
293 simgaP_germlines_fwr = sapply(bayesPDF_germlines_fwr,computeSigmaP)
294
295 simgaP_groups_cdr = sapply(bayesPDF_groups_cdr,computeSigmaP)
296 simgaP_groups_fwr = sapply(bayesPDF_groups_fwr,computeSigmaP)
297
298
299 #Format output
300
301 # Round expected mutation frequencies to 3 decimal places
302 matMutationInfo[germlinesOriginal[indelPos],] = NA
303 if(nrow(matMutationInfo)==1){
304 matMutationInfo[5:8] = round(matMutationInfo[,5:8]/sum(matMutationInfo[,5:8],na.rm=T),3)
305 }else{
306 matMutationInfo[,5:8] = t(round(apply(matMutationInfo[,5:8],1,function(x){ return(x/sum(x,na.rm=T)) }),3))
307 }
308
309 listPDFs = list()
310 nRows = length(unique(groups)) + length(unique(germlines)) + length(groups)
311
312 matOutput = matrix(NA,ncol=18,nrow=nRows)
313 rowNumb = 1
314 for(G in unique(groups)){
315 #print(G)
316 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])
317 listPDFs[[rowNumb]] = list("CDR"=bayesPDF_groups_cdr[[G]],"FWR"=bayesPDF_groups_fwr[[G]])
318 names(listPDFs)[rowNumb] = names(groups[groups==paste(G)])[1]
319 #if(names(groups)[which(groups==G)[1]]!="All sequences combined"){
320 gs = unique(germlines[groups==G])
321 rowNumb = rowNumb+1
322 if( !is.na(gs) ){
323 for( g in gs ){
324 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])
325 listPDFs[[rowNumb]] = list("CDR"=bayesPDF_germlines_cdr[[g]],"FWR"=bayesPDF_germlines_fwr[[g]])
326 names(listPDFs)[rowNumb] = names(germlines[germlines==paste(g)])[1]
327 rowNumb = rowNumb+1
328 indexesOfInterest = which(germlines==g)
329 numbSeqsOfInterest = length(indexesOfInterest)
330 rowNumb = seq(rowNumb,rowNumb+(numbSeqsOfInterest-1))
331 matOutput[rowNumb,] = matrix( c( rep("Sequence",numbSeqsOfInterest),
332 rownames(matInput)[indexesOfInterest],
333 c(matMutationInfo[indexesOfInterest,1:4]),
334 c(matMutationInfo[indexesOfInterest,5:8]),
335 c(bayes_cdr[indexesOfInterest,]),
336 c(bayes_fwr[indexesOfInterest,]),
337 c(simgaP_cdr[indexesOfInterest]),
338 c(simgaP_fwr[indexesOfInterest])
339 ), ncol=18, nrow=numbSeqsOfInterest,byrow=F)
340 increment=0
341 for( ioi in indexesOfInterest){
342 listPDFs[[min(rowNumb)+increment]] = list("CDR"=bayesPDF_cdr[[ioi]] , "FWR"=bayesPDF_fwr[[ioi]])
343 names(listPDFs)[min(rowNumb)+increment] = rownames(matInput)[ioi]
344 increment = increment + 1
345 }
346 rowNumb=max(rowNumb)+1
347
348 }
349 }
350 }
351 colsToFormat = 11:18
352 matOutput[,colsToFormat] = formatC( matrix(as.numeric(matOutput[,colsToFormat]), nrow=nrow(matOutput), ncol=length(colsToFormat)) , digits=3)
353 matOutput[matOutput== " NaN"] = NA
354
355
356
357 colnames(matOutput) = c("Type", "ID", "Observed_CDR_R", "Observed_CDR_S", "Observed_FWR_R", "Observed_FWR_S",
358 "Expected_CDR_R", "Expected_CDR_S", "Expected_FWR_R", "Expected_FWR_S",
359 paste( rep(testName,6), rep(c("Sigma","CIlower","CIupper"),2),rep(c("CDR","FWR"),each=3), sep="_"),
360 paste( rep(testName,2), rep("P",2),c("CDR","FWR"), sep="_")
361 )
362 fileName = paste(outputPath,outputID,".txt",sep="")
363 write.table(matOutput,file=fileName,quote=F,sep="\t",row.names=T,col.names=NA)
364 fileName = paste(outputPath,outputID,".RData",sep="")
365 save(listPDFs,file=fileName)
366
367 indelWarning = FALSE
368 if(sum(indelPos)>0){
369 indelWarning = "<P>Warning: The following sequences have either gaps and/or deletions, and have been ommited from the analysis.";
370 indelWarning = paste( indelWarning , "<UL>", sep="" )
371 for(indels in names(indelPos)[indelPos]){
372 indelWarning = paste( indelWarning , "<LI>", indels, "</LI>", sep="" )
373 }
374 indelWarning = paste( indelWarning , "</UL></P>", sep="" )
375 }
376
377 cloneWarning = FALSE
378 if(clonal==1){
379 if(sum(matInputErrors)>0){
380 cloneWarning = "<P>Warning: The following clones have sequences of unequal length.";
381 cloneWarning = paste( cloneWarning , "<UL>", sep="" )
382 for(clone in names(matInputErrors)[matInputErrors]){
383 cloneWarning = paste( cloneWarning , "<LI>", names(germlines)[as.numeric(clone)], "</LI>", sep="" )
384 }
385 cloneWarning = paste( cloneWarning , "</UL></P>", sep="" )
386 }
387 }
388 cat(paste("Success",outputID,indelWarning,cloneWarning,sep="|"))