# HG changeset patch
# User davidvanzessen
# Date 1505739679 14400
# Node ID 7ffd0fba8cf4b085dc037936903471d3475fd7db
# Parent  75853bceec005adbcea81ed4f7a15204bf46e3f6
Uploaded

diff -r 75853bceec00 -r 7ffd0fba8cf4 RScript.r
--- a/RScript.r	Tue Jan 17 07:24:44 2017 -0500
+++ b/RScript.r	Mon Sep 18 09:01:19 2017 -0400
@@ -18,7 +18,23 @@
 #require(xtable)
 cat("<tr><td>Reading input</td></tr>", file=logfile, append=T)
 dat = read.table(inFile, header=T, sep="\t", dec=".", fill=T, stringsAsFactors=F)
-dat = dat[,c("Patient",  "Receptor", "Sample", "Cell_Count", "Clone_Molecule_Count_From_Spikes", "Log10_Frequency", "Total_Read_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence", "Clone_Sequence")]
+
+needed_cols = c("Patient",  "Receptor", "Sample", "Cell_Count", "Clone_Molecule_Count_From_Spikes", "Log10_Frequency", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence", "Clone_Sequence")
+if(!all(needed_cols %in% names(dat))){
+	cat("Missing column(s):<br />", file=logfile, append=F)
+	missing_cols = needed_cols[!(needed_cols %in% names(dat))]
+	for(col in missing_cols){
+		cat(paste(col, "<br />"), file=logfile, append=T)
+	}
+	stop("Not all columns are present")
+}
+
+if(!("Total_Read_Count" %in% names(dat))){
+	dat$Total_Read_Count = 0
+}
+
+dat = dat[,c("Patient",  "Receptor", "Sample", "Cell_Count", "Clone_Molecule_Count_From_Spikes", "Log10_Frequency", "Total_Read_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence", "Clone_Sequence")] 
+
 dat$dsPerM = 0
 dat = dat[!is.na(dat$Patient),]
 dat$Related_to_leukemia_clone = F
@@ -76,7 +92,7 @@
 Titles = factor(Titles, levels=Titles)
 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
 
-single_patients = data.frame("Patient" = character(0),"Sample" = character(0), "on" = character(0), "Clone_Sequence" = character(0), "Frequency" = numeric(0), "normalized_read_count" = numeric(0), "V_Segment_Major_Gene" = character(0), "J_Segment_Major_Gene" = character(0), "Rearrangement" = character(0))
+single_patients = dat[NULL,]
 
 patient.merge.list = list() #cache the 'both' table, 2x speedup for more memory...
 patient.merge.list.second = list()
@@ -96,25 +112,20 @@
   onx = paste(on, ".x", sep="")
   ony = paste(on, ".y", sep="")
   splt = split(x, x$Sample, drop=T)
+  print(splt)
   type="pair"
   if(length(splt) == 1){
     print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample"))
-    splt[[2]] = data.frame("Patient" = character(0), "Receptor" = character(0), "Sample" = character(0), "Cell_Count" = numeric(0), "Clone_Molecule_Count_From_Spikes" = numeric(0), "Log10_Frequency" = numeric(0), "Total_Read_Count" = numeric(0), "dsMol_per_1e6_cells" = numeric(0), "J_Segment_Major_Gene" = character(0), "V_Segment_Major_Gene" = character(0), "Clone_Sequence" = character(0), "CDR3_Sense_Sequence" = character(0), "Related_to_leukemia_clone" = logical(0), "Frequency"= numeric(0), "normalized_read_count" = numeric(0), "paste" = character(0))
+    splt[[2]] = splt[[1]][NULL,]
     type="single"
   }
   patient1 = splt[[1]]
   patient2 = splt[[2]]
   
-  threshholdIndex = which(colnames(product) == "interval")
-  V_SegmentIndex = which(colnames(product) == "V_Segments")
-  J_SegmentIndex = which(colnames(product) == "J_Segments")
-  titleIndex = which(colnames(product) == "Titles")
-  sampleIndex = which(colnames(x) == "Sample")
-  patientIndex = which(colnames(x) == "Patient")
-  oneSample = paste(patient1[1,sampleIndex], sep="")
-  twoSample = paste(patient2[1,sampleIndex], sep="")
-  patient = paste(x[1,patientIndex])
-
+  oneSample = patient1[1,"Sample"]
+  twoSample = patient2[1,"Sample"]
+  patient = x[1,"Patient"]
+  
   switched = F
   if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
     tmp = twoSample
@@ -139,95 +150,74 @@
   }
   
   scatterplot_data_columns = c("Patient", "Sample", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge")
-  #scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns])
   scatterplot_data = patient1[NULL,scatterplot_data_columns]
-  #scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),]
-  #scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both"))
   scatterplot.data.type.factor = c(oneSample, twoSample, paste(c(oneSample, twoSample), "In Both"))
-  #scatterplot_data$type = factor(x=NULL, levels=scatterplot.data.type.factor)
   scatterplot_data$type = character(0)
   scatterplot_data$link = numeric(0)
   scatterplot_data$on = character(0)
   
-  #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge") #merge alles 'fuzzy'
   patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")[NULL,] #blegh
-
+  
   cs.exact.matches = patient1[patient1$Clone_Sequence %in% patient2$Clone_Sequence,]$Clone_Sequence
-
+  
   start.time = proc.time()
   merge.list = c()
-
+  
   if(patient %in% names(patient.merge.list)){
     patientMerge = patient.merge.list[[patient]]
     merge.list[["second"]] = patient.merge.list.second[[patient]]
     scatterplot_data = scatter_locus_data_list[[patient]]
     cat(paste("<td>", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (fetched from cache)</td></tr>", sep=""), file=logfile, append=T)
-
+    
     print(names(patient.merge.list))
   } else {
     #fuzzy matching here...
-    #merge.list = patientMerge$merge
-
-    #patient1.fuzzy = patient1[!(patient1$merge %in% merge.list),]
-    #patient2.fuzzy = patient2[!(patient2$merge %in% merge.list),]
-
+    
     patient1.fuzzy = patient1
     patient2.fuzzy = patient2
-
-    #patient1.fuzzy$merge = paste(patient1.fuzzy$V_Segment_Major_Gene, patient1.fuzzy$J_Segment_Major_Gene, patient1.fuzzy$CDR3_Sense_Sequence)
-    #patient2.fuzzy$merge = paste(patient2.fuzzy$V_Segment_Major_Gene, patient2.fuzzy$J_Segment_Major_Gene, patient2.fuzzy$CDR3_Sense_Sequence)
-
-    #patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J, patient1.fuzzy$CDR3_Sense_Sequence)
-    #patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J, patient2.fuzzy$CDR3_Sense_Sequence)
-
+    
     patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J)
     patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J)
-
-    #merge.freq.table = data.frame(table(c(patient1.fuzzy[!duplicated(patient1.fuzzy$merge),"merge"], patient2.fuzzy[!duplicated(patient2.fuzzy$merge),"merge"]))) #also remove?
-    #merge.freq.table.gt.1 = merge.freq.table[merge.freq.table$Freq > 1,]
-
-    #patient1.fuzzy = patient1.fuzzy[patient1.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
-    #patient2.fuzzy = patient2.fuzzy[patient2.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
-
+    
     patient.fuzzy = rbind(patient1.fuzzy, patient2.fuzzy)
     patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),]
-
+    
     merge.list = list()
-
+    
     merge.list[["second"]] = vector()
-	
-	link.count = 1
-	
+    
+    link.count = 1
+    
     while(nrow(patient.fuzzy) > 1){
       first.merge = patient.fuzzy[1,"merge"]
       first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"]
       first.sample = patient.fuzzy[1,"Sample"]
       merge.filter = first.merge == patient.fuzzy$merge
-
+      
       #length.filter = nchar(patient.fuzzy$Clone_Sequence) - nchar(first.clone.sequence) <= 9
-
+      
       first.sample.filter = first.sample == patient.fuzzy$Sample
       second.sample.filter = first.sample != patient.fuzzy$Sample
-
+      
       #first match same sample, sum to a single row, same for other sample
       #then merge rows like 'normal'
-
+      
       sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence)
-
-
-
+      
+      
+      
       #match.filter = merge.filter & grepl(first.clone.sequence, patient.fuzzy$Clone_Sequence) & length.filter & sample.filter
       first.match.filter = merge.filter & sequence.filter & first.sample.filter
       second.match.filter = merge.filter & sequence.filter & second.sample.filter
-
+      
       first.rows = patient.fuzzy[first.match.filter,]
       second.rows = patient.fuzzy[second.match.filter,]
-
+      
       first.rows.v = table(first.rows$V_Segment_Major_Gene)
       first.rows.v = names(first.rows.v[which.max(first.rows.v)])
       first.rows.j = table(first.rows$J_Segment_Major_Gene)
       first.rows.j = names(first.rows.j[which.max(first.rows.j)])
-
+      
       first.sum = data.frame(merge = first.clone.sequence,
                              Patient = patient,
                              Receptor = first.rows[1,"Receptor"],
@@ -249,61 +239,65 @@
                              normalized_read_count = sum(first.rows$normalized_read_count),
                              paste = first.rows[1,"paste"],
                              min_cell_paste = first.rows[1,"min_cell_paste"])
-
+      
       if(nrow(second.rows) > 0){
         second.rows.v = table(second.rows$V_Segment_Major_Gene)
         second.rows.v = names(second.rows.v[which.max(second.rows.v)])
         second.rows.j = table(second.rows$J_Segment_Major_Gene)
         second.rows.j = names(second.rows.j[which.max(second.rows.j)])
-
+        
         second.sum = data.frame(merge = first.clone.sequence,
-                               Patient = patient,
-                               Receptor = second.rows[1,"Receptor"],
-                               Sample = second.rows[1,"Sample"],
-                               Cell_Count = second.rows[1,"Cell_Count"],
-                               Clone_Molecule_Count_From_Spikes = sum(second.rows$Clone_Molecule_Count_From_Spikes),
-                               Log10_Frequency = log10(sum(second.rows$Frequency)),
-                               Total_Read_Count = sum(second.rows$Total_Read_Count),
-                               dsPerM = sum(second.rows$dsPerM),
-                               J_Segment_Major_Gene = second.rows.j,
-                               V_Segment_Major_Gene = second.rows.v,
-                               Clone_Sequence = first.clone.sequence,
-                               CDR3_Sense_Sequence = second.rows[1,"CDR3_Sense_Sequence"],
-                               Related_to_leukemia_clone = F,
-                               Frequency = sum(second.rows$Frequency),
-                               locus_V = second.rows[1,"locus_V"],
-                               locus_J = second.rows[1,"locus_J"],
-                               min_cell_count = second.rows[1,"min_cell_count"],
-                               normalized_read_count = sum(second.rows$normalized_read_count),
-                               paste = second.rows[1,"paste"],
-                               min_cell_paste = second.rows[1,"min_cell_paste"])
-
+                                Patient = patient,
+                                Receptor = second.rows[1,"Receptor"],
+                                Sample = second.rows[1,"Sample"],
+                                Cell_Count = second.rows[1,"Cell_Count"],
+                                Clone_Molecule_Count_From_Spikes = sum(second.rows$Clone_Molecule_Count_From_Spikes),
+                                Log10_Frequency = log10(sum(second.rows$Frequency)),
+                                Total_Read_Count = sum(second.rows$Total_Read_Count),
+                                dsPerM = sum(second.rows$dsPerM),
+                                J_Segment_Major_Gene = second.rows.j,
+                                V_Segment_Major_Gene = second.rows.v,
+                                Clone_Sequence = first.clone.sequence,
+                                CDR3_Sense_Sequence = second.rows[1,"CDR3_Sense_Sequence"],
+                                Related_to_leukemia_clone = F,
+                                Frequency = sum(second.rows$Frequency),
+                                locus_V = second.rows[1,"locus_V"],
+                                locus_J = second.rows[1,"locus_J"],
+                                min_cell_count = second.rows[1,"min_cell_count"],
+                                normalized_read_count = sum(second.rows$normalized_read_count),
+                                paste = second.rows[1,"paste"],
+                                min_cell_paste = second.rows[1,"min_cell_paste"])
+        
+        #print(names(patientMerge))
+        #print(merge(first.sum, second.sum, by="merge"))		   
         patientMerge = rbind(patientMerge, merge(first.sum, second.sum, by="merge"))
+        #print("test2")
         patient.fuzzy = patient.fuzzy[!(first.match.filter | second.match.filter),]
-
+        
         hidden.clone.sequences = c(first.rows[-1,"Clone_Sequence"], second.rows[second.rows$Clone_Sequence != first.clone.sequence,"Clone_Sequence"])
         merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
-
+        
         tmp.rows = rbind(first.rows, second.rows)
+        #print("test3")
         tmp.rows = tmp.rows[order(nchar(tmp.rows$Clone_Sequence)),]
         
         
         #add to the scatterplot data
         scatterplot.row = first.sum[,scatterplot_data_columns]
-		scatterplot.row$type = paste(first.sum[,"Sample"], "In Both")
-		scatterplot.row$link = link.count
-		scatterplot.row$on = onShort
-		
-		scatterplot_data = rbind(scatterplot_data, scatterplot.row)
+        scatterplot.row$type = paste(first.sum[,"Sample"], "In Both")
+        scatterplot.row$link = link.count
+        scatterplot.row$on = onShort
+        
+        scatterplot_data = rbind(scatterplot_data, scatterplot.row)
         
         scatterplot.row = second.sum[,scatterplot_data_columns]
-		scatterplot.row$type = paste(second.sum[,"Sample"], "In Both")
-		scatterplot.row$link = link.count
-		scatterplot.row$on = onShort
-		
-		scatterplot_data = rbind(scatterplot_data, scatterplot.row)    
-		
-		#write some information about the match to a log file
+        scatterplot.row$type = paste(second.sum[,"Sample"], "In Both")
+        scatterplot.row$link = link.count
+        scatterplot.row$on = onShort
+        
+        scatterplot_data = rbind(scatterplot_data, scatterplot.row)    
+        
+        #write some information about the match to a log file
         if (nrow(first.rows) > 1 | nrow(second.rows) > 1) {
           cat(paste("<tr><td>", patient, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="multiple_matches.html", append=T)
         } else {
@@ -314,7 +308,7 @@
             #cat(paste("<tr><td>", patient, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="single_matches.html", append=T)
           }
         }
-
+        
       } else if(nrow(first.rows) > 1) {
         if(patient1[1,"Sample"] == first.sample){
           patient1 = patient1[!(patient1$Clone_Sequence %in% first.rows$Clone_Sequence),]
@@ -323,31 +317,31 @@
           patient2 = patient2[!(patient2$Clone_Sequence %in% first.rows$Clone_Sequence),]
           patient2 = rbind(patient2, first.sum)
         }
-
+        
         hidden.clone.sequences = c(first.rows[-1,"Clone_Sequence"])
         merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
-
+        
         patient.fuzzy = patient.fuzzy[-first.match.filter,]
         
         #add to the scatterplot data
         scatterplot.row = first.sum[,scatterplot_data_columns]
-		scatterplot.row$type = first.sum[,"Sample"]
-		scatterplot.row$link = link.count
-		scatterplot.row$on = onShort
-		
-		scatterplot_data = rbind(scatterplot_data, scatterplot.row)
-
+        scatterplot.row$type = first.sum[,"Sample"]
+        scatterplot.row$link = link.count
+        scatterplot.row$on = onShort
+        
+        scatterplot_data = rbind(scatterplot_data, scatterplot.row)
+        
         cat(paste("<tr bgcolor='#DDF'><td>", patient, " row ", 1:nrow(first.rows), "</td><td>", first.rows$Sample, ":</td><td>", first.rows$Clone_Sequence, "</td><td>", first.rows$normalized_read_count, "</td></tr>", sep=""), file="single_matches.html", append=T)
       } else {
         patient.fuzzy = patient.fuzzy[-1,]
         
         #add to the scatterplot data
         scatterplot.row = first.sum[,scatterplot_data_columns]
-		scatterplot.row$type = first.sum[,"Sample"]
-		scatterplot.row$link = link.count
-		scatterplot.row$on = onShort
-		
-		scatterplot_data = rbind(scatterplot_data, scatterplot.row)
+        scatterplot.row$type = first.sum[,"Sample"]
+        scatterplot.row$link = link.count
+        scatterplot.row$on = onShort
+        
+        scatterplot_data = rbind(scatterplot_data, scatterplot.row)
       }
       link.count = link.count + 1    
     }
@@ -360,10 +354,10 @@
     scatter_locus_data_list[[patient]] <<- scatterplot_data
     cat(paste("<td>", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (finding both took ", (proc.time() - start.time)[[3]], "s)</td></tr>", sep=""), file=logfile, append=T)
   }
-
+  
   patient1 = patient1[!(patient1$Clone_Sequence %in% patient.merge.list.second[[patient]]),]
   patient2 = patient2[!(patient2$Clone_Sequence %in% patient.merge.list.second[[patient]]),]
-
+  
   
   patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony])
   #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony])
@@ -377,9 +371,9 @@
   
   #for(iter in 1){
   for(iter in 1:length(product[,1])){
-    threshhold = product[iter,threshholdIndex]
-    V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
-    J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
+    threshhold = product[iter,"interval"]
+    V_Segment = paste(".*", as.character(product[iter,"V_Segments"]), ".*", sep="")
+    J_Segment = paste(".*", as.character(product[iter,"J_Segments"]), ".*", sep="")
     #both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,onx] > threshhold & patientMerge[,ony] > threshhold) #both higher than threshold
     both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold) #highest of both is higher than threshold
     one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$merge %in% patientMerge[both,]$merge))
@@ -392,50 +386,50 @@
     locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
     locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count))
     #threshhold = 0
-    if(threshhold != 0){
+    if(threshhold != 0 | T){
       if(sum(one) > 0){
         dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
         colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
-        filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
+        filenameOne = paste(oneSample, "_", product[iter, "Titles"], "_", threshhold, sep="")
         write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
       }
       if(sum(two) > 0){
         dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
         colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
-        filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
+        filenameTwo = paste(twoSample, "_", product[iter, "Titles"], "_", threshhold, sep="")
         write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
       }
     } else {
       scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
       if(nrow(scatterplot_locus_data) > 0){
-        scatterplot_locus_data$Rearrangement = product[iter, titleIndex]
+        scatterplot_locus_data$Rearrangement = product[iter, "Titles"]
       }
       
       
-            
+      
       p = NULL
-      print(paste("nrow scatterplot_locus_data", nrow(scatterplot_locus_data)))
+      #print(paste("nrow scatterplot_locus_data", nrow(scatterplot_locus_data)))
       if(nrow(scatterplot_locus_data) != 0){
         if(on == "normalized_read_count"){
-		  write.table(scatterplot_locus_data, file=paste(oneSample, twoSample, product[iter, titleIndex], "scatterplot_locus_data.txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
+          write.table(scatterplot_locus_data, file=paste(oneSample, twoSample, product[iter, "Titles"], "scatterplot_locus_data.txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
           scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
           p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), normalized_read_count, group=link)) + geom_line() + scale_y_log10(breaks=scales,labels=scales, limits=c(1,1e6)) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
         } else {
           p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), Frequency, group=link)) + geom_line() + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
         }
         p = p + geom_point(aes(colour=type), position="dodge")
-        p = p + xlab("In one or both samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex]))
+        p = p + xlab("In one or both samples") + ylab(onShort) + ggtitle(paste(patient1[1,"Patient"], patient1[1,"Sample"], patient2[1,"Sample"], onShort, product[iter, "Titles"]))
       } else {
-        p = ggplot(NULL, aes(x=c("In one", "In Both"),y=0)) + geom_blank(NULL) + xlab("In one or both of the samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex]))
+        p = ggplot(NULL, aes(x=c("In one", "In Both"),y=0)) + geom_blank(NULL) + xlab("In one or both of the samples") + ylab(onShort) + ggtitle(paste(patient1[1,"Patient"], patient1[1,"Sample"], patient2[1,"Sample"], onShort, product[iter, titleIndex]))
       }
-      png(paste(patient1[1,patientIndex], "_", patient1[1,sampleIndex], "_", patient2[1,sampleIndex], "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep=""))
+      png(paste(patient1[1,"Patient"], "_", patient1[1,"Sample"], "_", patient2[1,"Sample"], "_", onShort, "_", product[iter, "Titles"],"_scatter.png", sep=""))
       print(p)
       dev.off()
     }
     if(sum(both) > 0){
       dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
       colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample))
-      filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
+      filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, "Titles"], "_", threshhold, sep="")
       write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
     } 
   }
@@ -487,20 +481,27 @@
   print(plt)
   dev.off()
 }
+
 if(length(patients) > 0){
 	cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
 
 	interval = intervalFreq
 	intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
 	product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
-	lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
+	for (current_patient in patients){
+		print(paste("Started working", unique(current_patient$Patient), "Frequency analysis"))
+		patientCountOnColumn(current_patient, product=product, interval=interval, on="Frequency", appendtxt=T)
+	}
 
 	cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
 
 	interval = intervalReads
 	intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
 	product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
-	lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count")
+	for (current_patient in patients){
+		print(paste("Started working on ", unique(current_patient$Patient), "Read Count analysis"))
+		patientCountOnColumn(current_patient, product=product, interval=interval, on="normalized_read_count")
+	}
 }
 if(nrow(single_patients) > 0){
 	scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
@@ -936,18 +937,20 @@
 		}
 		p = NULL
 		if(nrow(scatterplot_locus_data) != 0){
-		if(on == "normalized_read_count"){
-		scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
-		p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales, limits=c(1, 1e6))
+			filename.scatter = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_scatter_", threshhold, sep="")
+			write.table(scatterplot_locus_data, file=paste(filename.scatter, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
+			if(on == "normalized_read_count"){
+			scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
+			p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales, limits=c(1, 1e6))
+			} else {
+			p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + expand_limits(y=c(0,100))
+			#p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100))
+			}
+			p = p + geom_point(aes(colour=type), position="jitter")
+			p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
 		} else {
-		p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + expand_limits(y=c(0,100))
-		#p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100))
+			p = ggplot(NULL, aes(x=c("In one", "In multiple"),y=0)) + geom_blank(NULL) + xlab("In two or in three of the samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
 		}
-		p = p + geom_point(aes(colour=type), position="jitter")
-		p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
-	  } else {
-		p = ggplot(NULL, aes(x=c("In one", "In multiple"),y=0)) + geom_blank(NULL) + xlab("In two or in three of the samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
-	  }
 		png(paste(label1, "_", label2, "_", label3, "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep=""))
 		print(p)
 	  dev.off()