annotate RaGOO/make_agp.py @ 13:b9a3aeb162ab draft default tip

Uploaded
author dereeper
date Mon, 26 Jul 2021 18:22:37 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
13
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
1 #!/usr/bin/env python
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
2 import sys
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
3 import argparse
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
4
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
5 if __name__ == "__main__":
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
6 parser = argparse.ArgumentParser(description='Produce an AGP v2.0 file describing the ordering and orienting performed by RaGOO.')
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
7 parser.add_argument("orderings", metavar="<orderings.fofn>", type=str, help="file of orderings files ($ls ragoo_output/orderings/* > orderings.fofn)")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
8 parser.add_argument("fai", metavar="<contigs.fasta.fai>", type=str, help="index file for contigs ($samtools faidx contigs.fasta)")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
9 parser.add_argument("gap_len", metavar="gap_size", type=int, help="Gap size used for pseudomolecule padding.")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
10
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
11 # Get the command line arguments
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
12 args = parser.parse_args()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
13 orderings_fofn = args.orderings
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
14 fai_file = args.fai
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
15 gap_len = args.gap_len
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
16
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
17 # Save the contig orderings
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
18 orderings = dict()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
19 with open(orderings_fofn, 'r') as f1:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
20 for line1 in f1:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
21 chrom = line1[line1.rfind("/")+1:line1.rfind("_orderings.txt")] + "_RaGOO"
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
22 orderings[chrom] = []
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
23 with open(line1.rstrip(), 'r') as f2:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
24 for line2 in f2:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
25 L1 = line2.rstrip().split('\t')
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
26 # Append (ctg_header, orientation)
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
27 orderings[chrom].append((L1[0], L1[1]))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
28
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
29 # Get contig lengths
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
30 ctg_lens = dict()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
31 with open(fai_file, 'r') as f:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
32 for line in f:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
33 L1 = line.split('\t')
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
34 ctg_lens[L1[0]] = int(L1[1])
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
35
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
36 # Write the final AGP file
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
37 current_pos = 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
38 idx = 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
39 chrom_idx = 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
40 out_buff = []
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
41
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
42 sys.stdout.write("## AGP-version 2.0\n")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
43 sys.stdout.write("## AGP constructed by RaGOO\n")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
44 for chrom in orderings:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
45 if orderings[chrom]:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
46 for ctg, strand in orderings[chrom]:
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
47 start = current_pos
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
48 end = current_pos + ctg_lens[ctg] - 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
49 line_buff = list()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
50
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
51 # Write the line for the contig
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
52 line_buff.append(chrom)
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
53 line_buff.append(str(start))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
54 line_buff.append(str(end))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
55 line_buff.append(str(idx))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
56 line_buff.append("W")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
57 line_buff.append(ctg)
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
58 line_buff.append("1")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
59 line_buff.append(str(end - start + 1))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
60 line_buff.append(strand)
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
61 out_buff.append("\t".join(line_buff))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
62
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
63 # Write the line for the gap
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
64 idx += 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
65 current_pos += ctg_lens[ctg]
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
66 start = current_pos
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
67 end = current_pos + gap_len - 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
68 line_buff = list()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
69
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
70 line_buff.append(chrom)
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
71 line_buff.append(str(start))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
72 line_buff.append(str(end))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
73 line_buff.append(str(idx))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
74 line_buff.append("N")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
75 line_buff.append(str(gap_len))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
76 line_buff.append("scaffold")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
77 line_buff.append("yes")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
78 line_buff.append("align_genus")
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
79 out_buff.append("\t".join(line_buff))
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
80
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
81 idx += 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
82 current_pos += gap_len
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
83
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
84 # Pop the last gap since we don't end with padding
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
85 out_buff.pop()
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
86 idx =1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
87 chrom_idx += 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
88 current_pos = 1
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
89
b9a3aeb162ab Uploaded
dereeper
parents:
diff changeset
90 sys.stdout.write("\n".join(out_buff) + "\n")