annotate dante_gff_to_dna.py @ 18:039c45c01b47 draft

Uploaded
author petr-novak
date Tue, 24 Sep 2019 07:27:56 -0400
parents 1a766f9f623d
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
1 #!/usr/bin/env python3
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
2
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
3 import argparse
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
4 import time
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
5 import os
10
d0431a839606 Uploaded
petr-novak
parents: 0
diff changeset
6 import textwrap
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
7 from collections import defaultdict
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
8 from Bio import SeqIO
10
d0431a839606 Uploaded
petr-novak
parents: 0
diff changeset
9 import configuration
17
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
10 from dante_gff_output_filtering import parse_gff_line
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
11 t_nt_seqs_extraction = time.time()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
12
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
13
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
14 def str2bool(v):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
15 if v.lower() in ('yes', 'true', 't', 'y', '1'):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
16 return True
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
17 elif v.lower() in ('no', 'false', 'f', 'n', '0'):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
18 return False
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
19 else:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
20 raise argparse.ArgumentTypeError('Boolean value expected')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
21
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
22
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
23 def check_file_start(gff_file):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
24 count_comment = 0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
25 with open(gff_file, "r") as gff_all:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
26 line = gff_all.readline()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
27 while line.startswith("#"):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
28 line = gff_all.readline()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
29 count_comment += 1
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
30 return count_comment, line
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
31
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
32
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
33 def extract_nt_seqs(DNA_SEQ, DOM_GFF, OUT_DIR, CLASS_TBL, EXTENDED):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
34 ''' Extract nucleotide sequences of protein domains found by DANTE from input DNA seq.
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
35 Sequences are saved in fasta files separately for each transposon lineage.
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
36 Sequences extraction is based on position of Best_Hit alignment reported by LASTAL.
17
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
37 The positions can be extended (optional) based on what part of database domain was aligned
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
38 (Best_Hit_DB_Pos attribute).
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
39 The strand orientation needs to be considered in extending and extracting the sequence itself
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
40 '''
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
41 [count_comment, first_line] = check_file_start(DOM_GFF)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
42 unique_classes = get_unique_classes(CLASS_TBL)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
43 files_dict = defaultdict(str)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
44 domains_counts_dict = defaultdict(int)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
45 allSeqs = SeqIO.to_dict(SeqIO.parse(DNA_SEQ, 'fasta'))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
46 with open(DOM_GFF, "r") as domains:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
47 for comment_idx in range(count_comment):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
48 next(domains)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
49 seq_id_stored = first_line.split("\t")[0]
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
50 allSeqs = SeqIO.to_dict(SeqIO.parse(DNA_SEQ, 'fasta'))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
51 seq_nt = allSeqs[seq_id_stored]
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
52 for line in domains:
17
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
53 gff_line = parse_gff_line(line)
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
54 elem_type = gff_line['attributes']['Final_Classification']
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
55 if elem_type == configuration.AMBIGUOUS_TAG:
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
56 continue # skip ambiguous classification
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
57 seq_id = gff_line['seqid']
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
58 dom_type = gff_line['attributes']['Name']
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
59 strand = gff_line['strand']
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
60 align_nt_start = int(gff_line['attributes']['Best_Hit'].split(":")[
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
61 -1].split("-")[0])
17
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
62 align_nt_end = int(gff_line['attributes']['Best_Hit'].split(":")[
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
63 -1].split("-")[1].split("[")[0])
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
64 if seq_id != seq_id_stored:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
65 seq_id_stored = seq_id
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
66 seq_nt = allSeqs[seq_id_stored]
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
67 if EXTENDED:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
68 ## which part of database sequence was aligned
17
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
69 db_part = gff_line['attributes']['Best_Hit_DB_Pos']
1a766f9f623d Uploaded
petr-novak
parents: 10
diff changeset
70 ## db_part = line.split("\t")[8].split(";")[4].split("=")[1]
0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
71 ## datatabse seq length
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
72 dom_len = int(db_part.split("of")[1])
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
73 ## start of alignment on database seq
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
74 db_start = int(db_part.split("of")[0].split(":")[0])
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
75 ## end of alignment on database seq
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
76 db_end = int(db_part.split("of")[0].split(":")[1])
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
77 ## number of nucleotides missing in the beginning
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
78 dom_nt_prefix = (db_start - 1) * 3
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
79 ## number of nucleotides missing in the end
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
80 dom_nt_suffix = (dom_len - db_end) * 3
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
81 if strand == "+":
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
82 dom_nt_start = align_nt_start - dom_nt_prefix
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
83 dom_nt_end = align_nt_end + dom_nt_suffix
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
84 ## reverse extending for - strand
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
85 else:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
86 dom_nt_start = align_nt_start - dom_nt_suffix
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
87 dom_nt_end = align_nt_end + dom_nt_prefix
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
88 ## correction for domain after extending having negative starting positon
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
89 dom_nt_start = max(1, dom_nt_start)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
90 else:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
91 dom_nt_start = align_nt_start
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
92 dom_nt_end = align_nt_end
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
93 full_dom_nt = seq_nt.seq[dom_nt_start - 1:dom_nt_end]
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
94 ## for - strand take reverse complement of the extracted sequence
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
95 if strand == "-":
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
96 full_dom_nt = full_dom_nt.reverse_complement()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
97 full_dom_nt = str(full_dom_nt)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
98 ## report when domain classified to the last level and no Ns in extracted seq
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
99 if elem_type in unique_classes and "N" not in full_dom_nt:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
100 # lineages reported in separate fasta files
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
101 if not elem_type in files_dict:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
102 files_dict[elem_type] = os.path.join(
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
103 OUT_DIR, "{}.fasta".format(elem_type.split("|")[
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
104 -1].replace("/", "_")))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
105 with open(files_dict[elem_type], "a") as out_nt_seq:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
106 out_nt_seq.write(">{}:{}-{}|{}[{}]\n{}\n".format(
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
107 seq_nt.id, dom_nt_start, dom_nt_end, dom_type,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
108 elem_type, textwrap.fill(full_dom_nt,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
109 configuration.FASTA_LINE)))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
110 domains_counts_dict[elem_type] += 1
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
111 return domains_counts_dict
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
112
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
113
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
114 def get_unique_classes(CLASS_TBL):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
115 ''' Get all the lineages of current domains classification table to check if domains are classified to the last level.
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
116 Only the sequences of unambiguous and completely classified domains will be extracted.
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
117 '''
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
118 unique_classes = []
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
119 with open(CLASS_TBL, "r") as class_tbl:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
120 for line in class_tbl:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
121 line_class = "|".join(line.rstrip().split("\t")[1:])
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
122 if line_class not in unique_classes:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
123 unique_classes.append(line_class)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
124 return unique_classes
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
125
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
126
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
127 def write_domains_stat(domains_counts_dict, OUT_DIR):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
128 ''' Report counts of domains for individual lineages
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
129 '''
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
130 total = 0
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
131 with open(
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
132 os.path.join(OUT_DIR,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
133 configuration.EXTRACT_DOM_STAT), "w") as dom_stat:
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
134 for domain, count in domains_counts_dict.items():
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
135 dom_stat.write(";{}:{}\n".format(domain, count))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
136 total += count
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
137 dom_stat.write(";TOTAL:{}\n".format(total))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
138
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
139
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
140 def main(args):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
141
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
142 DNA_SEQ = args.input_dna
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
143 DOM_GFF = args.domains_gff
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
144 OUT_DIR = args.out_dir
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
145 CLASS_TBL = args.classification
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
146 EXTENDED = args.extended
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
147
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
148 if not os.path.exists(OUT_DIR):
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
149 os.makedirs(OUT_DIR)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
150
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
151 domains_counts_dict = extract_nt_seqs(DNA_SEQ, DOM_GFF, OUT_DIR, CLASS_TBL,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
152 EXTENDED)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
153 write_domains_stat(domains_counts_dict, OUT_DIR)
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
154
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
155 print("ELAPSED_TIME_EXTRACTION = {} s\n".format(time.time() -
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
156 t_nt_seqs_extraction))
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
157
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
158
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
159 if __name__ == "__main__":
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
160
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
161 # Command line arguments
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
162 parser = argparse.ArgumentParser()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
163 parser.add_argument('-i',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
164 '--input_dna',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
165 type=str,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
166 required=True,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
167 help='path to input DNA sequence')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
168 parser.add_argument('-d',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
169 '--domains_gff',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
170 type=str,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
171 required=True,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
172 help='GFF file of protein domains')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
173 parser.add_argument('-cs',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
174 '--classification',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
175 type=str,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
176 required=True,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
177 help='protein domains classification file')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
178 parser.add_argument('-out',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
179 '--out_dir',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
180 type=str,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
181 default=configuration.EXTRACT_OUT_DIR,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
182 help='output directory')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
183 parser.add_argument(
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
184 '-ex',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
185 '--extended',
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
186 type=str2bool,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
187 default=True,
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
188 help=
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
189 'extend the domains edges if not the whole datatabase sequence was aligned')
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
190 args = parser.parse_args()
77d9f2ecb28a Uploaded
petr-novak
parents:
diff changeset
191 main(args)