# HG changeset patch # User davidvanzessen # Date 1464771804 14400 # Node ID f0e8dac22c6e8a16178bc5cca37c3752451027ba # Parent 480fdd383fdbec050e5d02fc9948e4d522724688 Uploaded diff -r 480fdd383fdb -r f0e8dac22c6e merge_and_filter.r --- a/merge_and_filter.r Tue May 31 08:30:50 2016 -0400 +++ b/merge_and_filter.r Wed Jun 01 05:03:24 2016 -0400 @@ -8,12 +8,13 @@ hotspotsfile = args[5] gene_identification_file= args[6] output = args[7] -unmatchedfile = args[8] -method=args[9] -functionality=args[10] -unique_type=args[11] -filter_unique=args[12] -class_filter=args[13] +before.unique.file = args[8] +unmatchedfile = args[9] +method=args[10] +functionality=args[11] +unique_type=args[12] +filter_unique=args[13] +class_filter=args[14] summ = read.table(summaryfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") sequences = read.table(sequencesfile, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") @@ -109,9 +110,21 @@ print(paste("Number of sequences in result after", unique_type, "filtering:", nrow(result))) #remove the sequences that have an 'n' (or 'N') in the FR2, FR3, CDR1 and CDR2 regions. -sequences = sequences[grepl("n|N", sequences$FR2.IMGT) | grepl("n|N", sequences$FR3.IMGT) | grepl("n|N", sequences$CDR1.IMGT) | grepl("n|N", sequences$CDR2.IMGT),] +sequences = sequences[,c("Sequence.ID", "FR1.IMGT", "CDR1.IMGT", "FR2.IMGT", "CDR2.IMGT", "FR3.IMGT", "CDR3.IMGT")] +names(sequences) = c("Sequence.ID", "FR1.IMGT.seq", "CDR1.IMGT.seq", "FR2.IMGT.seq", "CDR2.IMGT.seq", "FR3.IMGT.seq", "CDR3.IMGT.seq") +result = merge(result, sequences, by="Sequence.ID", all.x=T) + +print(paste("Number of sequences in result after merging with sequences:", nrow(result))) -result = result[!(result$Sequence.ID %in% sequences$Sequence.ID),] +print(paste("Number of N in CDR1:", sum(grepl("n|N", result$CDR1.IMGT.seq)))) +print(paste("Number of N in FR2:", sum(grepl("n|N", result$FR2.IMGT.seq)))) +print(paste("Number of N in CDR2:", sum(grepl("n|N", result$CDR2.IMGT.seq)))) +print(paste("Number of N in FR3:", sum(grepl("n|N", result$FR3.IMGT.seq)))) + +print(paste("Number of sequences with N in CDR1 or FR2 or CDR2 or FR3:", sum(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq)))) + + +result = result[!(grepl("n|N", result$FR2.IMGT.seq) | grepl("n|N", result$FR3.IMGT.seq) | grepl("n|N", result$CDR1.IMGT.seq) | grepl("n|N", result$CDR2.IMGT.seq)),] print(paste("Number of sequences in result after n filtering:", nrow(result))) @@ -132,15 +145,15 @@ result$JGene = gsub("^Homsap ", "", result$J.GENE.and.allele) result$JGene = gsub("[*].*", "", result$JGene) -write.table(result, "before_unique_filter.txt", sep="\t",quote=F,row.names=F,col.names=T) +write.table(result, before.unique.file, sep="\t", quote=F,row.names=F,col.names=T) if(filter_unique != "no"){ #clmns = names(result) if(grepl("_c", filter_unique)){ - result$unique.def = paste(result$CDR1.IMGT, result$FR2.IMGT, result$CDR2.IMGT, result$FR3.IMGT, result$CDR3.IMGT, result$best_match) + result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq, result$best_match) } else { - result$unique.def = paste(result$CDR1.IMGT, result$FR2.IMGT, result$CDR2.IMGT, result$FR3.IMGT, result$CDR3.IMGT) + result$unique.def = paste(result$CDR1.IMGT.seq, result$FR2.IMGT.seq, result$CDR2.IMGT.seq, result$FR3.IMGT.seq, result$CDR3.IMGT.seq) } #fltr = result$unique.def %in% result.filtered$unique.def diff -r 480fdd383fdb -r f0e8dac22c6e sequence_overview.r --- a/sequence_overview.r Tue May 31 08:30:50 2016 -0400 +++ b/sequence_overview.r Wed Jun 01 05:03:24 2016 -0400 @@ -2,36 +2,27 @@ args <- commandArgs(trailingOnly = TRUE) -gene.matches = args[1] -sequence.file = args[2] -merged.file = args[3] -outputdir = args[4] -gene.classes = unlist(strsplit(args[5], ",")) -hotspot.analysis.sum.file = args[6] +input.file = args[1] +outputdir = args[2] +gene.classes = unlist(strsplit(args[3], ",")) +hotspot.analysis.sum.file = args[4] NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") NTsum.file = paste(outputdir, "ntsum.txt", sep="/") main.html = "index.html" setwd(outputdir) -genes = read.table(gene.matches, header=T, sep="\t", fill=T) -sequences = read.table(sequence.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") -merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") +merged = read.table(input.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="") -dat = merge(sequences, genes, by="Sequence.ID") - -dat = dat[dat$Sequence.ID %in% merged$Sequence.ID,] +merged$seq_conc = paste(merged$CDR1.IMGT.seq, merged$FR2.IMGT.seq, merged$CDR2.IMGT.seq, merged$FR3.IMGT.seq, merged$CDR3.IMGT.seq) -dat$seq_conc = paste(dat$CDR1.IMGT, dat$FR2.IMGT, dat$CDR2.IMGT, dat$FR3.IMGT, dat$CDR3.IMGT) -#dat$seq_conc = paste(dat$CDR1.IMGT, dat$CDR2.IMGT, dat$CDR3.IMGT) - -IDs = dat[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] +IDs = merged[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] IDs$best_match = as.character(IDs$best_match) #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")]) -dat = data.frame(table(dat$seq_conc, dat$Functionality)) +dat = data.frame(table(merged$seq_conc, merged$Functionality)) #dat = dat[dat$Freq > 1,] @@ -124,18 +115,16 @@ -NToverview = genes[,c("Sequence.ID", "best_match")] -sequences$seq = paste(sequences$CDR1.IMGT, sequences$FR2.IMGT, sequences$CDR2.IMGT, sequences$FR3.IMGT, sep="_") - -NToverview = merge(NToverview, sequences[,c("Sequence.ID", "seq")], by="Sequence.ID") - -NToverview = NToverview[NToverview$Sequence.ID %in% merged$Sequence.ID,] +NToverview = merged +NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_") NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq)) NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq)) NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq)) NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq)) +print(sum(colSums(NToverview[,c("A", "C", "T", "G")]))) + #Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T)) #NToverview = rbind(NToverview, NTsum) diff -r 480fdd383fdb -r f0e8dac22c6e sequence_overview.r.old --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sequence_overview.r.old Wed Jun 01 05:03:24 2016 -0400 @@ -0,0 +1,203 @@ +library(reshape2) + +args <- commandArgs(trailingOnly = TRUE) + +gene.matches = args[1] +sequence.file = args[2] +merged.file = args[3] +outputdir = args[4] +gene.classes = unlist(strsplit(args[5], ",")) +hotspot.analysis.sum.file = args[6] +NToverview.file = paste(outputdir, "ntoverview.txt", sep="/") +NTsum.file = paste(outputdir, "ntsum.txt", sep="/") +main.html = "index.html" + +setwd(outputdir) + +genes = read.table(gene.matches, header=T, sep="\t", fill=T) +sequences = read.table(sequence.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") +merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="") +hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="") + +dat = merge(sequences, genes, by="Sequence.ID") + +dat = dat[dat$Sequence.ID %in% merged$Sequence.ID,] + +dat$seq_conc = paste(dat$CDR1.IMGT, dat$FR2.IMGT, dat$CDR2.IMGT, dat$FR3.IMGT, dat$CDR3.IMGT) +#dat$seq_conc = paste(dat$CDR1.IMGT, dat$CDR2.IMGT, dat$CDR3.IMGT) + +IDs = dat[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")] +IDs$best_match = as.character(IDs$best_match) + +#dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")]) + +dat = data.frame(table(dat$seq_conc, dat$Functionality)) + +#dat = dat[dat$Freq > 1,] + +names(dat) = c("seq_conc", "Functionality", "Freq") + +dat$seq_conc = factor(dat$seq_conc) + +dat = dat[order(as.character(dat$seq_conc)),] + +#writing html from R... +td = function(val) { paste("", val, "", sep="") } +tr = function(val) { capture.output(cat("", td(val), "", sep="")) } +make.link = function(id, clss, val) { paste("", val, "", sep="") } +tbl = function(df) { res = ""; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "
"); } + +cat("", file=main.html, append=F) +cat("", file=main.html, append=T) +cat("", file=main.html, append=T) + +for(i in 1:nrow(dat)){ + ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca1",] + ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca2",] + + cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg1",] + cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg2",] + cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg3",] + cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg4",] + + cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cm",] + + classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm)) + + classes.sum = sum(classes) + + if(classes.sum == 1){ + next + } + + id = as.numeric(dat[i,"seq_conc"]) + + functionality = dat[i,"Functionality"] + + if(nrow(ca1) > 0){ + cat(tbl(ca1), file=paste("ca1_", id, ".html", sep="")) + } + + if(nrow(ca2) > 0){ + cat(tbl(ca2), file=paste("ca2_", id, ".html", sep="")) + } + + if(nrow(cg1) > 0){ + cat(tbl(cg1), file=paste("cg1_", id, ".html", sep="")) + } + + if(nrow(cg2) > 0){ + cat(tbl(cg2), file=paste("cg2_", id, ".html", sep="")) + } + + if(nrow(cg3) > 0){ + cat(tbl(cg3), file=paste("cg3_", id, ".html", sep="")) + } + + if(nrow(cg4) > 0){ + cat(tbl(cg4), file=paste("cg4_", id, ".html", sep="")) + } + + if(nrow(cm) > 0){ + cat(tbl(cm), file=paste("cm_", id, ".html", sep="")) + } + + ca1.html = make.link(id, "ca1", nrow(ca1)) + ca2.html = make.link(id, "ca2", nrow(ca2)) + + cg1.html = make.link(id, "cg1", nrow(cg1)) + cg2.html = make.link(id, "cg2", nrow(cg2)) + cg3.html = make.link(id, "cg3", nrow(cg3)) + cg4.html = make.link(id, "cg4", nrow(cg4)) + + cm.html = make.link(id, "cm", nrow(cm)) + + rw = c(as.character(dat[i,"seq_conc"]), as.character(functionality), ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html) + + cat(tr(rw), file=main.html, append=T) +} + +cat("
CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once
SequenceFunctionalityca1ca2cg1cg2cg3cg4cm
", file=main.html, append=T) + + +#ACGT overview + + + +NToverview = genes[,c("Sequence.ID", "best_match")] +sequences$seq = paste(sequences$CDR1.IMGT, sequences$FR2.IMGT, sequences$CDR2.IMGT, sequences$FR3.IMGT, sep="_") + +NToverview = merge(NToverview, sequences[,c("Sequence.ID", "seq")], by="Sequence.ID") + +NToverview = NToverview[NToverview$Sequence.ID %in% merged$Sequence.ID,] + +NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq)) +NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq)) +NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq)) +NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq)) + +#Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T)) + +#NToverview = rbind(NToverview, NTsum) + +NTresult = data.frame(nt=c("A", "C", "T", "G")) + +for(clazz in gene.classes){ + NToverview.sub = NToverview[grepl(clazz, paste("^", NToverview$best_match, sep="")),] + new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G)) + new.col.y = sum(new.col.x) + new.col.z = round(new.col.x / new.col.y * 100, 2) + + tmp = names(NTresult) + NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) + names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep="")) +} + +new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G)) +new.col.y = sum(new.col.x) +new.col.z = round(new.col.x / new.col.y * 100, 2) + +tmp = names(NTresult) +NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z)) +names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep="")) + +names(hotspot.analysis.sum) = names(NTresult) + +hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult) + +print(hotspot.analysis.sum) + +write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0") + + +write.table(NToverview, NToverview.file, quote=F, sep="\t", row.names=F, col.names=T) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff -r 480fdd383fdb -r f0e8dac22c6e wrapper.sh --- a/wrapper.sh Tue May 31 08:30:50 2016 -0400 +++ b/wrapper.sh Wed Jun 01 05:03:24 2016 -0400 @@ -76,7 +76,7 @@ echo "---------------- merge_and_filter.r ----------------" echo "---------------- merge_and_filter.r ----------------
" >> $output -Rscript $dir/merge_and_filter.r $PWD/summary.txt $PWD/sequences.txt $PWD/mutationanalysis.txt $PWD/mutationstats.txt $PWD/hotspots.txt $outdir/identified_genes.txt $outdir/merged.txt $outdir/unmatched.txt $method $functionality $unique ${filter_unique} ${class_filter} 2>&1 +Rscript $dir/merge_and_filter.r $PWD/summary.txt $PWD/sequences.txt $PWD/mutationanalysis.txt $PWD/mutationstats.txt $PWD/hotspots.txt $outdir/identified_genes.txt $outdir/merged.txt $outdir/before_unique_filter.txt $outdir/unmatched.txt $method $functionality $unique ${filter_unique} ${class_filter} 2>&1 echo "---------------- mutation_analysis.r ----------------" echo "---------------- mutation_analysis.r ----------------
" >> $output @@ -104,8 +104,8 @@ mkdir $outdir/sequence_overview -Rscript $dir/sequence_overview.r $outdir/identified_genes.txt $PWD/sequences.txt $outdir/merged.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt 2>&1 -#Rscript $dir/sequence_overview.r $outdir/before_unique_filter.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt 2>&1 +#Rscript $dir/sequence_overview.r $outdir/identified_genes.txt $PWD/sequences.txt $outdir/merged.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt 2>&1 +Rscript $dir/sequence_overview.r $outdir/before_unique_filter.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt 2>&1 echo "" > $outdir/base_overview.html