-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathalignAccessoryGenomes
executable file
·178 lines (137 loc) · 6.35 KB
/
alignAccessoryGenomes
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#!/usr/bin/env python3
# encoding: utf-8
'''
alignAccessoryGenomes -- Run nucmer to find alignments of all accessory genome contigs against themselves to find potential regions of non-homologous recombination
alignAccessoryGenomes
@author: Darrin Lemmer
@copyright: 2018 TGen North. All rights reserved.
@license: ACADEMIC AND RESEARCH LICENSE -- see ../LICENSE
@contact: [email protected]
'''
import subprocess
import sys
import re
import os
import argparse
import skbio
__version__ = '0.1.0'
__date__ = '2018-10-30'
__updated__ = '2018-10-30'
def expand_path(path):
user_match = re.match('^(~)(.*)$', path)
if user_match:
path = os.path.expanduser(path)
return os.path.abspath(path)
def _parse_delta_line(current_line, current_contigs, output_handle):
line_match = re.match(r'^>(.*)_NODE[^ ]+ (.*)_NODE[^ ]+ (\d+) (\d+)\s*$', current_line)
if line_match:
current_contigs = ( line_match.group(1), line_match.group(2) )
if (current_contigs is None or current_contigs[0] != current_contigs[1]):
output_handle.write(current_line)
return current_contigs
def parse_delta_file(delta_filename, output_filename):
current_contigs = None
delta_handle = open(delta_filename, 'r')
output_handle = open(output_filename, 'w')
for current_line in delta_handle:
current_contigs = _parse_delta_line(current_line, current_contigs, output_handle)
output_handle.close()
delta_handle.close()
def remove_unfit_contigs(fasta, min_length, min_coverage):
from skbio import DNA
for seq in skbio.io.registry.read(fasta, format='fasta', constructor=DNA):
contig_match = re.match(r'^(.*)_NODE_\d+_length_(\d+)_cov_(\d+)\.[^ ]+$', seq.metadata['id'])
if contig_match:
name = contig_match.group(1)
length = int(contig_match.group(2))
coverage = int(contig_match.group(3))
if length >= min_length and coverage >= min_coverage:
yield seq
else:
yield seq
def run_nucmer(fasta, breaklength):
output = subprocess.getoutput("nucmer -b %i -p pairwise %s %s" % (breaklength, fasta, fasta))
return output
class CLIError(Exception):
'''Generic exception to raise and log different fatal errors.'''
def __init__(self, msg):
super(CLIError).__init__(type(self))
self.msg = "E: %s" % msg
def __str__(self):
return self.msg
def __unicode__(self):
return self.msg
def main(argv=None):
'''Command line options.'''
if argv is None:
argv = sys.argv
else:
sys.argv.extend(argv)
program_name = os.path.basename(sys.argv[0])
program_version = "v%s" % __version__
program_build_date = str(__updated__)
program_version_message = '%%(prog)s %s (%s)' % (program_version, program_build_date)
if __name__ == '__main__':
program_shortdesc = __import__('__main__').__doc__.split("\n")[1]
else:
program_shortdesc = __doc__.split("\n")[1]
#program_shortdesc = __import__('__main__').__doc__.split("\n")[1]
program_license = '''%s
Created by TGen North on %s.
Copyright 2018 TGen North. All rights reserved.
Available for academic and research use only under a license
from The Translational Genomics Research Institute (TGen)
that is free for non-commercial use.
Distributed on an "AS IS" basis without warranties
or conditions of any kind, either express or implied.
USAGE
''' % (program_shortdesc, str(__date__))
try:
# Setup argument parser
parser = argparse.ArgumentParser(description=program_license, formatter_class=argparse.RawDescriptionHelpFormatter)
required_group = parser.add_argument_group("required arguments")
optional_group = parser.add_argument_group("optional arguments")
required_group.add_argument("-f", "--fasta", required=True, help="path to concatenation of all accessory genome contigs. [REQUIRED]")
#required_group.add_argument("-b", "--bam-dir", required=True, dest="bdir", metavar="DIR", help="directory of bam files to analyze. [REQUIRED]")
optional_group.add_argument("-o", "--out-dir", dest="odir", metavar="DIR", help="directory to write output files to. [default: `pwd`]")
optional_group.add_argument("-l", "--min-length", dest="minlen", default=5000, type=int, help="minimum length alignment to keep. [default: 5000]")
optional_group.add_argument("-c", "--min-coverage", dest="mincov", default=2, type=int, help="minimum average coverage required. [default: 2]")
#optional_group.add_argument("-s", "--submitter", dest="job_manager", default="SLURM", help="cluster job submitter to use (PBS, SLURM, SGE, none). [default: SLURM]")
#optional_group.add_argument("--submitter-args", dest="sargs", metavar="ARGS", help="additional arguments to pass to the job submitter, enclosed in \"\".")
optional_group.add_argument("-V", "--version", action="version", version=program_version_message)
# Process arguments
args = parser.parse_args()
#bam_dir = args.bdir
out_dir = args.odir
min_length = args.minlen
min_coverage = args.mincov
fasta = expand_path(args.fasta)
if not out_dir:
out_dir = os.getcwd()
out_dir = expand_path(out_dir)
#bam_dir = expand_path(bam_dir)
#if os.path.exists(out_dir):
# response = input(
# "\nOutput folder %s already exists!\nFiles in it may be overwritten!\nShould we continue anyway [N]? " % out_dir)
# if not re.match('^[Yy]', response):
# print("Operation cancelled!")
# quit()
#else:
if not os.path.exists(out_dir):
os.makedirs(out_dir)
except KeyboardInterrupt:
### handle keyboard interrupt ###
return 0
except Exception as e:
indent = len(program_name) * " "
sys.stderr.write(program_name + ": " + repr(e) + "\n")
sys.stderr.write(indent + " for help use --help\n")
return 2
trimmed_fasta = os.path.join(out_dir, "accessory_genomes_trimmed.fasta")
skbio.io.registry.write(remove_unfit_contigs(fasta, min_length, min_coverage), 'fasta', trimmed_fasta)
os.chdir(out_dir)
sys.stdout.write(run_nucmer(trimmed_fasta, 50))
parse_delta_file("pairwise.delta", "pairwise_filtered.delta")
return 0
if __name__ == "__main__":
main()