Mercurial > repos > davidvanzessen > shm_csr

changeset 7:ad9be244b104 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 07 Nov 2016 03:04:07 -0500
parents 2ddb9a21f635
children 3968d04b5724
files sequence_overview.r shm_csr.py shm_csr.r wrapper.sh
diffstat 4 files changed, 57 insertions(+), 17 deletions(-) [+]
line wrap: on
line diff
--- a/sequence_overview.r	Tue Nov 01 10:48:38 2016 -0400
+++ b/sequence_overview.r	Mon Nov 07 03:04:07 2016 -0500
@@ -10,6 +10,8 @@
 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
 main.html = "index.html"
+empty.region.filter = args[6]
+
 
 setwd(outputdir)
 
@@ -19,19 +21,21 @@
 
 #before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]
 
-before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
+if(empty.region.filter == "leader"){
+	before.unique$seq_conc = paste(before.unique$FR1.IMGT.seq, before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
+} else if(empty.region.filter == "FR1"){
+	before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
+} else if(empty.region.filter == "CDR1"){
+	before.unique$seq_conc = paste(before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
+} else if(empty.region.filter == "FR2"){
+	before.unique$seq_conc = paste(before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq)
+}
 
 IDs = before.unique[,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(before.unique$seq_conc))
 
-dat = data.frame(table(before.unique$seq_conc))
-#dat = data.frame(table(merged$seq_conc, merged$Functionality))
-
-#dat = dat[dat$Freq > 1,]
-
-#names(dat) = c("seq_conc", "Functionality", "Freq")
 names(dat) = c("seq_conc", "Freq")
 
 dat$seq_conc = factor(dat$seq_conc)
@@ -67,7 +71,17 @@
 }
 
 cat("<table border='1' class='pure-table pure-table-striped'>", 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)
+
+if(empty.region.filter == "leader"){
+	cat("<caption>FR1+CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
+} else if(empty.region.filter == "FR1"){
+	cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
+} else if(empty.region.filter == "CDR1"){
+	cat("<caption>FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
+} else if(empty.region.filter == "FR2"){
+	cat("<caption>CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
+}
+
 cat("<tr>", file=main.html, append=T)
 cat("<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><th>un</th>", file=main.html, append=T)
 cat("<th>total CA</th><th>total CG</th><th>number of subclasses</th><th>present in both Ca and Cg</th><th>Ca1+Ca2</th>", file=main.html, append=T)
@@ -240,9 +254,18 @@
 
 #ACGT overview
 
-NToverview = merged[!grepl("^unmatched", merged$best_match),]
+#NToverview = merged[!grepl("^unmatched", merged$best_match),]
+NToverview = merged
 
-NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_")
+if(empty.region.filter == "leader"){
+	NToverview$seq = paste(NToverview$FR1.IMGT.seq, NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
+} else if(empty.region.filter == "FR1"){
+	NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
+} else if(empty.region.filter == "CDR1"){
+	NToverview$seq = paste(NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
+} else if(empty.region.filter == "FR2"){
+	NToverview$seq = paste(NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq)
+}
 
 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
--- a/shm_csr.py	Tue Nov 01 10:48:38 2016 -0400
+++ b/shm_csr.py	Mon Nov 07 03:04:07 2016 -0500
@@ -247,9 +247,12 @@
 	
 
 def get_xyz(lst, gene, f, fname):
+	
 	x = int(round(f(lst)))
 	y = valuedic[gene + "_" + fname]
 	z = str(round(x / float(y) * 100, 1)) if y != 0 else "0"
+	if gene == "unmatched":
+		print x, y, z
 	return (str(x), str(y), z)
 
 dic = {"RGYW": RGYWCount, "WRCY": WRCYCount, "WA": WACount, "TW": TWCount}
@@ -271,8 +274,8 @@
 				else:
 					x, y, z = get_xyz([curr[x] for x in [y for y, z in genedic.iteritems() if geneMatcher.match(z)]], gene, func, fname)
 					o.write("," + x + "," + y + "," + z)
-
 			x, y, z = get_xyz([y for x, y in curr.iteritems() if not genedic[x].startswith("unmatched")], "total", func, fname)
+			#x, y, z = get_xyz([y for x, y in curr.iteritems()], "total", func, fname)
 			o.write("," + x + "," + y + "," + z + "\n")
 
 
--- a/shm_csr.r	Tue Nov 01 10:48:38 2016 -0400
+++ b/shm_csr.r	Mon Nov 07 03:04:07 2016 -0500
@@ -119,9 +119,9 @@
 if(empty.region.filter == "FR1") {
 	regions = c("CDR1", "FR2", "CDR2", "FR3")
 } else if (empty.region.filter == "CDR1") {
-	regions = c("FR2", "CDR2", "FR3", "CDR3")
+	regions = c("FR2", "CDR2", "FR3")
 } else if (empty.region.filter == "FR2") {
-	regions = c("CDR2", "FR3", "CDR3")
+	regions = c("CDR2", "FR3")
 }
 
 sum_by_row = function(x, columns) { sum(as.numeric(x[columns]), na.rm=T) }
@@ -229,11 +229,25 @@
 		matrx[9,z] = round(matrx[9,x] / matrx[9,y], digits=1)
 
 		if(fname == "sum"){
-			matrx[10,x] = round(f(rowSums(tmp[,c("FR2.IMGT.Nb.of.nucleotides", "FR3.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
+			
+			regions.fr = regions[grepl("FR", regions)]
+			regions.fr = paste(regions.fr, ".IMGT.Nb.of.nucleotides", sep="")
+			regions.cdr = regions[grepl("CDR", regions)]
+			regions.cdr = paste(regions.cdr, ".IMGT.Nb.of.nucleotides", sep="")
+			
+			if(length(regions.fr) > 1){ #in case there is only on FR region (rowSums needs >1 column)
+				matrx[10,x] = round(f(rowSums(tmp[,regions.fr], na.rm=T)), digits=1)
+			} else {
+				matrx[10,x] = round(f(tmp[,regions.fr], na.rm=T), digits=1)
+			}
 			matrx[10,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
 			matrx[10,z] = round(matrx[10,x] / matrx[10,y] * 100, digits=1)
 
-			matrx[11,x] = round(f(rowSums(tmp[,c("CDR1.IMGT.Nb.of.nucleotides", "CDR2.IMGT.Nb.of.nucleotides")], na.rm=T)), digits=1)
+			if(length(regions.cdr) > 1){ #in case there is only on CDR region
+				matrx[11,x] = round(f(rowSums(tmp[,regions.cdr], na.rm=T)), digits=1)
+			} else {
+				matrx[11,x] = round(f(tmp[,regions.cdr], na.rm=T), digits=1)
+			}
 			matrx[11,y] = round(f(tmp$VRegionNucleotides, na.rm=T), digits=1)
 			matrx[11,z] = round(matrx[11,x] / matrx[11,y] * 100, digits=1)
 		}
--- a/wrapper.sh	Tue Nov 01 10:48:38 2016 -0400
+++ b/wrapper.sh	Mon Nov 07 03:04:07 2016 -0500
@@ -197,7 +197,7 @@
 
 mkdir $outdir/sequence_overview
 
-Rscript $dir/sequence_overview.r $outdir/before_unique_filter.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/merged.txt $outdir/sequence_overview $classes $outdir/hotspot_analysis_sum.txt ${empty_region_filter} 2>&1
 
 echo "<table border='1'>" > $outdir/base_overview.html