Mercurial > repos > davidvanzessen > mutation_analysis

changeset 90:f0e8dac22c6e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Wed, 01 Jun 2016 05:03:24 -0400
parents 480fdd383fdb
children 5e237c243088
files merge_and_filter.r sequence_overview.r sequence_overview.r.old wrapper.sh
diffstat 4 files changed, 242 insertions(+), 37 deletions(-) [+]
line wrap: on
line diff
--- 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
--- 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)
--- /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("<td>", val, "</td>", sep="") }
+tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) }
+make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") }
+tbl = function(df) { res = "<table border='1'>"; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "</table>"); }
+
+cat("<table border='1'>", file=main.html, append=F)
+cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
+cat("<tr><th>Sequence</th><th>Functionality</th><th>ca1</th><th>ca2</th><th>cg1</th><th>cg2</th><th>cg3</th><th>cg4</th><th>cm</th></tr>", 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("</table>", 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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- 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 ----------------<br />" >> $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 ----------------<br />" >> $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 "<table border='1'>" > $outdir/base_overview.html