Mercurial > repos > dereeper > ragoo
comparison RaGOO/make_agp.py @ 13:b9a3aeb162ab draft default tip
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") |