annotate baseline/Baseline_Main.r @ 4:5ffd52fc35c4 draft

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