# HG changeset patch # User davidvanzessen # Date 1411045004 14400 # Node ID 68c6c7624ffc66fa63a2a8ff4bdb5234c9f27350 # Parent 8313c6cc65c50e4b0566e620c9f26c09ed00b901 Uploaded diff -r 8313c6cc65c5 -r 68c6c7624ffc RScript.r --- a/RScript.r Mon Sep 15 07:40:52 2014 -0400 +++ b/RScript.r Thu Sep 18 08:56:44 2014 -0400 @@ -34,7 +34,6 @@ cat("Removing duplicates", file=logfile, append=T) dat = dat[!duplicated(dat$paste),] patients = split(dat, dat$Patient, drop=T) -rm(dat) intervalReads = rev(c(0,10,25,50,100,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") @@ -91,8 +90,6 @@ read2Count = vector() locussum1 = vector() locussum2 = vector() - - print(patient) #for(iter in 1){ for(iter in 1:length(product[,1])){ threshhold = product[iter,threshholdIndex] @@ -184,11 +181,6 @@ 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))) -#patientFrequencyCount(patient1) -#lapply(patients[c(5,6,10)], FUN=patientFrequencyCount) -#lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) -#lapply(patients[c(6,7,8)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) -#lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) cat("Starting Cell Count analysis", file=logfile, append=T) @@ -196,11 +188,166 @@ 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))) -#patientResult = patientReadCount(patient1) -#lapply(patients[c(5,6,10)], FUN=patientReadCount) -#lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes") -#lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes") mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes") cat("", file=logfile, append=T) + +tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){ + onShort = "reads" + if(on == "Frequency"){ + onShort = "freq" + } + 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="") + + patientMerge = merge(patient1, patient2, by="Clone_Sequence") + patientMerge = merge(patientMerge, patient3, by="Clone_Sequence") + colnames(patientMerge)[32:length(colnames(patientMerge))] = paste(colnames(patientMerge)[32:length(colnames(patientMerge))], ".z", sep="") + res1 = vector() + res2 = vector() + res3 = 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[,paste(on, ".x", sep="")] > threshhold & patientMerge[,paste(on, ".y", sep="")] > threshhold & patientMerge[,paste(on, ".z", sep="")] > threshhold) + one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$Clone_Sequence %in% patientMerge[all,]$Clone_Sequence)) + two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$Clone_Sequence %in% patientMerge[all,]$Clone_Sequence)) + three = (grepl(V_Segment, patient3$V_Segment_Major_Gene) & grepl(J_Segment, patient3$J_Segment_Major_Gene) & patient3[,on] > threshhold & !(patient3$Clone_Sequence %in% patientMerge[all,]$Clone_Sequence)) + + 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)) + #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", "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", "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", "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(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", "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),"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)) + colnames(patientResult)[6] = oneSample + colnames(patientResult)[8] = twoSample + colnames(patientResult)[10] = threeSample + + 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("Normalized Read Count", oneSample), paste("Number of sequences", twoSample), paste("Normalized Read Count", twoSample), paste("Number of sequences", threeSample), paste("Normalized Read Count", 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() +} + +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 = dat[dat$Patient == "VanDongen_cALL_14696.1",] +two = dat[dat$Patient == "VanDongen_cALL_14696.2",] +three = dat[dat$Patient == "VanDongen_cALL_14696.3",] +tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="normalized_read_count") + +one = dat[dat$Sample == "16278_Left",] +two = dat[dat$Sample == "26402_Left",] +three = dat[dat$Sample == "26759_Left",] +tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="normalized_read_count") + +one = dat[dat$Sample == "16278_Right",] +two = dat[dat$Sample == "26402_Right",] +three = dat[dat$Sample == "26759_Right",] +tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="normalized_read_count") + + +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 = dat[dat$Patient == "VanDongen_cALL_14696.1",] +two = dat[dat$Patient == "VanDongen_cALL_14696.2",] +three = dat[dat$Patient == "VanDongen_cALL_14696.3",] +tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="Frequency", T) + +one = dat[dat$Sample == "16278_Left",] +two = dat[dat$Sample == "26402_Left",] +three = dat[dat$Sample == "26759_Left",] +tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="Frequency", T) + +one = dat[dat$Sample == "16278_Right",] +two = dat[dat$Sample == "26402_Right",] +three = dat[dat$Sample == "26759_Right",] +tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="Frequency", T) + + + diff -r 8313c6cc65c5 -r 68c6c7624ffc wrapper.sh --- a/wrapper.sh Mon Sep 15 07:40:52 2014 -0400 +++ b/wrapper.sh Thu Sep 18 08:56:44 2014 -0400 @@ -9,7 +9,6 @@ dir="$(cd "$(dirname "$0")" && pwd)" mkdir $outputDir -echo "testtestsetset" > $outputFile Rscript --verbose $dir/RScript.r $inputFile $outputDir $outputFile $min_freq $min_cells 2>&1 cp $dir/jquery-1.11.0.min.js $outputDir @@ -143,7 +142,6 @@ echo "" >> $html echo "" >> $html done < patients.txt -rm tmp.txt html="index.html" echo "" > $html @@ -170,6 +168,126 @@ done echo "Triplets:" >> $html +while read sample1 sample2 sample3 +do + sample1="$(echo ${sample1} | tr -d '\r' | tr -d '\n')" + sample2="$(echo ${sample2} | tr -d '\r' | tr -d '\n')" + sample3="$(echo ${sample3} | tr -d '\r' | tr -d '\n')" + patient="${sample1}_${sample2}_${sample3}" + echo "$patient" + html="${patient}.html" + echo "$patient" >> "index.html" + echo "$header" > $html + oldLocus="" + tail -n+2 ${patient}_freq.txt | sed "s/>//" > tmp.txt + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "" >> $html + echo "" >> $html + echo "" >> $html + while read locus j_segment v_segment cut_off_value all one read_count1 two read_count2 three read_count3 + do + if [ "$locus" != "$oldLocus" ] ; then + echo "" >> $html + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + echo "" >> $html + echo "" >> $html + if [ "$all" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + if [ "$one" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + if [ "$two" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + if [ "$three" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + echo "" >> $html + oldLocus="$locus" + done < tmp.txt + echo "
Ig/TCR gene rearrangement typeProximal gene segmentDistal gene segmentCut off valueNumber of sequences ${patient}_AllNumber of sequences_$sample1Read Count $sample1Number of sequences_$sample2Read Count $sample2Number of sequences_$sample3Read Count $sample3
$locus$v_segment$j_segment>$cut_off_value$all$all$one$one$read_count1$two$two$read_count2$three$three$read_count3
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + + tail -n+2 ${patient}_reads.txt | sed "s/>//" > tmp.txt + echo "
" >> $html + echo "
" >> $html + echo "" >> $html + echo "" >> $html + echo "" >> $html + while read locus j_segment v_segment cut_off_value all one read_count1 two read_count2 three read_count3 + do + if [ "$locus" != "$oldLocus" ] ; then + echo "" >> $html + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + echo "" >> $html + echo "" >> $html + if [ "$all" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + if [ "$one" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + if [ "$two" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + if [ "$three" != "0" ] && [ "$cut_off_value" != "0" ] ; then + echo "" >> $html + else + echo "" >> $html + fi + echo "" >> $html + echo "" >> $html + oldLocus="$locus" + done < tmp.txt + echo "
Ig/TCR gene rearrangement typeProximal gene segmentDistal gene segmentCut off valueNumber of sequences ${patient}_AllNumber of sequences_$sample1Read Count $sample1Number of sequences_$sample2Read Count $sample2Number of sequences_$sample3Read Count $sample3
$locus$v_segment$j_segment>$cut_off_value$all$all$one$one$read_count1$two$two$read_count2$three$three$read_count3
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "
" >> $html + echo "" >> $html +done < triplets.txt +rm tmp.txt + + +html="index.html" + echo "" >> $html echo "" >> $html