--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/imgt_loader.py	Mon Nov 09 07:21:09 2015 -0500
@@ -0,0 +1,141 @@
+import pandas as pd
+try:
+	pd.options.mode.chained_assignment = None  # default='warn'
+except:
+	pass
+import re
+import argparse
+import os
+
+def stop_err( msg, ret=1 ):
+    sys.stderr.write( msg )
+    sys.exit( ret )
+
+#docs.python.org/dev/library/argparse.html
+parser = argparse.ArgumentParser()
+parser.add_argument("--summ", help="The 1_Summary file from the imgt output")
+parser.add_argument("--aa", help="The 5_AA-Sequence file from the imgt output")
+parser.add_argument("--junction", help="The 6_Junction file from the imgt output")
+parser.add_argument("--output", help="Output file")
+
+args = parser.parse_args()
+
+old_summary_columns = [u'Sequence ID', u'JUNCTION frame', u'V-GENE and allele', u'D-GENE and allele', u'J-GENE and allele', u'CDR1-IMGT length', u'CDR2-IMGT length', u'CDR3-IMGT length', u'Orientation']
+old_sequence_columns = [u'CDR1-IMGT', u'CDR2-IMGT', u'CDR3-IMGT']
+old_junction_columns = [u'JUNCTION']
+
+added_summary_columns = [u'Functionality', u'V-REGION identity %', u'V-REGION identity nt', u'D-REGION reading frame', u'AA JUNCTION', u'Functionality comment', u'Sequence']
+added_sequence_columns = [u'FR1-IMGT', u'FR2-IMGT', u'FR3-IMGT', u'CDR3-IMGT', u'JUNCTION', u'J-REGION', u'FR4-IMGT']
+added_junction_columns = [u"P3'V-nt nb", u'N1-REGION-nt nb', u"P5'D-nt nb", u"P3'D-nt nb", u'N2-REGION-nt nb', u"P5'J-nt nb", u"3'V-REGION trimmed-nt nb", u"5'D-REGION trimmed-nt nb", u"3'D-REGION trimmed-nt nb", u"5'J-REGION trimmed-nt nb"]
+
+outFile = args.output
+
+#fSummary = pd.read_csv(triplets[0][0], sep="\t", low_memory=False)
+fSummary = pd.read_csv(args.summ, sep="\t", dtype=object)
+#fSequence = pd.read_csv(triplets[0][1], sep="\t", low_memory=False)
+fSequence = pd.read_csv(args.aa, sep="\t", dtype=object)
+#fJunction = pd.read_csv(triplets[0][2], sep="\t", low_memory=False)
+fJunction = pd.read_csv(args.junction, sep="\t", dtype=object)
+tmp = fSummary[["Sequence ID", "JUNCTION frame", "V-GENE and allele", "D-GENE and allele", "J-GENE and allele"]]
+
+tmp["CDR1 Seq"] = fSequence["CDR1-IMGT"]
+tmp["CDR1 Length"] = fSummary["CDR1-IMGT length"]
+
+tmp["CDR2 Seq"] = fSequence["CDR2-IMGT"]
+tmp["CDR2 Length"] = fSummary["CDR2-IMGT length"]
+
+tmp["CDR3 Seq"] = fSequence["CDR3-IMGT"]
+tmp["CDR3 Length"] = fSummary["CDR3-IMGT length"]
+
+tmp["CDR3 Seq DNA"] = fJunction["JUNCTION"]
+tmp["CDR3 Length DNA"] = '1'
+tmp["Strand"] = fSummary["Orientation"]
+tmp["CDR3 Found How"] = 'a'
+
+for col in added_summary_columns:
+    tmp[col] = fSummary[col]
+
+for col in added_sequence_columns:
+    tmp[col] = fSequence[col]
+
+for col in added_junction_columns:
+    tmp[col] = fJunction[col]
+
+outFrame = tmp
+
+outFrame.columns = [u'ID', u'VDJ Frame', u'Top V Gene', u'Top D Gene', u'Top J Gene', u'CDR1 Seq', u'CDR1 Length', u'CDR2 Seq', u'CDR2 Length', u'CDR3 Seq', u'CDR3 Length', u'CDR3 Seq DNA', u'CDR3 Length DNA', u'Strand', u'CDR3 Found How', u'Functionality', 'V-REGION identity %', 'V-REGION identity nt', 'D-REGION reading frame', 'AA JUNCTION', 'Functionality comment', 'Sequence', 'FR1-IMGT', 'FR2-IMGT', 'FR3-IMGT', 'CDR3-IMGT', 'JUNCTION', 'J-REGION', 'FR4-IMGT', 'P3V-nt nb', 'N1-REGION-nt nb', 'P5D-nt nb', 'P3D-nt nb', 'N2-REGION-nt nb', 'P5J-nt nb', '3V-REGION trimmed-nt nb', '5D-REGION trimmed-nt nb', '3D-REGION trimmed-nt nb', '5J-REGION trimmed-nt nb']
+
+"""
+IGHV[0-9]-[0-9ab]+-?[0-9]?D?
+TRBV[0-9]{1,2}-?[0-9]?-?[123]?
+IGKV[0-3]D?-[0-9]{1,2}
+IGLV[0-9]-[0-9]{1,2}
+TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?
+TRGV[234589]
+TRDV[1-3]
+
+IGHD[0-9]-[0-9ab]+
+TRBD[12]
+TRDD[1-3]
+
+IGHJ[1-6]
+TRBJ[12]-[1-7]
+IGKJ[1-5]
+IGLJ[12367]
+TRAJ[0-9]{1,2}
+TRGJP?[12]
+TRDJ[1-4]
+"""
+
+vPattern = [r"(IGHV[0-9]-[0-9ab]+-?[0-9]?D?)",
+						r"(TRBV[0-9]{1,2}-?[0-9]?-?[123]?)",
+						r"(IGKV[0-3]D?-[0-9]{1,2})",
+						r"(IGLV[0-9]-[0-9]{1,2})",
+						r"(TRAV[0-9]{1,2}(-[1-46])?(/DV[45678])?)",
+						r"(TRGV[234589])",
+						r"(TRDV[1-3])"]
+
+dPattern = [r"(IGHD[0-9]-[0-9ab]+)",
+						r"(TRBD[12])",
+						r"(TRDD[1-3])"]
+						
+jPattern = [r"(IGHJ[1-6])",
+						r"(TRBJ[12]-[1-7])",
+						r"(IGKJ[1-5])",
+						r"(IGLJ[12367])",
+						r"(TRAJ[0-9]{1,2})",
+						r"(TRGJP?[12])",
+						r"(TRDJ[1-4])"]
+
+vPattern = re.compile(r"|".join(vPattern))
+
+dPattern = re.compile(r"|".join(dPattern))
+
+jPattern = re.compile(r"|".join(jPattern))
+
+
+def filterGenes(s, pattern):
+    if type(s) is not str:
+        return "NA"
+    res = pattern.search(s)
+    if res:
+        return res.group(0)
+    return "NA"
+
+
+
+outFrame["Top V Gene"] = outFrame["Top V Gene"].apply(lambda x: filterGenes(x, vPattern))
+outFrame["Top D Gene"] = outFrame["Top D Gene"].apply(lambda x: filterGenes(x, dPattern))
+outFrame["Top J Gene"] = outFrame["Top J Gene"].apply(lambda x: filterGenes(x, jPattern))
+
+
+tmp = outFrame["VDJ Frame"]
+tmp = tmp.replace("in-frame", "In-frame")
+tmp = tmp.replace("null", "Out-of-frame")
+tmp = tmp.replace("out-of-frame", "Out-of-frame")
+outFrame["VDJ Frame"] = tmp
+outFrame["CDR3 Length DNA"] = outFrame["CDR3 Seq DNA"].map(str).map(len)
+safeLength = lambda x: len(x) if type(x) == str else 0
+#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows?
+#outFrame = outFrame[(outFrame["CDR3 Seq DNA"].map(safeLength) > 0) & (outFrame["Top V Gene"] != "NA") & (outFrame["Top D Gene"] != "NA") & (outFrame["Top J Gene"] != "NA")] #filter out weird rows?
+outFrame.to_csv(outFile, sep="\t", index=False, index_label="index")

--- a/mutation_analysis.r	Thu Nov 05 10:16:32 2015 -0500
+++ b/mutation_analysis.r	Mon Nov 09 07:21:09 2015 -0500
@@ -358,7 +358,7 @@
 dat$percentage_mutations = round(dat$VRegionMutations / dat$VRegionNucleotides * 100, 2)
 
 p = ggplot(dat, aes(best_match, percentage_mutations))
-p = p + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA) + geom_point(aes(colour=best_match), position="jitter")
+p = p + geom_point(aes(colour=best_match), position="jitter") + geom_boxplot(aes(middle=mean(percentage_mutations)), alpha=0.1, outlier.shape = NA)
 p = p + xlab("Subclass") + ylab("Frequency") + ggtitle("Frequency scatter plot")
 write.table(dat[,c("Sequence.ID", "best_match", "VRegionMutations", "VRegionNucleotides", "percentage_mutations")], "scatter.txt", sep="\t",quote=F,row.names=F,col.names=T)
 

--- a/mutation_analysis.xml	Thu Nov 05 10:16:32 2015 -0500
+++ b/mutation_analysis.xml	Mon Nov 09 07:21:09 2015 -0500
@@ -1,7 +1,7 @@
 <tool id="mutation_analysis_shm" name="Mutation Analysis" version="1.0">
 	<description></description>
 	<command interpreter="bash">
-		wrapper.sh $in_file $method $out_file $out_file.files_path ${in_file.name} ${include_fr1} $functionality $unique
+		wrapper.sh $in_file $method $out_file $out_file.files_path ${in_file.name} ${include_fr1} $functionality $unique $naive_output
 	</command>
 	<inputs>
 		<param name="in_file" type="data" label="IMGT zip file to be analysed" />
@@ -27,10 +27,18 @@
 			<option value="AA.JUNCTION">AA.JUNCTION</option>
 			<option value="none">Don't remove duplicates</option>
 		</param>
-		
+		<conditional name="naive_output_cond">
+			<param name="naive_output" type="select" label="Output a file for naive analysis?">
+				<option value="yes">Yes</option>
+				<option value="no" selected="true">No</option>
+			</param>
+		</conditional>
 	</inputs>
 	<outputs>
 		<data format="html" name="out_file" label = "Mutation analysis on ${in_file.name}"/>
+		<data format="tabular" name="naive_output" label = "Naive input data from ${in_file.name}" >
+		    <filter>naive_output_cond['naive_output'] == "yes"</filter>
+		</data>
 	</outputs>
 	<help>
 		Takes an IMGT zip (http://www.imgt.org/HighV-QUEST/search.action) file and creates a summarization of the mutation analysis.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/naive_output.r	Mon Nov 09 07:21:09 2015 -0500
@@ -0,0 +1,23 @@
+args <- commandArgs(trailingOnly = TRUE)
+
+naive.file = args[1]
+shm.file = args[2]
+output.file = args[3]
+
+naive = read.table(naive.file, sep="\t", header=T, quote="", fill=T)
+shm.merge = read.table(shm.file, sep="\t", header=T, quote="", fill=T)
+
+
+final = merge(naive, shm.merge[,c("Sequence.ID", "best_match")], by.x="ID", by.y="Sequence.ID")
+print(paste("nrow final:", nrow(final)))
+names(final)[names(final) == "best_match"] = "Sample"
+final.numeric = final[,sapply(final, is.numeric)]
+final.numeric[is.na(final.numeric)] = 0
+final[,sapply(final, is.numeric)] = final.numeric
+print(paste("nrow final:", nrow(final)))
+final2 = final
+final2$Sample = gsub("[0-9]", "", final2$Sample)
+final = rbind(final, final2)
+final$Replicate = 1
+
+write.table(final, output.file, quote=F, sep="\t", row.names=F, col.names=T)

--- a/wrapper.sh	Thu Nov 05 10:16:32 2015 -0500
+++ b/wrapper.sh	Mon Nov 09 07:21:09 2015 -0500
@@ -9,6 +9,7 @@
 include_fr1=$6
 functionality=$7
 unique=$8
+naive_output=$9
 mkdir $outdir
 
 type="`file $input`"
@@ -25,6 +26,7 @@
 
 cat $PWD/files/*/1_* > $PWD/summary.txt
 cat $PWD/files/*/3_* > $PWD/sequences.txt
+cat $PWD/files/*/5_* > $PWD/aa.txt
 cat $PWD/files/*/6_* > $PWD/junction.txt
 cat $PWD/files/*/7_* > $PWD/mutationanalysis.txt
 cat $PWD/files/*/8_* > $PWD/mutationstats.txt
@@ -134,5 +136,25 @@
 
 echo "</html>" >> $output
 
+
+#optional output for naive
+
+if [[ "$naive_output" != "None" ]]
+then
+	echo "naive_output: $naive_output"
+	python $dir/imgt_loader.py --summ $PWD/summary.txt --aa $PWD/aa.txt --junction $PWD/junction.txt --output $naive_output
+	Rscript $dir/naive_output.r $naive_output $outdir/merged.txt $naive_output 2>&1
+fi
+
+
+
+
+
+
+
+
+
+
+
 #rm $outdir/HS12RSS.txt
 #rm $outdir/HS23RSS.txt