added circular sequence fetching script for newer spades outputs

Author: rozovr

get_simple_cycs.py

@@ -62,17 +62,20 @@ def parse_user_input():
 cycs = set([])
 fp = open(fastg_name, 'r')
 
-simple_cycs_ofile = fp.name.replace(".fastg", ".simple_cycs.fasta")
+simple_cycs_ofile = fp.name.replace(".fastg", ".simple.cycs.fasta")
 f_cycs_out = open(simple_cycs_ofile, 'w')
-filt_graph_ofile = fp.name.replace(".fastg", ".simple_cycs_filtered.fastg")
+filt_graph_ofile = fp.name.replace(".fastg", ".simple.cycs_filtered.fastg")
 f_graph_out = open(filt_graph_ofile, 'w')
 
 
 #### potential issue: output fastg 
 #### may contain links to nodes that have been removed
 #### may be potential hazard downstream
-
+cnt = 0
 for name,seq,qual in readfq(fp):
+    cnt+=1
+    if cnt%2!=1:
+        continue
     orig_name = name
     name = re.sub('[:,]'," ", name[:-1])
     name_parts = name.split(" ")

paper/get_spades_cyc_rr_contigs.py

@@ -0,0 +1,91 @@
+import argparse, sys, re, os
+sys.path.insert(0, '../recycle/')
+from recycle.utils import *
+import numpy as np
+
+
+def parse_user_input():
+    parser = argparse.ArgumentParser(description='gets cycles from spades contigs output. Outputs these to separate fasta file')
+    parser.add_argument('-g','--graph', help='Input (SPAdes 3.60+) FASTG to process', required=True, type=str)
+    parser.add_argument('-p','--paths', help='Input (SPAdes 3.60+) contigs.paths file to process', required=True, type=str)
+    parser.add_argument('-s','--sequences', help='Input (SPAdes 3.60+) contigs.fasta sequences to process', required=True, type=str)
+    parser.add_argument('-m','--length', help='Minimum cycle length to keep (shorter cycles put in new graph file; default = 1000)', 
+        required=False, type=int, default=1000)
+    parser.add_argument('-k','--max_k', help='integer reflecting maximum k value used by the assembler',
+     required=True, type=int, default=55)
+    return parser.parse_args()
+
+def get_cyc_labels_from_paths_file(paths):
+	""" scans lines of paths file
+		when length greater than 2 and 
+		first node same as last, add label of path
+		to returned set; add only forward strand versions
+	"""
+	labels = set({})
+	fp = open(paths, 'r')
+	lines = fp.read()
+	fp.close()
+	lines = re.sub(';\n',",", lines)
+	lines = lines.split("\n")
+	last_label = ''
+	for ind in range(len(lines)):
+		line = lines[ind].rstrip()
+		if (ind % 4 == 2 or ind % 4 == 3):
+			continue
+		if ind%4==0:
+			last_label = line
+		else:
+			path = line.split(",")
+			if len(path)>1 and path[0]==path[-1]:
+				labels.add(last_label)
+	return labels	
+
+def update_paths_dict(paths_dict, labels, seqs):
+	fs = open(seqs, 'r')
+	for name,seq,qual in readfq(fs):
+		if name in labels:
+			paths_dict[name] = seq
+	return seqs
+
+
+
+args = parse_user_input()
+fastg = args.graph
+paths = args.paths
+seqs = args.sequences
+min_length = args.length
+max_k = args.max_k
+fg = open(fastg, 'r')
+files_dir = os.path.dirname(fg.name)
+
+G = get_fastg_digraph(fastg)
+SEQS = get_fastg_seqs_dict(fastg,G)
+long_self_loops = get_long_self_loops(G, min_length, SEQS)
+final_paths_dict = {}
+path_count = 0
+
+for nd in long_self_loops:
+    name = get_spades_type_name(path_count, nd,
+        SEQS, G, get_cov_from_spades_name(nd[0]))
+    final_paths_dict[name] = nd
+    path_count += 1
+
+path_labels_set = get_cyc_labels_from_paths_file(paths)
+update_paths_dict(final_paths_dict, path_labels_set, seqs)
+
+# output - fasta of cyc sequences
+fp = open(paths, 'r')
+(root,ext) = os.path.splitext(fp.name)
+fasta_ofile = root + ext.replace(".paths", ".cycs.fasta")
+f_cycs_fasta = open(fasta_ofile, 'w')
+
+for p in final_paths_dict.keys():
+    if p not in path_labels_set:
+    	# cycle = False b/c want last k-mer to match first
+    	seq = get_seq_from_path(final_paths_dict[p], SEQS, max_k_val=max_k, cycle=False)
+    else: 
+    	seq = final_paths_dict[p]
+    # 
+    # print ""
+    if len(seq)>=min_length:
+        f_cycs_fasta.write(">" + p + "\n" + seq + "\n")

paper/nucmer.sim.mk

@@ -1,14 +1,14 @@
 # creates a report on the number of plasmids reported
 # and their alignment relative to some reference
 REF_CNT = 100
-CV = 0.25
+CV = 0.5
 PARENT_DIR = /home/nasheran/rozovr/recycle_paper_data
 INPUT_DIR = $(PARENT_DIR)/ref_$(REF_CNT)
 # TWO_STEP_ASSEM = 0
 # REPEAT_RES_1 = 0
 # REPEAT_RES_2 = 0
 REF = $(PARENT_DIR)/plasmids_sim_ref_$(REF_CNT).cln.fasta
-INPUT = $(INPUT_DIR)/assembly_graph
+INPUT = $(INPUT_DIR)/contigs
 # ifeq ($(REPEAT_RES_1),0)
 # 	INPUT1 = before_rr
 # else
@@ -36,7 +36,7 @@ all:  $(INPUT).cycs.dbl.fasta $(INPUT).nucmer.delta $(INPUT).nucmer$(CV).summary
 
 # 1) concat output cycles to self:
 $(INPUT).cycs.dbl.fasta: $(INPUT).cycs.fasta
-	python ~/recycle/paper/concat_seq_to_self.py -i $(INPUT).cycs.fasta -m 0
+	python ~/recycle/paper/concat_seq_to_self.py -i $(INPUT).cycs.fasta -m 55
 
 # 2) run nucmer
 $(INPUT).nucmer.delta: $(INPUT).cycs.dbl.fasta

paper/tally_hits.py

@@ -1,7 +1,9 @@
+import argparse, sys
+sys.path.insert(0, '../recycle/')
 from recycle.utils import *
-import argparse
 import numpy as np
 
+
 def count_property_range_hits(prop, node_dict, hits):
 	""" picks which values to use in tuples based on property
 		counts vals having min_val <= val < max_val

recycle/utils.py

@@ -80,11 +80,9 @@ def get_fastg_seqs_dict(fastg_name, G):
     """ 
     fp = open(fastg_name, 'r')
     seqs = {}
-    # nodes = G.nodes()
     for name,seq,qual in readfq(fp):
         name_parts = re.sub('[:,]'," ", name[:-1]).split()
         node = name_parts[0]
-        # if node in nodes:
         seqs[node] = seq
     return seqs