Mercurial > repos > davidvanzessen > clonal_sequences_in_paired_samples

args <- commandArgs(trailingOnly = TRUE)

inFile = args[1]
outDir = args[2]
logfile = args[3]
min_freq = as.numeric(args[4])
min_cells = as.numeric(args[5])
mergeOn = args[6]

cat("<html><table><tr><td>Starting analysis</td></tr>", file=logfile, append=F)

library(ggplot2)
library(reshape2)
library(data.table)
library(grid)
library(parallel)
#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[!is.na(dat$Patient),]
dat$Related_to_leukemia_clone = F

setwd(outDir)
cat("<tr><td>Selecting first V/J Genes</td></tr>", file=logfile, append=T)
dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))

cat("<tr><td>Calculating Frequency</td></tr>", file=logfile, append=T)

dat$Frequency = ((10^dat$Log10_Frequency)*100)

dat = dat[dat$Frequency >= min_freq,]

triplets = dat[grepl("VanDongen_cALL_14696", dat$Patient) | grepl("(16278)|(26402)|(26759)", dat$Sample),]

cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T)

dat$locus_V = substring(dat$V_Segment_Major_Gene, 0, 4)
dat$locus_J = substring(dat$J_Segment_Major_Gene, 0, 4)
min_cell_count = data.frame(data.table(dat)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("Patient", "locus_V", "locus_J")])

dat$min_cell_paste = paste(dat$Patient, dat$locus_V, dat$locus_J)
min_cell_count$min_cell_paste = paste(min_cell_count$Patient, min_cell_count$locus_V, min_cell_count$locus_J)

min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]

dat = merge(dat, min_cell_count, by="min_cell_paste")

dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * dat$min_cell_count / 2, digits=2) #??????????????????????????????????? wel of geen / 2

dat = dat[dat$normalized_read_count >= min_cells,]

dat$paste = paste(dat$Sample, dat$Clone_Sequence)

cat("<tr><td>Adding duplicate V+J+CDR3 sequences</td></tr>", file=logfile, append=T)
#remove duplicate V+J+CDR3, add together numerical values
dat= data.frame(data.table(dat)[, list(Receptor=unique(.SD$Receptor),
                                        Cell_Count=unique(.SD$Cell_Count),
                                        Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes),
                                        Total_Read_Count=sum(.SD$Total_Read_Count),
                                        dsPerM=ifelse("dsPerM" %in% names(dat), sum(.SD$dsPerM), 0),
                                        Related_to_leukemia_clone=all(.SD$Related_to_leukemia_clone),
                                        Frequency=sum(.SD$Frequency),
                                        locus_V=unique(.SD$locus_V),
                                        locus_J=unique(.SD$locus_J),
                                        min_cell_count=unique(.SD$min_cell_count),
                                        normalized_read_count=sum(.SD$normalized_read_count),
                                        Log10_Frequency=sum(.SD$Log10_Frequency),
                                        Clone_Sequence=.SD$Clone_Sequence[1],
                                        min_cell_paste=.SD$min_cell_paste[1],
                                        paste=unique(.SD$paste)), by=c("Patient", "Sample", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "CDR3_Sense_Sequence")])


patients = split(dat, dat$Patient, drop=T)
intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000))
intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5))
V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
Titles = c("Total", "IGH-Vh-Jh", "IGH-Dh-Jh", "Vk-Jk", "Vk-Kde" , "Intron-Kde", "TCRG", "TCRD-Vd-Dd", "TCRD-Dd-Dd", "TCRB-Vb-Jb")
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))

patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
  if (!is.data.frame(x) & is.list(x)){
    x = x[[1]]
  }
  x$Sample = factor(x$Sample, levels=unique(x$Sample))
  onShort = "reads"
  if(on == "Frequency"){
    onShort = "freq"
  }
  onx = paste(on, ".x", sep="")
  ony = paste(on, ".y", sep="")
  splt = split(x, x$Sample, drop=T)
  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))
    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])

  switched = F
  if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
    tmp = twoSample
    twoSample = oneSample
    oneSample = tmp
    tmp = patient1
    patient1 = patient2
    patient2 = tmp
    switched = T
  }
  if(appendtxt){
    cat(paste(patient, oneSample, twoSample, type, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
  }
  cat(paste("<tr><td>", patient, "</td></tr>", sep=""), file=logfile, append=T)

  if(mergeOn == "Clone_Sequence"){
    patient1$merge = paste(patient1$Clone_Sequence)
    patient2$merge = paste(patient2$Clone_Sequence)
  } else {
    patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
    patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
  }

  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 = scatterplot_data[!duplicated(scatterplot_data$merge),]
  scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both"))
  scatterplot_data$on = onShort

  patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")
  patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony])
  res1 = vector()
  res2 = vector()
  resBoth = vector()
  read1Count = vector()
  read2Count = vector()
  locussum1 = vector()
  locussum2 = vector()

  #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="")
    #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))
    two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$merge %in% patientMerge[both,]$merge))
    read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count))
    read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count))
    res1 = append(res1, sum(one))
    res2 = append(res2, sum(two))
    resBoth = append(resBoth, sum(both))
    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(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="")
        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="")
        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]
      }
      in_both = (scatterplot_locus_data$merge %in% patientMerge[both,]$merge)
      if(any(in_both)){
        scatterplot_locus_data[in_both,]$type = "In Both"
      }
      in_one = (scatterplot_locus_data$merge %in% patient1$merge)
      not_in_one = !in_one
      if(any(not_in_one)){
				scatterplot_locus_data[not_in_one,]$type = twoSample
      }
      if(type == "single"){
        single_patients <<- rbind(single_patients, scatterplot_locus_data)
      }
      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) + expand_limits(y=c(0,1000000))
        } else {
          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 both samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex]))
      } 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]))
      }
      png(paste(patient1[1,patientIndex], "_", patient1[1,sampleIndex], "_", patient2[1,sampleIndex], "_", onShort, "_", product[iter, titleIndex],"_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="")
      write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
    }
  }
  patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "Both"=resBoth, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "Sum"=res1 + res2 + resBoth, "percentage" = round((resBoth/(res1 + res2 + resBoth)) * 100, digits=2), "Locus_sum1"=locussum1, "Locus_sum2"=locussum2)
  if(sum(is.na(patientResult$percentage)) > 0){
    patientResult[is.na(patientResult$percentage),]$percentage = 0
  }
  colnames(patientResult)[6] = oneSample
  colnames(patientResult)[8] = twoSample
  colnamesBak = colnames(patientResult)
  colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", paste("Number of sequences ", patient, "_Both", sep=""), paste("Number of sequences", oneSample, sep=""), paste("Normalized Read Count", oneSample), paste("Number of sequences", twoSample, sep=""), paste("Normalized Read Count", twoSample), paste("Sum number of sequences", patient), paste("Percentage of sequences ", patient, "_Both", sep=""), paste("Locus Sum", oneSample), paste("Locus Sum", twoSample))
  write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
  colnames(patientResult) = colnamesBak

  patientResult$Locus = factor(patientResult$Locus, Titles)
  patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))

  plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
  plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
  plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
  plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
  plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
  png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
  print(plt)
  dev.off()
  #(t,r,b,l)
  plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
  plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
  plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
  plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
  plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
  png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
  print(plt)
  dev.off()

  patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
  patientResult$relativeValue = patientResult$value * 10
  patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
  plt = ggplot(patientResult)
  plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
  plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
  plt = plt + geom_text(data=patientResult[patientResult$variable == oneSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=-0.2)
  plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.8)
  plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
  png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
  print(plt)
  dev.off()
}

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)

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

cat("</table></html>", file=logfile, append=T)

scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
p = ggplot(single_patients, aes(Rearrangement, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales) + expand_limits(y=c(0,1000000))
p = p + geom_point(aes(colour=type), position="jitter")
p = p + xlab("In one or both samples") + ylab("Reads")
p = p + facet_grid(.~Patient) + ggtitle("Scatterplot of the reads of the patients with a single sample")
png("singles_reads_scatterplot.png", width=640 * length(unique(single_patients$Patient)), height=1080)
print(p)
dev.off()

p = ggplot(single_patients, aes(Rearrangement, 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 both samples") + ylab("Frequency")
p = p + facet_grid(.~Patient) + ggtitle("Scatterplot of the frequency of the patients with a single sample")
png("singles_freq_scatterplot.png", width=640 * length(unique(single_patients$Patient)), height=1080)
print(p)
dev.off()

tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){
  onShort = "reads"
  if(on == "Frequency"){
    onShort = "freq"
  }
  onx = paste(on, ".x", sep="")
  ony = paste(on, ".y", sep="")
  onz = paste(on, ".z", sep="")
  type="triplet"

  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(patient1) == "Sample")
  patientIndex = which(colnames(patient1) == "Patient")
  oneSample = paste(patient1[1,sampleIndex], sep="")
  twoSample = paste(patient2[1,sampleIndex], sep="")
  threeSample = paste(patient3[1,sampleIndex], sep="")

  if(mergeOn == "Clone_Sequence"){
    patient1$merge = paste(patient1$Clone_Sequence)
		patient2$merge = paste(patient2$Clone_Sequence)
		patient3$merge = paste(patient3$Clone_Sequence)

  } else {
		patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
		patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
		patient3$merge = paste(patient3$V_Segment_Major_Gene, patient3$J_Segment_Major_Gene, patient3$CDR3_Sense_Sequence)
  }

  patientMerge = merge(patient1, patient2, by="merge")
  patientMerge = merge(patientMerge, patient3, by="merge")
  colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="")
  patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
  patientMerge12 = merge(patient1, patient2, by="merge")
  patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony])
  patientMerge13 = merge(patient1, patient3, by="merge")
  patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony])
  patientMerge23 = merge(patient2, patient3, by="merge")
  patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony])


  scatterplot_data_columns = c("Clone_Sequence", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge")
  scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns], patient3[,scatterplot_data_columns])
  scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),]
  scatterplot_data$type = factor(x="In one", levels=c("In one", "In two", "In three", "In multiple"))

  res1 = vector()
  res2 = vector()
  res3 = vector()
  res12 = vector()
  res13 = vector()
  res23 = vector()
  resAll = vector()
  read1Count = vector()
  read2Count = vector()
  read3Count = vector()

  if(appendTriplets){
    cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3)
  }
  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="")
    #all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,onx] > threshhold & patientMerge[,ony] > threshhold & patientMerge[,onz] > threshhold)
    all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold)

    one_two = (grepl(V_Segment, patientMerge12$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge12$J_Segment_Major_Gene.x) & patientMerge12$thresholdValue > threshhold & !(patientMerge12$merge %in% patientMerge[all,]$merge))
    one_three = (grepl(V_Segment, patientMerge13$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge13$J_Segment_Major_Gene.x) & patientMerge13$thresholdValue > threshhold & !(patientMerge13$merge %in% patientMerge[all,]$merge))
    two_three = (grepl(V_Segment, patientMerge23$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge23$J_Segment_Major_Gene.x) & patientMerge23$thresholdValue > threshhold & !(patientMerge23$merge %in% patientMerge[all,]$merge))

    one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$merge %in% patientMerge[all,]$merge) & !(patient1$merge %in% patientMerge12[one_two,]$merge) & !(patient1$merge %in% patientMerge13[one_three,]$merge))
    two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$merge %in% patientMerge[all,]$merge) & !(patient2$merge %in% patientMerge12[one_two,]$merge) & !(patient2$merge %in% patientMerge23[two_three,]$merge))
    three = (grepl(V_Segment, patient3$V_Segment_Major_Gene) & grepl(J_Segment, patient3$J_Segment_Major_Gene) & patient3[,on] > threshhold & !(patient3$merge %in% patientMerge[all,]$merge) & !(patient3$merge %in% patientMerge13[one_three,]$merge) & !(patient3$merge %in% patientMerge23[two_three,]$merge))

    read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x))
    read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y))
    read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z))
    res1 = append(res1, sum(one))
    res2 = append(res2, sum(two))
    res3 = append(res3, sum(three))
    resAll = append(resAll, sum(all))
    res12 = append(res12, sum(one_two))
    res13 = append(res13, sum(one_three))
    res23 = append(res23, sum(two_three))
    #threshhold = 0
    if(threshhold != 0){
      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(label1, "_", product[iter, titleIndex], "_", 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(label2, "_", product[iter, titleIndex], "_", threshhold, sep="")
        write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
      }
      if(sum(three) > 0){
        dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
        colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
        filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
        write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
      }
      if(sum(one_two) > 0){
        dfOne_two = patientMerge12[one_two,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(dfOne_two) = 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))
        filenameOne_two = paste(label1, "_", label2, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
        write.table(dfOne_two, file=paste(filenameOne_two, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
      }
      if(sum(one_three) > 0){
        dfOne_three = patientMerge13[one_three,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(dfOne_three) = 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", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
        filenameOne_three = paste(label1, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
        write.table(dfOne_three, file=paste(filenameOne_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
      }
      if(sum(two_three) > 0){
        dfTwo_three = patientMerge23[two_three,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(dfTwo_three) = c(paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample),"Clone_Sequence", paste("Proximal segment", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
        filenameTwo_three = paste(label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
        write.table(dfTwo_three, file=paste(filenameTwo_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
      }
    } else { #scatterplot data
      scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
      in_two = (scatterplot_locus_data$merge %in% patientMerge12[one_two,]$merge) | (scatterplot_locus_data$merge %in% patientMerge13[one_three,]$merge) | (scatterplot_locus_data$merge %in% patientMerge23[two_three,]$merge)
      if(sum(in_two) > 0){
				scatterplot_locus_data[in_two,]$type = "In two"
      }
      in_three = (scatterplot_locus_data$merge %in% patientMerge[all,]$merge)
      if(sum(in_three)> 0){
				scatterplot_locus_data[in_three,]$type = "In three"
      }
      not_in_one = scatterplot_locus_data$type != "In one"
      if(sum(not_in_one) > 0){
				scatterplot_locus_data[not_in_one,]$type = "In multiple"
      }
      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) + expand_limits(y=c(0,1000000))
        } else {
          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(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()
    }
    if(sum(all) > 0){
      dfAll = patientMerge[all,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", "V_Segment_Major_Gene.z", "J_Segment_Major_Gene.z", "normalized_read_count.z", "Frequency.z", "Related_to_leukemia_clone.z")]
      colnames(dfAll) = 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), paste("Proximal segment", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
      filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
      write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
    }
  }
  #patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "All"=resAll, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "tmp3"=res3, "read_count3"=round(read3Count))
  patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "All"=resAll, "tmp1"=res1, "tmp2"=res2, "tmp3"=res3, "tmp12"=res12, "tmp13"=res13, "tmp23"=res23)
  colnames(patientResult)[6] = oneSample
  colnames(patientResult)[7] = twoSample
  colnames(patientResult)[8] = threeSample
  colnames(patientResult)[9] = paste(oneSample, twoSample, sep="_")
  colnames(patientResult)[10] = paste(oneSample, twoSample, sep="_")
  colnames(patientResult)[11] = paste(oneSample, twoSample, sep="_")

  colnamesBak = colnames(patientResult)
  colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", "Number of sequences All", paste("Number of sequences", oneSample), paste("Number of sequences", twoSample), paste("Number of sequences", threeSample), paste("Number of sequences", oneSample, twoSample), paste("Number of sequences", oneSample, threeSample), paste("Number of sequences", twoSample, threeSample))
  write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
  colnames(patientResult) = colnamesBak

  patientResult$Locus = factor(patientResult$Locus, Titles)
  patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))

  plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")])
  plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a")
  plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0)
  plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All")
  plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
  png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080)
  print(plt)
  dev.off()

  fontSize = 4

  bak = patientResult
  patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2)
  patientResult$relativeValue = patientResult$value * 10
  patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
  plt = ggplot(patientResult)
  plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
  plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
  plt = plt + geom_text(data=patientResult[patientResult$variable == oneSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=-0.7, size=fontSize)
  plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.4, size=fontSize)
  plt = plt + geom_text(data=patientResult[patientResult$variable == threeSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=1.5, size=fontSize)
  plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample")
  png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080)
  print(plt)
  dev.off()
}

triplets$uniqueID = "ID"

triplets[grepl("16278_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
triplets[grepl("26402_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
triplets[grepl("26759_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"

triplets[grepl("16278_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
triplets[grepl("26402_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
triplets[grepl("26759_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"

triplets[grepl("14696", triplets$Patient),]$uniqueID = "14696"

triplets$locus_V = substring(triplets$V_Segment_Major_Gene, 0, 4)
triplets$locus_J = substring(triplets$J_Segment_Major_Gene, 0, 4)
min_cell_count = data.frame(data.table(triplets)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("uniqueID", "locus_V", "locus_J")])

triplets$min_cell_paste = paste(triplets$uniqueID, triplets$locus_V, triplets$locus_J)
min_cell_count$min_cell_paste = paste(min_cell_count$uniqueID, min_cell_count$locus_V, min_cell_count$locus_J)

min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]

triplets = merge(triplets, min_cell_count, by="min_cell_paste")

triplets$normalized_read_count = round(triplets$Clone_Molecule_Count_From_Spikes / triplets$Cell_Count * triplets$min_cell_count / 2, digits=2) #??????????????????????????????????? wel of geen / 2

triplets = triplets[triplets$normalized_read_count >= min_cells,]

column_drops = c("locus_V", "locus_J", "min_cell_count", "min_cell_paste")

triplets = triplets[,!(colnames(triplets) %in% column_drops)]

#remove duplicate V+J+CDR3, add together numerical values
triplets = data.frame(data.table(triplets)[, list(Receptor=unique(.SD$Receptor),
                                                 Cell_Count=unique(.SD$Cell_Count),
                                                 Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes),
                                                 Total_Read_Count=sum(.SD$Total_Read_Count),
                                                 dsPerM=ifelse("dsPerM" %in% names(dat), sum(.SD$dsPerM), 0),
                                                 Related_to_leukemia_clone=all(.SD$Related_to_leukemia_clone),
                                                 Frequency=sum(.SD$Frequency),
                                                 normalized_read_count=sum(.SD$normalized_read_count),
                                                 Log10_Frequency=sum(.SD$Log10_Frequency),
                                                 Clone_Sequence=.SD$Clone_Sequence[1]), by=c("Patient", "Sample", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "CDR3_Sense_Sequence")])


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

one = triplets[triplets$Sample == "14696_reg_BM",]
two = triplets[triplets$Sample == "24536_reg_BM",]
three = triplets[triplets$Sample == "24062_reg_BM",]
tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="normalized_read_count", T)

one = triplets[triplets$Sample == "16278_Left",]
two = triplets[triplets$Sample == "26402_Left",]
three = triplets[triplets$Sample == "26759_Left",]
tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="normalized_read_count", T)

one = triplets[triplets$Sample == "16278_Right",]
two = triplets[triplets$Sample == "26402_Right",]
three = triplets[triplets$Sample == "26759_Right",]
tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="normalized_read_count", 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)))

one = triplets[triplets$Sample == "14696_reg_BM",]
two = triplets[triplets$Sample == "24536_reg_BM",]
three = triplets[triplets$Sample == "24062_reg_BM",]
tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="Frequency", F)

one = triplets[triplets$Sample == "16278_Left",]
two = triplets[triplets$Sample == "26402_Left",]
three = triplets[triplets$Sample == "26759_Left",]
tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="Frequency", F)

one = triplets[triplets$Sample == "16278_Right",]
two = triplets[triplets$Sample == "26402_Right",]
three = triplets[triplets$Sample == "26759_Right",]
tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="Frequency", F)
author	davidvanzessen
date	Tue, 02 Jun 2015 05:33:58 -0400
parents	ce8bd23d0335
children	642b4593f0a4