# HG changeset patch # User davidvanzessen # Date 1465216014 14400 # Node ID b869a126e2c4b3d199b935ee0f8521b393eb11f7 # Parent 5e237c2430881ffce6914a85c18de86b92d4ae18 Uploaded diff -r 5e237c243088 -r b869a126e2c4 gene_identification.py --- a/gene_identification.py Fri Jun 03 10:35:52 2016 -0400 +++ b/gene_identification.py Mon Jun 06 08:26:54 2016 -0400 @@ -160,9 +160,9 @@ chunksInCM = len(compiledregex["cm"]) requiredChunkPercentage = 0.7 varsInCA = float(len(ca1.keys()) * 2) -varsInCG = float(len(cg1.keys()) * 2) - 1 # -1 because the sliding window doesn't hit the first nt twice +varsInCG = float(len(cg1.keys()) * 2) - 2 # -1 because the sliding window doesn't hit the first nt twice varsInCM = 0 -requiredVarPercentage = 0.7 + first = True diff -r 5e237c243088 -r b869a126e2c4 sequence_overview.r --- a/sequence_overview.r Fri Jun 03 10:35:52 2016 -0400 +++ b/sequence_overview.r Mon Jun 06 08:26:54 2016 -0400 @@ -66,7 +66,7 @@ id = as.numeric(dat[i,"seq_conc"]) - functionality = paste(unique(allc[,"Functionality"], sep=",")) + functionality = paste(unique(allc[,"Functionality"]), collapse=",") if(nrow(ca1) > 0){ cat(tbl(ca1), file=paste("ca1_", id, ".html", sep="")) @@ -106,8 +106,7 @@ 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) - #print(rw) + rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html) cat(tr(rw), file=main.html, append=T) } @@ -117,8 +116,6 @@ #ACGT overview - - NToverview = merged NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_") diff -r 5e237c243088 -r b869a126e2c4 sequence_overview.r.old --- a/sequence_overview.r.old Fri Jun 03 10:35:52 2016 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,203 +0,0 @@ -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
Sequence	Functionality	ca1	ca2	cg1	cg2	cg3	cg4	cm

", 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 5e237c243088 -r b869a126e2c4 wrapper.sh --- a/wrapper.sh Fri Jun 03 10:35:52 2016 -0400 +++ b/wrapper.sh Mon Jun 06 08:26:54 2016 -0400 @@ -55,7 +55,6 @@ echo "${BLASTN_DIR}" -echo "identification ($method)" echo "---------------- identification ($method) ----------------" echo "---------------- identification ($method) ----------------
" >> $output @@ -109,7 +108,7 @@ echo "" > $outdir/base_overview.html -while read ID class seq A C G T +while IFS=$'\t' read ID class seq A C G T do echo "" >> $outdir/base_overview.html done < $outdir/sequence_overview/ntoverview.txt

$ID

$seq

$class