Mercurial > repos > dereeper > ragoo

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

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