Mercurial > repos > davidvanzessen > mutation_analysis

changeset 39:7377bf7e632d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 02 Nov 2015 04:52:26 -0500
parents 6b6dbbcc771d
children e022c21f8c47
files merge_and_filter.r mutation_analysis.py wrapper.sh
diffstat 3 files changed, 129 insertions(+), 114 deletions(-) [+]
line wrap: on
line diff
--- a/merge_and_filter.r	Tue Aug 04 08:55:14 2015 -0400
+++ b/merge_and_filter.r	Mon Nov 02 04:52:26 2015 -0500
@@ -2,17 +2,19 @@
 
 
 summaryfile = args[1]
-mutationanalysisfile = args[2]
-mutationstatsfile = args[3]
-hotspotsfile = args[4]
-gene_identification_file= args[5]
-output = args[6]
-unmatchedfile = args[7]
-method=args[8]
-functionality=args[9]
-unique_type=args[10]
+junctionfile = args[2]
+mutationanalysisfile = args[3]
+mutationstatsfile = args[4]
+hotspotsfile = args[5]
+gene_identification_file= args[6]
+output = args[7]
+unmatchedfile = args[8]
+method=args[9]
+functionality=args[10]
+unique_type=args[11]
 
 summ = read.table(summaryfile, header=T, sep="\t", fill=T, stringsAsFactors=F)
+junctions = read.table(junctionfile, header=T, sep="\t", fill=T, stringsAsFactors=F)
 mutationanalysis = read.table(mutationanalysisfile, header=T, sep="\t", fill=T, stringsAsFactors=F)
 mutationstats = read.table(mutationstatsfile, header=T, sep="\t", fill=T, stringsAsFactors=F)
 hotspots = read.table(hotspotsfile, header=T, sep="\t", fill=T, stringsAsFactors=F)
@@ -95,5 +97,11 @@
 print(paste("Number of rows in result:", nrow(result)))
 print(paste("Number of rows in unmatched:", nrow(unmatched)))
 
+
+#remove the sequences that have an 'n' (or 'N') in the junction.
+junctions = junctions[grepl("n|N", junctions$JUNCTION),]
+
+result = result[!(result$Sequence.ID %in% junctions$Sequence.ID),]
+
 write.table(x=result, file=output, sep="\t",quote=F,row.names=F,col.names=T)
 write.table(x=unmatched, file=unmatchedfile, sep="\t",quote=F,row.names=F,col.names=T)
--- a/mutation_analysis.py	Tue Aug 04 08:55:14 2015 -0400
+++ b/mutation_analysis.py	Mon Nov 02 04:52:26 2015 -0500
@@ -2,7 +2,8 @@
 import argparse
 
 parser = argparse.ArgumentParser()
-parser.add_argument("--input", help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation")
+parser.add_argument("--input",
+                    help="The '7_V-REGION-mutation-and-AA-change-table' and '10_V-REGION-mutation-hotspots' merged together, with an added 'best_match' annotation")
 parser.add_argument("--genes", help="The genes available in the 'best_match' column")
 parser.add_argument("--includefr1", help="The genes available in the 'best_match' column")
 parser.add_argument("--output", help="Output file")
@@ -13,7 +14,7 @@
 genes = str(args.genes).split(",")
 print "includefr1 =", args.includefr1
 include_fr1 = True if args.includefr1 == "yes" else False
-outfile = args.output 
+outfile = args.output
 
 genedic = dict()
 
@@ -28,58 +29,61 @@
 fr2Index = 0
 cdr2Index = 0
 fr3Index = 0
-first=True
+first = True
 IDlist = []
 mutationList = []
 
 with open(infile, 'r') as i:
-  for line in i:
-		if first:
-			linesplt = line.split("\t")
-			IDIndex = linesplt.index("Sequence.ID")
-			best_matchIndex = linesplt.index("best_match")
-			fr1Index = linesplt.index("FR1.IMGT")
-			cdr1Index = linesplt.index("CDR1.IMGT")
-			fr2Index = linesplt.index("FR2.IMGT")
-			cdr2Index = linesplt.index("CDR2.IMGT")
-			fr3Index = linesplt.index("FR3.IMGT")
-			first = False
-			continue
-		linecount += 1
-		linesplt = line.split("\t")
-		ID = linesplt[IDIndex]
-		genedic[ID] = linesplt[best_matchIndex]
-		mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if x] if include_fr1 else []
-		mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x]
-		mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x]
-		mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x]
-		mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"]
-		mutationdic[ID + "_FR3"] = [mutationMatcher.match(x).groups() for x in linesplt[fr3Index].split("|") if x]
-		
-		mutationList += mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
-		
-		IDlist += [ID]
+    for line in i:
+        if first:
+            linesplt = line.split("\t")
+            IDIndex = linesplt.index("Sequence.ID")
+            best_matchIndex = linesplt.index("best_match")
+            fr1Index = linesplt.index("FR1.IMGT")
+            cdr1Index = linesplt.index("CDR1.IMGT")
+            fr2Index = linesplt.index("FR2.IMGT")
+            cdr2Index = linesplt.index("CDR2.IMGT")
+            fr3Index = linesplt.index("FR3.IMGT")
+            first = False
+            continue
+        linecount += 1
+        linesplt = line.split("\t")
+        ID = linesplt[IDIndex]
+        genedic[ID] = linesplt[best_matchIndex]
+        mutationdic[ID + "_FR1"] = [mutationMatcher.match(x).groups() for x in linesplt[fr1Index].split("|") if
+                                    x] if include_fr1 else []
+        mutationdic[ID + "_CDR1"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr1Index].split("|") if x]
+        mutationdic[ID + "_FR2"] = [mutationMatcher.match(x).groups() for x in linesplt[fr2Index].split("|") if x]
+        mutationdic[ID + "_CDR2"] = [mutationMatcher.match(x).groups() for x in linesplt[cdr2Index].split("|") if x]
+        mutationdic[ID + "_FR2-CDR2"] = mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"]
+        mutationdic[ID + "_FR3"] = [mutationMatcher.match(x).groups() for x in linesplt[fr3Index].split("|") if x]
 
+        mutationList += mutationdic[ID + "_FR1"] + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[
+            ID + "_CDR2"] + mutationdic[ID + "_FR3"]
 
-AA_mutation = [0] * (int(max(mutationList, key=lambda i:int(i[4]) if i[4] else 0)[4]) + 1) #[4] is the position of the AA mutation, None if silent
+        IDlist += [ID]
+
+AA_mutation = [0] * (int(max(mutationList, key=lambda i: int(i[4]) if i[4] else 0)[
+                             4]) + 1)  # [4] is the position of the AA mutation, None if silent
 
 for mutation in mutationList:
-	if mutation[4]: #if non silent mutation
-		AA_mutation[int(mutation[4])] += 1
+    if mutation[4]:  # if non silent mutation
+        AA_mutation[int(mutation[4])] += 1
 
 aa_mutations_file = outfile[:outfile.rindex("/")] + "/aa_mutations.txt"
 with open(aa_mutations_file, 'w') as o:
-	o.write(",".join([str(x) for x in AA_mutation]) + "\n")
-    
+    o.write(",".join([str(x) for x in AA_mutation]) + "\n")
+
 if linecount == 0:
-	print "No data, exiting"
-	with open(outfile, 'w') as o:
-		o.write("RGYW (%)," + ("0,0,0\n" * len(genes)))
-		o.write("WRCY (%)," + ("0,0,0\n" * len(genes)))
-		o.write("WA (%)," + ("0,0,0\n" * len(genes)))
-		o.write("TW (%)," + ("0,0,0\n" * len(genes)))
-	import sys
-	sys.exit()
+    print "No data, exiting"
+    with open(outfile, 'w') as o:
+        o.write("RGYW (%)," + ("0,0,0\n" * len(genes)))
+        o.write("WRCY (%)," + ("0,0,0\n" * len(genes)))
+        o.write("WA (%)," + ("0,0,0\n" * len(genes)))
+        o.write("TW (%)," + ("0,0,0\n" * len(genes)))
+    import sys
+
+    sys.exit()
 
 hotspotMatcher = re.compile("[actg]+,(\d+)-(\d+)\((.*)\)")
 RGYWCount = {g: 0 for g in genes}
@@ -94,78 +98,80 @@
 atagcctIndex = 0
 first = True
 with open(infile, 'r') as i:
-	for line in i:
-		if first:
-			linesplt = line.split("\t")
-			ataIndex = linesplt.index("X.a.t.a")
-			tatIndex = linesplt.index("t.a.t.")
-			aggctatIndex = linesplt.index("X.a.g.g.c.t..a.t.")
-			atagcctIndex = linesplt.index("X.a.t..a.g.c.c.t.")
-			first = False
-			continue
-		linesplt = line.split("\t")
-		gene = linesplt[best_matchIndex]
-		ID = linesplt[IDIndex]
-		RGYW = [(int(x),int(y),z) for (x,y,z) in [hotspotMatcher.match(x).groups() for x in linesplt[aggctatIndex].split("|") if x]]
-		WRCY = [(int(x),int(y),z) for (x,y,z) in [hotspotMatcher.match(x).groups() for x in linesplt[atagcctIndex].split("|") if x]]
-		WA = [(int(x),int(y),z) for (x,y,z) in [hotspotMatcher.match(x).groups() for x in linesplt[ataIndex].split("|") if x]]
-		TW = [(int(x),int(y),z) for (x,y,z) in [hotspotMatcher.match(x).groups() for x in linesplt[tatIndex].split("|") if x]]
-		RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0,0,0,0
-		
-		mutationList = (mutationdic[ID + "_FR1"] if include_fr1 else []) + mutationdic[ID + "_CDR1"] + mutationdic[ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
-		for mutation in mutationList:
-			frm, where, to, AAfrm, AAwhere, AAto, junk = mutation
-			mutation_in_RGYW = any([(start <= int(where) <= end) for (start,end,region) in RGYW])
-			mutation_in_WRCY = any([(start <= int(where) <= end) for (start,end,region) in WRCY])
-			mutation_in_WA = any([(start <= int(where) <= end) for (start,end,region) in WA])
-			mutation_in_TW = any([(start <= int(where) <= end) for (start,end,region) in TW])
-			
-			in_how_many_motifs = sum([mutation_in_RGYW, mutation_in_WRCY, mutation_in_WA, mutation_in_TW])
-			
-			if in_how_many_motifs > 0:
-				RGYWCount[ID] += (1.0 * int(mutation_in_RGYW)) / in_how_many_motifs
-				WRCYCount[ID] += (1.0 * int(mutation_in_WRCY)) / in_how_many_motifs
-				WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs
-				TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs
+    for line in i:
+        if first:
+            linesplt = line.split("\t")
+            ataIndex = linesplt.index("X.a.t.a")
+            tatIndex = linesplt.index("t.a.t.")
+            aggctatIndex = linesplt.index("X.a.g.g.c.t..a.t.")
+            atagcctIndex = linesplt.index("X.a.t..a.g.c.c.t.")
+            first = False
+            continue
+        linesplt = line.split("\t")
+        gene = linesplt[best_matchIndex]
+        ID = linesplt[IDIndex]
+        RGYW = [(int(x), int(y), z) for (x, y, z) in
+                [hotspotMatcher.match(x).groups() for x in linesplt[aggctatIndex].split("|") if x]]
+        WRCY = [(int(x), int(y), z) for (x, y, z) in
+                [hotspotMatcher.match(x).groups() for x in linesplt[atagcctIndex].split("|") if x]]
+        WA = [(int(x), int(y), z) for (x, y, z) in
+              [hotspotMatcher.match(x).groups() for x in linesplt[ataIndex].split("|") if x]]
+        TW = [(int(x), int(y), z) for (x, y, z) in
+              [hotspotMatcher.match(x).groups() for x in linesplt[tatIndex].split("|") if x]]
+        RGYWCount[ID], WRCYCount[ID], WACount[ID], TWCount[ID] = 0, 0, 0, 0
+
+        mutationList = (mutationdic[ID + "_FR1"] if include_fr1 else []) + mutationdic[ID + "_CDR1"] + mutationdic[
+            ID + "_FR2"] + mutationdic[ID + "_CDR2"] + mutationdic[ID + "_FR3"]
+        for mutation in mutationList:
+            frm, where, to, AAfrm, AAwhere, AAto, junk = mutation
+            mutation_in_RGYW = any([(start <= int(where) <= end) for (start, end, region) in RGYW])
+            mutation_in_WRCY = any([(start <= int(where) <= end) for (start, end, region) in WRCY])
+            mutation_in_WA = any([(start <= int(where) <= end) for (start, end, region) in WA])
+            mutation_in_TW = any([(start <= int(where) <= end) for (start, end, region) in TW])
+
+            in_how_many_motifs = sum([mutation_in_RGYW, mutation_in_WRCY, mutation_in_WA, mutation_in_TW])
+
+            if in_how_many_motifs > 0:
+                RGYWCount[ID] += (1.0 * int(mutation_in_RGYW)) / in_how_many_motifs
+                WRCYCount[ID] += (1.0 * int(mutation_in_WRCY)) / in_how_many_motifs
+                WACount[ID] += (1.0 * int(mutation_in_WA)) / in_how_many_motifs
+                TWCount[ID] += (1.0 * int(mutation_in_TW)) / in_how_many_motifs
 
 directory = outfile[:outfile.rfind("/") + 1]
 value = 0
 valuedic = dict()
 for gene in genes:
-	with open(directory + gene + "_value.txt", 'r') as v:
-		valuedic[gene] = int(v.readlines()[0].rstrip())
+    with open(directory + gene + "_value.txt", 'r') as v:
+        valuedic[gene] = int(v.readlines()[0].rstrip())
 with open(directory + "total_value.txt", 'r') as v:
-	valuedic["total"] = int(v.readlines()[0].rstrip())
+    valuedic["total"] = int(v.readlines()[0].rstrip())
 
 dic = {"RGYW": RGYWCount, "WRCY": WRCYCount, "WA": WACount, "TW": TWCount}
 arr = ["RGYW", "WRCY", "WA", "TW"]
-with open(outfile, 'w') as o:	
-	for typ in arr:
-		o.write(typ + " (%)")
-		curr = dic[typ]
-		for gene in genes:
-			geneMatcher = re.compile(".*" + gene + ".*")
-			if valuedic[gene] is 0:
-				o.write(",0,0,0")
-			else:
-				x = int(round(sum([curr[x] for x in [y for y,z in genedic.iteritems() if geneMatcher.match(z)]])))
-				y = valuedic[gene]
-				z = str(round(x / float(valuedic[gene]) * 100, 1))
-				o.write("," + str(x) + "," + str(y) + "," + z)
-		#for total
-		x = int(round(sum([y for x,y in curr.iteritems()])))
-		y = valuedic["total"]
-		z = str(round(x / float(valuedic["total"]) * 100, 1))
-		o.write("," + str(x) + "," + str(y) + "," + z + "\n")
+with open(outfile, 'w') as o:
+    for typ in arr:
+        o.write(typ + " (%)")
+        curr = dic[typ]
+        for gene in genes:
+            geneMatcher = re.compile(".*" + gene + ".*")
+            if valuedic[gene] is 0:
+                o.write(",0,0,0")
+            else:
+                x = int(round(sum([curr[x] for x in [y for y, z in genedic.iteritems() if geneMatcher.match(z)]])))
+                y = valuedic[gene]
+                z = str(round(x / float(valuedic[gene]) * 100, 1))
+                o.write("," + str(x) + "," + str(y) + "," + z)
+        # for total
+        x = int(round(sum([y for x, y in curr.iteritems()])))
+        y = valuedic["total"]
+        z = str(round(x / float(valuedic["total"]) * 100, 1))
+        o.write("," + str(x) + "," + str(y) + "," + z + "\n")
 
 
-#for testing
+# for testing
 seq_motif_file = outfile[:outfile.rindex("/")] + "/motif_per_seq.txt"
-first = True
 with open(seq_motif_file, 'w') as o:
-	for ID in IDlist:
-		if first:
-			o.write("ID\tRGYWC\tWRCY\tWA\tTW\n")
-			first = False
-			continue
-		o.write(ID + "\t" + str(round(RGYWCount[ID], 2)) + "\t" + str(round(WRCYCount[ID], 2)) + "\t" + str(round(WACount[ID], 2)) + "\t" + str(round(TWCount[ID], 2)) + "\n")
+    o.write("ID\tRGYWC\tWRCY\tWA\tTW\n")
+    for ID in IDlist:
+        o.write(ID + "\t" + str(round(RGYWCount[ID], 2)) + "\t" + str(round(WRCYCount[ID], 2)) + "\t" + str(
+            round(WACount[ID], 2)) + "\t" + str(round(TWCount[ID], 2)) + "\n")
--- a/wrapper.sh	Tue Aug 04 08:55:14 2015 -0400
+++ b/wrapper.sh	Mon Nov 02 04:52:26 2015 -0500
@@ -24,6 +24,7 @@
 fi
 
 cat $PWD/files/*/1_* > $PWD/summary.txt
+cat $PWD/files/*/6_* > $PWD/junction.txt
 cat $PWD/files/*/7_* > $PWD/mutationanalysis.txt
 cat $PWD/files/*/8_* > $PWD/mutationstats.txt
 cat $PWD/files/*/10_* > $PWD/hotspots.txt
@@ -53,7 +54,7 @@
 
 
 echo "merging"
-Rscript $dir/merge_and_filter.r $PWD/summary.txt $PWD/mutationanalysis.txt $PWD/mutationstats.txt $PWD/hotspots.txt $outdir/identified_genes.txt $outdir/merged.txt $outdir/unmatched.txt $method $functionality $unique
+Rscript $dir/merge_and_filter.r $PWD/summary.txt $PWD/junction.txt $PWD/mutationanalysis.txt $PWD/mutationstats.txt $PWD/hotspots.txt $outdir/identified_genes.txt $outdir/merged.txt $outdir/unmatched.txt $method $functionality $unique
 
 genes="ca,ca1,ca2,cg,cg1,cg2,cg3,cg4,cm"
 echo "R mutation analysis"