annotate srst2.py @ 0:6f870ed59b6e draft

Uploaded
author nml
date Mon, 06 Feb 2017 12:31:04 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1 #!/usr/bin/env python
6f870ed59b6e Uploaded
nml
parents:
diff changeset
2
6f870ed59b6e Uploaded
nml
parents:
diff changeset
3 # SRST2 - Short Read Sequence Typer (v2)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
4 # Python Version 2.7.5
6f870ed59b6e Uploaded
nml
parents:
diff changeset
5 #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
6 # Authors - Michael Inouye (minouye@unimelb.edu.au), Harriet Dashnow (h.dashnow@gmail.com),
6f870ed59b6e Uploaded
nml
parents:
diff changeset
7 # Kathryn Holt (kholt@unimelb.edu.au), Bernie Pope (bjpope@unimelb.edu.au)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
8 #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
9 # see LICENSE.txt for the license
6f870ed59b6e Uploaded
nml
parents:
diff changeset
10 #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
11 # Dependencies:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
12 # bowtie2 http://bowtie-bio.sourceforge.net/bowtie2/index.shtml version 2.1.0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
13 # SAMtools http://samtools.sourceforge.net Version: 0.1.18 (Version: 0.1.19 DOES NOT WORK - loss of edge coverage)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
14 # SciPy http://www.scipy.org/install.html
6f870ed59b6e Uploaded
nml
parents:
diff changeset
15 #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
16 # Git repository: https://github.com/katholt/srst2/
6f870ed59b6e Uploaded
nml
parents:
diff changeset
17 # README: https://github.com/katholt/srst2/blob/master/README.md
6f870ed59b6e Uploaded
nml
parents:
diff changeset
18 # Questions or feature requests: https://github.com/katholt/srst2/issues
6f870ed59b6e Uploaded
nml
parents:
diff changeset
19 # Manuscript: http://biorxiv.org/content/early/2014/06/26/006627
6f870ed59b6e Uploaded
nml
parents:
diff changeset
20
6f870ed59b6e Uploaded
nml
parents:
diff changeset
21
6f870ed59b6e Uploaded
nml
parents:
diff changeset
22 from argparse import (ArgumentParser, FileType)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
23 import logging
6f870ed59b6e Uploaded
nml
parents:
diff changeset
24 from subprocess import call, check_output, CalledProcessError, STDOUT
6f870ed59b6e Uploaded
nml
parents:
diff changeset
25 import os, sys, re, collections, operator
6f870ed59b6e Uploaded
nml
parents:
diff changeset
26 from scipy.stats import binom_test, linregress
6f870ed59b6e Uploaded
nml
parents:
diff changeset
27 from math import log
6f870ed59b6e Uploaded
nml
parents:
diff changeset
28 from itertools import groupby
6f870ed59b6e Uploaded
nml
parents:
diff changeset
29 from operator import itemgetter
6f870ed59b6e Uploaded
nml
parents:
diff changeset
30 from collections import OrderedDict
6f870ed59b6e Uploaded
nml
parents:
diff changeset
31 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
32 from version import srst2_version
6f870ed59b6e Uploaded
nml
parents:
diff changeset
33 except:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
34 srst2_version = "version unknown"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
35
6f870ed59b6e Uploaded
nml
parents:
diff changeset
36 edge_a = edge_z = 2
6f870ed59b6e Uploaded
nml
parents:
diff changeset
37
6f870ed59b6e Uploaded
nml
parents:
diff changeset
38
6f870ed59b6e Uploaded
nml
parents:
diff changeset
39 def parse_args():
6f870ed59b6e Uploaded
nml
parents:
diff changeset
40 "Parse the input arguments, use '-h' for help."
6f870ed59b6e Uploaded
nml
parents:
diff changeset
41
6f870ed59b6e Uploaded
nml
parents:
diff changeset
42 parser = ArgumentParser(description='SRST2 - Short Read Sequence Typer (v2)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
43
6f870ed59b6e Uploaded
nml
parents:
diff changeset
44 # version number of srst2, print and then exit
6f870ed59b6e Uploaded
nml
parents:
diff changeset
45 parser.add_argument('--version', action='version', version='%(prog)s ' + srst2_version)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
46
6f870ed59b6e Uploaded
nml
parents:
diff changeset
47 # Read inputs
6f870ed59b6e Uploaded
nml
parents:
diff changeset
48 parser.add_argument(
6f870ed59b6e Uploaded
nml
parents:
diff changeset
49 '--input_se', nargs='+',type=str, required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
50 help='Single end read file(s) for analysing (may be gzipped)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
51 parser.add_argument(
6f870ed59b6e Uploaded
nml
parents:
diff changeset
52 '--input_pe', nargs='+', type=str, required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
53 help='Paired end read files for analysing (may be gzipped)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
54 parser.add_argument(
6f870ed59b6e Uploaded
nml
parents:
diff changeset
55 '--forward', type=str, required=False, default="_1",
6f870ed59b6e Uploaded
nml
parents:
diff changeset
56 help='Designator for forward reads (only used if NOT in MiSeq format sample_S1_L001_R1_001.fastq.gz; otherwise default is _1, i.e. expect forward reads as sample_1.fastq.gz)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
57 parser.add_argument(
6f870ed59b6e Uploaded
nml
parents:
diff changeset
58 '--reverse', type=str, required=False, default="_2",
6f870ed59b6e Uploaded
nml
parents:
diff changeset
59 help='Designator for reverse reads (only used if NOT in MiSeq format sample_S1_L001_R2_001.fastq.gz; otherwise default is _2, i.e. expect forward reads as sample_2.fastq.gz')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
60 parser.add_argument('--read_type', type=str, choices=['q', 'qseq', 'f'], default='q',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
61 help='Read file type (for bowtie2; default is q=fastq; other options: qseq=solexa, f=fasta).')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
62
6f870ed59b6e Uploaded
nml
parents:
diff changeset
63 # MLST parameters
6f870ed59b6e Uploaded
nml
parents:
diff changeset
64 parser.add_argument('--mlst_db', type=str, required=False, nargs=1, help='Fasta file of MLST alleles (optional)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
65 parser.add_argument('--mlst_delimiter', type=str, required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
66 help='Character(s) separating gene name from allele number in MLST database (default "-", as in arcc-1)', default="-")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
67 parser.add_argument('--mlst_definitions', type=str, required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
68 help='ST definitions for MLST scheme (required if mlst_db supplied and you want to calculate STs)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
69 parser.add_argument('--mlst_max_mismatch', type=str, required=False, default = "10",
6f870ed59b6e Uploaded
nml
parents:
diff changeset
70 help='Maximum number of mismatches per read for MLST allele calling (default 10)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
71
6f870ed59b6e Uploaded
nml
parents:
diff changeset
72 # Gene database parameters
6f870ed59b6e Uploaded
nml
parents:
diff changeset
73 parser.add_argument('--gene_db', type=str, required=False, nargs='+', help='Fasta file/s for gene databases (optional)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
74 parser.add_argument('--no_gene_details', action="store_false", required=False, help='Switch OFF verbose reporting of gene typing')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
75 parser.add_argument('--gene_max_mismatch', type=str, required=False, default = "10",
6f870ed59b6e Uploaded
nml
parents:
diff changeset
76 help='Maximum number of mismatches per read for gene detection and allele calling (default 10)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
77
6f870ed59b6e Uploaded
nml
parents:
diff changeset
78 # Cutoffs for scoring/heuristics
6f870ed59b6e Uploaded
nml
parents:
diff changeset
79 parser.add_argument('--min_coverage', type=float, required=False, help='Minimum %%coverage cutoff for gene reporting (default 90)',default=90)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
80 parser.add_argument('--max_divergence', type=float, required=False, help='Maximum %%divergence cutoff for gene reporting (default 10)',default=10)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
81 parser.add_argument('--min_depth', type=float, required=False, help='Minimum mean depth to flag as dubious allele call (default 5)',default=5)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
82 parser.add_argument('--min_edge_depth', type=float, required=False, help='Minimum edge depth to flag as dubious allele call (default 2)',default=2)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
83 parser.add_argument('--prob_err', type=float, default=0.01, help='Probability of sequencing error (default 0.01)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
84
6f870ed59b6e Uploaded
nml
parents:
diff changeset
85 # Mapping parameters for bowtie2
6f870ed59b6e Uploaded
nml
parents:
diff changeset
86 parser.add_argument('--stop_after', type=str, required=False, help='Stop mapping after this number of reads have been mapped (otherwise map all)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
87 parser.add_argument('--other', type=str, help='Other arguments to pass to bowtie2.', required=False)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
88
6f870ed59b6e Uploaded
nml
parents:
diff changeset
89 # Samtools parameters
6f870ed59b6e Uploaded
nml
parents:
diff changeset
90 parser.add_argument('--mapq', type=int, default=1, help='Samtools -q parameter (default 1)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
91 parser.add_argument('--baseq', type=int, default=20, help='Samtools -Q parameter (default 20)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
92
6f870ed59b6e Uploaded
nml
parents:
diff changeset
93 # Reporting options
6f870ed59b6e Uploaded
nml
parents:
diff changeset
94 parser.add_argument('--output', type=str, required=True, help='Prefix for srst2 output files')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
95 parser.add_argument('--log', action="store_true", required=False, help='Switch ON logging to file (otherwise log to stdout)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
96 parser.add_argument('--save_scores', action="store_true", required=False, help='Switch ON verbose reporting of all scores')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
97 parser.add_argument('--report_new_consensus', action="store_true", required=False, help='If a matching alleles is not found, report the consensus allele. Note, only SNP differences are considered, not indels.')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
98 parser.add_argument('--report_all_consensus', action="store_true", required=False, help='Report the consensus allele for the most likely allele. Note, only SNP differences are considered, not indels.')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
99
6f870ed59b6e Uploaded
nml
parents:
diff changeset
100 # Run options
6f870ed59b6e Uploaded
nml
parents:
diff changeset
101 parser.add_argument('--use_existing_pileup', action="store_true", required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
102 help='Use existing pileups if available, otherwise they will be generated') # to facilitate testing of rescoring from pileups
6f870ed59b6e Uploaded
nml
parents:
diff changeset
103 parser.add_argument('--use_existing_scores', action="store_true", required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
104 help='Use existing scores files if available, otherwise they will be generated') # to facilitate testing of reporting from scores
6f870ed59b6e Uploaded
nml
parents:
diff changeset
105 parser.add_argument('--keep_interim_alignment', action="store_true", required=False, default=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
106 help='Keep interim files (sam & unsorted bam), otherwise they will be deleted after sorted bam is created') # to facilitate testing of sam processing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
107 # parser.add_argument('--keep_final_alignment', action="store_true", required=False, default=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
108 # help='Keep interim files (sam & unsorted bam), otherwise they will be deleted after sorted bam is created') # to facilitate testing of sam processing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
109
6f870ed59b6e Uploaded
nml
parents:
diff changeset
110 # Compile previous output files
6f870ed59b6e Uploaded
nml
parents:
diff changeset
111 parser.add_argument('--prev_output', nargs='+', type=str, required=False,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
112 help='SRST2 results files to compile (any new results from this run will also be incorporated)')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
113
6f870ed59b6e Uploaded
nml
parents:
diff changeset
114 return parser.parse_args()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
115
6f870ed59b6e Uploaded
nml
parents:
diff changeset
116
6f870ed59b6e Uploaded
nml
parents:
diff changeset
117 # Exception to raise if the command we try to run fails for some reason
6f870ed59b6e Uploaded
nml
parents:
diff changeset
118 class CommandError(Exception):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
119 pass
6f870ed59b6e Uploaded
nml
parents:
diff changeset
120
6f870ed59b6e Uploaded
nml
parents:
diff changeset
121 def run_command(command, **kwargs):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
122 'Execute a shell command and check the exit status and any O/S exceptions'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
123 command_str = ' '.join(command)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
124 logging.info('Running: {}'.format(command_str))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
125 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
126 exit_status = call(command, **kwargs)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
127 except OSError as e:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
128 message = "Command '{}' failed due to O/S error: {}".format(command_str, str(e))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
129 raise CommandError({"message": message})
6f870ed59b6e Uploaded
nml
parents:
diff changeset
130 if exit_status != 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
131 message = "Command '{}' failed with non-zero exit status: {}".format(command_str, exit_status)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
132 raise CommandError({"message": message})
6f870ed59b6e Uploaded
nml
parents:
diff changeset
133
6f870ed59b6e Uploaded
nml
parents:
diff changeset
134
6f870ed59b6e Uploaded
nml
parents:
diff changeset
135 def bowtie_index(fasta_files):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
136 'Build a bowtie2 index from the given input fasta(s)'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
137
6f870ed59b6e Uploaded
nml
parents:
diff changeset
138 # check that both bowtie and samtools have the right versions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
139 check_command_version(['bowtie2', '--version'],
6f870ed59b6e Uploaded
nml
parents:
diff changeset
140 'bowtie2-align version 2.1.0',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
141 'bowtie',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
142 '2.1.0')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
143
6f870ed59b6e Uploaded
nml
parents:
diff changeset
144 for fasta in fasta_files:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
145 built_index = fasta + '.1.bt2'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
146 if os.path.exists(built_index):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
147 logging.info('Index for {} is already built...'.format(fasta))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
148 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
149 logging.info('Building bowtie2 index for {}...'.format(fasta))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
150 run_command(['bowtie2-build', fasta, fasta])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
151
6f870ed59b6e Uploaded
nml
parents:
diff changeset
152
6f870ed59b6e Uploaded
nml
parents:
diff changeset
153 def modify_bowtie_sam(raw_bowtie_sam,max_mismatch):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
154 # fix sam flags for comprehensive pileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
155 with open(raw_bowtie_sam) as sam, open(raw_bowtie_sam + '.mod', 'w') as sam_mod:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
156 for line in sam:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
157 if not line.startswith('@'):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
158 fields = line.split('\t')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
159 flag = int(fields[1])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
160 flag = (flag - 256) if (flag & 256) else flag
6f870ed59b6e Uploaded
nml
parents:
diff changeset
161 m = re.search("NM:i:(\d+)\s",line)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
162 if m != None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
163 num_mismatch = m.group(1)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
164 if int(num_mismatch) <= int(max_mismatch):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
165 sam_mod.write('\t'.join([fields[0], str(flag)] + fields[2:]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
166 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
167 logging.info('Excluding read from SAM file due to missing NM (num mismatches) field: ' + fields[0])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
168 num_mismatch = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
169 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
170 sam_mod.write(line)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
171 return(raw_bowtie_sam,raw_bowtie_sam + '.mod')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
172
6f870ed59b6e Uploaded
nml
parents:
diff changeset
173
6f870ed59b6e Uploaded
nml
parents:
diff changeset
174 def parse_fai(fai_file,db_type,delimiter):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
175 'Get sequence lengths for reference alleles - important for scoring'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
176 'Get gene names also, required if no MLST definitions provided'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
177 size = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
178 gene_clusters = [] # for gene DBs, this is cluster ID
6f870ed59b6e Uploaded
nml
parents:
diff changeset
179 allele_symbols = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
180 gene_cluster_symbols = {} # key = cluster ID, value = gene symbol (for gene DBs)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
181 unique_allele_symbols = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
182 unique_gene_symbols = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
183 delimiter_check = [] # list of names that may violate the MLST delimiter supplied
6f870ed59b6e Uploaded
nml
parents:
diff changeset
184 with open(fai_file) as fai:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
185 for line in fai:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
186 fields = line.split('\t')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
187 name = fields[0] # full allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
188 size[name] = int(fields[1]) # store length
6f870ed59b6e Uploaded
nml
parents:
diff changeset
189 if db_type!="mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
190 allele_info = name.split()[0].split("__")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
191 if len(allele_info) > 2:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
192 gene_cluster = allele_info[0] # ID number for the cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
193 cluster_symbol = allele_info[1] # gene name for the cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
194 name = allele_info[2] # specific allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
195 if gene_cluster in gene_cluster_symbols:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
196 if gene_cluster_symbols[gene_cluster] != cluster_symbol:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
197 unique_gene_symbols = False # already seen this cluster symbol
6f870ed59b6e Uploaded
nml
parents:
diff changeset
198 logging.info( "Non-unique:" + gene_cluster + ", " + cluster_symbol)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
199 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
200 gene_cluster_symbols[gene_cluster] = cluster_symbol
6f870ed59b6e Uploaded
nml
parents:
diff changeset
201 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
202 # treat as unclustered database, use whole header
6f870ed59b6e Uploaded
nml
parents:
diff changeset
203 gene_cluster = cluster_symbol = name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
204 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
205 gene_cluster = name.split(delimiter)[0] # accept gene clusters raw for mlst
6f870ed59b6e Uploaded
nml
parents:
diff changeset
206 # check if the delimiter makes sense
6f870ed59b6e Uploaded
nml
parents:
diff changeset
207 parts = name.split(delimiter)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
208 if len(parts) != 2:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
209 delimiter_check.append(name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
210 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
211 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
212 x = int(parts[1])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
213 except:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
214 delimiter_check.append(name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
215
6f870ed59b6e Uploaded
nml
parents:
diff changeset
216 # check if we have seen this allele name before
6f870ed59b6e Uploaded
nml
parents:
diff changeset
217 if name in allele_symbols:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
218 unique_allele_symbols = False # already seen this allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
219 allele_symbols.append(name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
220
6f870ed59b6e Uploaded
nml
parents:
diff changeset
221 # record gene (cluster):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
222 if gene_cluster not in gene_clusters:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
223 gene_clusters.append(gene_cluster)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
224
6f870ed59b6e Uploaded
nml
parents:
diff changeset
225 if len(delimiter_check) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
226 print "Warning! MLST delimiter is " + delimiter + " but these genes may violate the pattern and cause problems:"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
227 print ",".join(delimiter_check)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
228
6f870ed59b6e Uploaded
nml
parents:
diff changeset
229 return size, gene_clusters, unique_gene_symbols, unique_allele_symbols, gene_cluster_symbols
6f870ed59b6e Uploaded
nml
parents:
diff changeset
230
6f870ed59b6e Uploaded
nml
parents:
diff changeset
231
6f870ed59b6e Uploaded
nml
parents:
diff changeset
232 def read_pileup_data(pileup_file, size, prob_err, consensus_file = ""):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
233 with open(pileup_file) as pileup:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
234 prob_success = 1 - prob_err # Set by user, default is prob_err = 0.01
6f870ed59b6e Uploaded
nml
parents:
diff changeset
235 hash_alignment = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
236 hash_max_depth = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
237 hash_edge_depth = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
238 max_depth = 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
239 avg_depth_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
240 next_to_del_depth_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
241 coverage_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
242 mismatch_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
243 indel_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
244 missing_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
245 size_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
246
6f870ed59b6e Uploaded
nml
parents:
diff changeset
247 # Split all lines in the pileup by whitespace
6f870ed59b6e Uploaded
nml
parents:
diff changeset
248 pileup_split = ( x.split() for x in pileup )
6f870ed59b6e Uploaded
nml
parents:
diff changeset
249 # Group the split lines based on the first field (allele)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
250 for allele, lines in groupby(pileup_split, itemgetter(0)):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
251
6f870ed59b6e Uploaded
nml
parents:
diff changeset
252 # Reset variables for new allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
253 allele_line = 1 # Keep track of line for this allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
254 exp_nuc_num = 0 # Expected position in ref allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
255 allele_size = size[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
256 total_depth = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
257 depth_a = depth_z = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
258 position_depths = [0] * allele_size # store depths in case required for penalties; then we don't need to track total_missing_bases
6f870ed59b6e Uploaded
nml
parents:
diff changeset
259 hash_alignment[allele] = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
260 total_missing_bases = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
261 total_mismatch = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
262 ins_poscount = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
263 del_poscount = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
264 next_to_del_depth = 99999
6f870ed59b6e Uploaded
nml
parents:
diff changeset
265 consensus_seq = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
266
6f870ed59b6e Uploaded
nml
parents:
diff changeset
267 for fields in lines:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
268 # Parse this line and store details required for scoring
6f870ed59b6e Uploaded
nml
parents:
diff changeset
269 nuc_num = int(fields[1]) # Actual position in ref allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
270 exp_nuc_num += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
271 allele_line += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
272 nuc = fields[2]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
273 nuc_depth = int(fields[3])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
274 position_depths[nuc_num-1] = nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
275 if len(fields) <= 5:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
276 aligned_bases = ''
6f870ed59b6e Uploaded
nml
parents:
diff changeset
277 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
278 aligned_bases = fields[4]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
279
6f870ed59b6e Uploaded
nml
parents:
diff changeset
280 # Missing bases (pileup skips basepairs)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
281 if nuc_num > exp_nuc_num:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
282 total_missing_bases += abs(exp_nuc_num - nuc_num)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
283 exp_nuc_num = nuc_num
6f870ed59b6e Uploaded
nml
parents:
diff changeset
284 if nuc_depth == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
285 total_missing_bases += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
286
6f870ed59b6e Uploaded
nml
parents:
diff changeset
287 # Calculate depths for this position
6f870ed59b6e Uploaded
nml
parents:
diff changeset
288 if nuc_num <= edge_a:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
289 depth_a += nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
290 if abs(nuc_num - allele_size) < edge_z:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
291 depth_z += nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
292 if nuc_depth > max_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
293 hash_max_depth[allele] = nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
294 max_depth = nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
295
6f870ed59b6e Uploaded
nml
parents:
diff changeset
296 total_depth = total_depth + nuc_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
297
6f870ed59b6e Uploaded
nml
parents:
diff changeset
298 # Parse aligned bases list for this position in the pileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
299 num_match = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
300 ins_readcount = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
301 del_readcount = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
302 nuc_counts = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
303
6f870ed59b6e Uploaded
nml
parents:
diff changeset
304 i = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
305 while i < len(aligned_bases):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
306
6f870ed59b6e Uploaded
nml
parents:
diff changeset
307 if aligned_bases[i] == "^":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
308 # Signifies start of a read, next char is mapping quality (skip it)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
309 i += 2
6f870ed59b6e Uploaded
nml
parents:
diff changeset
310 continue
6f870ed59b6e Uploaded
nml
parents:
diff changeset
311
6f870ed59b6e Uploaded
nml
parents:
diff changeset
312 if aligned_bases[i] == "+":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
313 i += int(aligned_bases[i+1]) + 2 # skip to next read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
314 ins_readcount += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
315 continue
6f870ed59b6e Uploaded
nml
parents:
diff changeset
316
6f870ed59b6e Uploaded
nml
parents:
diff changeset
317 if aligned_bases[i] == "-":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
318 i += int(aligned_bases[i+1]) + 2 # skip to next read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
319 continue
6f870ed59b6e Uploaded
nml
parents:
diff changeset
320
6f870ed59b6e Uploaded
nml
parents:
diff changeset
321 if aligned_bases[i] == "*":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
322 i += 1 # skip to next read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
323 del_readcount += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
324 continue
6f870ed59b6e Uploaded
nml
parents:
diff changeset
325
6f870ed59b6e Uploaded
nml
parents:
diff changeset
326 if aligned_bases[i] == "." or aligned_bases[i] == ",":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
327 num_match += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
328 i += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
329 continue
6f870ed59b6e Uploaded
nml
parents:
diff changeset
330
6f870ed59b6e Uploaded
nml
parents:
diff changeset
331 elif aligned_bases[i].upper() in "ATCG":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
332 this_nuc = aligned_bases[i].upper()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
333 if this_nuc not in nuc_counts:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
334 nuc_counts[this_nuc] = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
335 nuc_counts[this_nuc] += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
336
6f870ed59b6e Uploaded
nml
parents:
diff changeset
337 i += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
338
6f870ed59b6e Uploaded
nml
parents:
diff changeset
339 # Save the most common nucleotide at this position
6f870ed59b6e Uploaded
nml
parents:
diff changeset
340 consensus_nuc = nuc # by default use reference nucleotide
6f870ed59b6e Uploaded
nml
parents:
diff changeset
341 max_freq = num_match # Number of bases matching the reference
6f870ed59b6e Uploaded
nml
parents:
diff changeset
342 for nucleotide in nuc_counts:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
343 if nuc_counts[nucleotide] > max_freq:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
344 consensus_nuc = nucleotide
6f870ed59b6e Uploaded
nml
parents:
diff changeset
345 max_freq = nuc_counts[nucleotide]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
346 consensus_seq += (consensus_nuc)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
347
6f870ed59b6e Uploaded
nml
parents:
diff changeset
348 # Calculate details of this position for scoring and reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
349
6f870ed59b6e Uploaded
nml
parents:
diff changeset
350 # mismatches and indels
6f870ed59b6e Uploaded
nml
parents:
diff changeset
351 num_mismatch = nuc_depth - num_match
6f870ed59b6e Uploaded
nml
parents:
diff changeset
352 if num_mismatch > num_match:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
353 total_mismatch += 1 # record as mismatch (could be a snp or deletion)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
354 if del_readcount > num_match:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
355 del_poscount += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
356 if ins_readcount > nuc_depth / 2:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
357 ins_poscount += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
358
6f870ed59b6e Uploaded
nml
parents:
diff changeset
359 # Hash for later processing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
360 hash_alignment[allele].append((num_match, num_mismatch, prob_success)) # snp or deletion
6f870ed59b6e Uploaded
nml
parents:
diff changeset
361 if ins_readcount > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
362 hash_alignment[allele].append((nuc_depth - ins_readcount, ins_readcount, prob_success)) # penalize for any insertion calls at this position
6f870ed59b6e Uploaded
nml
parents:
diff changeset
363
6f870ed59b6e Uploaded
nml
parents:
diff changeset
364 # Determine the consensus sequence if required
6f870ed59b6e Uploaded
nml
parents:
diff changeset
365 if consensus_file != "":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
366 if consensus_file.split(".")[-2] == "new_consensus_alleles":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
367 consensus_type = "variant"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
368 elif consensus_file.split(".")[-2] == "all_consensus_alleles":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
369 consensus_type = "consensus"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
370 with open(consensus_file, "a") as consensus_outfile:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
371 consensus_outfile.write(">{0}.{1} {2}\n".format(allele, consensus_type, pileup_file.split(".")[1].split("__")[1]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
372 outstring = consensus_seq + "\n"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
373 consensus_outfile.write(outstring)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
374
6f870ed59b6e Uploaded
nml
parents:
diff changeset
375 # Finished reading pileup for this allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
376
6f870ed59b6e Uploaded
nml
parents:
diff changeset
377 # Check for missing bases at the end of the allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
378 if nuc_num < allele_size:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
379 total_missing_bases += abs(allele_size - nuc_num)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
380 # determine penalty based on coverage of last 2 bases
6f870ed59b6e Uploaded
nml
parents:
diff changeset
381 penalty = float(position_depths[nuc_num-1] + position_depths[nuc_num-2])/2
6f870ed59b6e Uploaded
nml
parents:
diff changeset
382 m = min(position_depths[nuc_num-1],position_depths[nuc_num-2])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
383 hash_alignment[allele].append((0, penalty, prob_success))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
384 if next_to_del_depth > m:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
385 next_to_del_depth = m # keep track of lowest near-del depth for reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
386
6f870ed59b6e Uploaded
nml
parents:
diff changeset
387 # Calculate allele summary stats and save
6f870ed59b6e Uploaded
nml
parents:
diff changeset
388 avg_depth = round(total_depth / float(allele_line),3)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
389 avg_a = depth_a / float(edge_a) # Avg depth at 5' end, num basepairs determined by edge_a
6f870ed59b6e Uploaded
nml
parents:
diff changeset
390 avg_z = depth_z / float(edge_z) # 3'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
391 hash_max_depth[allele] = max_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
392 hash_edge_depth[allele] = (avg_a, avg_z)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
393 min_penalty = max(5, int(avg_depth))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
394 coverage_allele[allele] = 100*(allele_size - total_missing_bases - del_poscount)/float(allele_size) # includes in-read deletions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
395 mismatch_allele[allele] = total_mismatch - del_poscount # snps only
6f870ed59b6e Uploaded
nml
parents:
diff changeset
396 indel_allele[allele] = del_poscount + ins_poscount # insertions or deletions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
397 missing_allele[allele] = total_missing_bases # truncated bases
6f870ed59b6e Uploaded
nml
parents:
diff changeset
398 size_allele[allele] = allele_size
6f870ed59b6e Uploaded
nml
parents:
diff changeset
399
6f870ed59b6e Uploaded
nml
parents:
diff changeset
400 # Penalize truncations or large deletions (i.e. positions not covered in pileup)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
401 j = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
402 while j < (len(position_depths)-2):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
403 # note end-of-seq truncations are dealt with above)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
404 if position_depths[j]==0 and position_depths[j+1]!=0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
405 penalty = float(position_depths[j+1]+position_depths[j+2])/2 # mean of next 2 bases
6f870ed59b6e Uploaded
nml
parents:
diff changeset
406 hash_alignment[allele].append((0, penalty, prob_success))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
407 m = min(position_depths[nuc_num-1],position_depths[nuc_num-2])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
408 if next_to_del_depth > m:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
409 next_to_del_depth = m # keep track of lowest near-del depth for reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
410 j += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
411
6f870ed59b6e Uploaded
nml
parents:
diff changeset
412 # Store depth info for reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
413 avg_depth_allele[allele] = avg_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
414 if next_to_del_depth == 99999:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
415 next_to_del_depth = "NA"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
416 next_to_del_depth_allele[allele] = next_to_del_depth
6f870ed59b6e Uploaded
nml
parents:
diff changeset
417
6f870ed59b6e Uploaded
nml
parents:
diff changeset
418 return hash_alignment, hash_max_depth, hash_edge_depth, avg_depth_allele, coverage_allele, mismatch_allele, indel_allele, missing_allele, size_allele, next_to_del_depth_allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
419
6f870ed59b6e Uploaded
nml
parents:
diff changeset
420
6f870ed59b6e Uploaded
nml
parents:
diff changeset
421 def score_alleles(args, mapping_files_pre, hash_alignment, hash_max_depth, hash_edge_depth,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
422 avg_depth_allele, coverage_allele, mismatch_allele, indel_allele, missing_allele,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
423 size_allele, next_to_del_depth_allele, run_type):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
424
6f870ed59b6e Uploaded
nml
parents:
diff changeset
425 if args.save_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
426 scores_output = file(mapping_files_pre + '.scores', 'w')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
427 scores_output.write("Allele\tScore\tAvg_depth\tEdge1_depth\tEdge2_depth\tPercent_coverage\tSize\tMismatches\tIndels\tTruncated_bases\tDepthNeighbouringTruncation\tmaxMAF\tLeastConfident_Rate\tLeastConfident_Mismatches\tLeastConfident_Depth\tLeastConfident_Pvalue\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
428
6f870ed59b6e Uploaded
nml
parents:
diff changeset
429 scores = {} # key = allele, value = score
6f870ed59b6e Uploaded
nml
parents:
diff changeset
430 mix_rates = {} # key = allele, value = highest minor allele frequency, 0 -> 0.5
6f870ed59b6e Uploaded
nml
parents:
diff changeset
431
6f870ed59b6e Uploaded
nml
parents:
diff changeset
432 for allele in hash_alignment:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
433 if (run_type == "mlst") or (coverage_allele[allele] > args.min_coverage):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
434 pvals = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
435 min_pval = 1.0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
436 min_pval_data = (999,999) # (mismatch, depth) for position with lowest p-value
6f870ed59b6e Uploaded
nml
parents:
diff changeset
437 mix_rate = 0 # highest minor allele frequency 0 -> 0.5
6f870ed59b6e Uploaded
nml
parents:
diff changeset
438 for nuc_info in hash_alignment[allele]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
439 if nuc_info is not None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
440 match, mismatch, prob_success = nuc_info
6f870ed59b6e Uploaded
nml
parents:
diff changeset
441 if match > 0 or mismatch > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
442 if mismatch == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
443 p_value = 1.0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
444 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
445 p_value = binom_test([match, mismatch], None, prob_success)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
446 # Weight pvalue by (depth/max_depth)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
447 max_depth = hash_max_depth[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
448 weight = (match + mismatch) / float(max_depth)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
449 p_value *= weight
6f870ed59b6e Uploaded
nml
parents:
diff changeset
450 if p_value < min_pval:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
451 min_pval = p_value
6f870ed59b6e Uploaded
nml
parents:
diff changeset
452 min_pval_data = (mismatch,match + mismatch)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
453 if p_value > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
454 p_value = -log(p_value, 10)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
455 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
456 p_value = 1000
6f870ed59b6e Uploaded
nml
parents:
diff changeset
457 pvals.append(p_value)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
458 mismatch_prop = float(match)/float(match+mismatch)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
459 if min(mismatch_prop, 1-mismatch_prop) > mix_rate:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
460 mix_rate = min(mismatch_prop, 1-mismatch_prop)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
461 # Fit linear model to observed Pval distribution vs expected Pval distribution (QQ plot)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
462 pvals.sort(reverse=True)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
463 len_obs_pvals = len(pvals)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
464 exp_pvals = range(1, len_obs_pvals + 1)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
465 exp_pvals2 = [-log(float(ep) / (len_obs_pvals + 1), 10) for ep in exp_pvals]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
466
6f870ed59b6e Uploaded
nml
parents:
diff changeset
467 # Slope is score
6f870ed59b6e Uploaded
nml
parents:
diff changeset
468 slope, _intercept, _r_value, _p_value, _std_err = linregress(exp_pvals2, pvals)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
469
6f870ed59b6e Uploaded
nml
parents:
diff changeset
470 # Store all scores for later processing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
471 scores[allele] = slope
6f870ed59b6e Uploaded
nml
parents:
diff changeset
472 mix_rates[allele] = mix_rate
6f870ed59b6e Uploaded
nml
parents:
diff changeset
473
6f870ed59b6e Uploaded
nml
parents:
diff changeset
474 # print scores for each allele, if requested
6f870ed59b6e Uploaded
nml
parents:
diff changeset
475 if args.save_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
476 if allele in hash_edge_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
477 start_depth, end_depth = hash_edge_depth[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
478 edge_depth_str = str(start_depth) + '\t' + str(end_depth)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
479 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
480 edge_depth_str = "NA\tNA"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
481 this_depth = avg_depth_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
482 this_coverage = coverage_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
483 this_mismatch = mismatch_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
484 this_indel = indel_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
485 this_missing = missing_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
486 this_size = size_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
487 this_next_to_del_depth = next_to_del_depth_allele.get(allele, "NA")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
488 scores_output.write('\t'.join([allele, str(slope), str(this_depth), edge_depth_str,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
489 str(this_coverage), str(this_size), str(this_mismatch), str(this_indel), str(this_missing), str(this_next_to_del_depth), str(mix_rate), str(float(min_pval_data[0])/min_pval_data[1]),str(min_pval_data[0]),str(min_pval_data[1]),str(min_pval)]) + '\n')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
490
6f870ed59b6e Uploaded
nml
parents:
diff changeset
491 if args.save_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
492 scores_output.close()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
493
6f870ed59b6e Uploaded
nml
parents:
diff changeset
494 return(scores,mix_rates)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
495
6f870ed59b6e Uploaded
nml
parents:
diff changeset
496 # Check that an acceptable version of a command is installed
6f870ed59b6e Uploaded
nml
parents:
diff changeset
497 # Exits the program if it can't be found.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
498 # - command_list is the command to run to determine the version.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
499 # - version_identifier is the unique string we look for in the stdout of the program.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
500 # - command_name is the name of the command to show in error messages.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
501 # - required_version is the version number to show in error messages.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
502 def check_command_version(command_list, version_identifier, command_name, required_version):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
503 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
504 command_stdout = check_output(command_list, stderr=STDOUT)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
505 except OSError as e:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
506 logging.error("Failed command: {}".format(' '.join(command_list)))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
507 logging.error(str(e))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
508 logging.error("Could not determine the version of {}.".format(command_name))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
509 logging.error("Do you have {} installed in your PATH?".format(command_name))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
510 exit(-1)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
511 except CalledProcessError as e:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
512 # some programs such as samtools return a non-zero exit status
6f870ed59b6e Uploaded
nml
parents:
diff changeset
513 # when you ask for the version (sigh). We ignore it here.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
514 command_stdout = e.output
6f870ed59b6e Uploaded
nml
parents:
diff changeset
515
6f870ed59b6e Uploaded
nml
parents:
diff changeset
516 if version_identifier not in command_stdout:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
517 logging.error("Incorrect version of {} installed.".format(command_name))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
518 logging.error("{} version {} is required by SRST2.".format(command_name, required_version))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
519 exit(-1)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
520
6f870ed59b6e Uploaded
nml
parents:
diff changeset
521
6f870ed59b6e Uploaded
nml
parents:
diff changeset
522 def run_bowtie(mapping_files_pre,sample_name,fastqs,args,db_name,db_full_path):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
523
6f870ed59b6e Uploaded
nml
parents:
diff changeset
524 print "Starting mapping with bowtie2"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
525
6f870ed59b6e Uploaded
nml
parents:
diff changeset
526 # check that both bowtie and samtools have the right versions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
527 check_command_version(['bowtie2', '--version'],
6f870ed59b6e Uploaded
nml
parents:
diff changeset
528 'bowtie2-align version 2.1.0',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
529 'bowtie',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
530 '2.1.0')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
531
6f870ed59b6e Uploaded
nml
parents:
diff changeset
532 check_command_version(['samtools'],
6f870ed59b6e Uploaded
nml
parents:
diff changeset
533 'Version: 0.1.18',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
534 'samtools',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
535 '0.1.18')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
536
6f870ed59b6e Uploaded
nml
parents:
diff changeset
537 command = ['bowtie2']
6f870ed59b6e Uploaded
nml
parents:
diff changeset
538
6f870ed59b6e Uploaded
nml
parents:
diff changeset
539 if len(fastqs)==1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
540 # single end
6f870ed59b6e Uploaded
nml
parents:
diff changeset
541 command += ['-U', fastqs[0]]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
542 elif len(fastqs)==2:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
543 # paired end
6f870ed59b6e Uploaded
nml
parents:
diff changeset
544 command += ['-1', fastqs[0], '-2', fastqs[1]]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
545
6f870ed59b6e Uploaded
nml
parents:
diff changeset
546 sam = mapping_files_pre + ".sam"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
547 logging.info('Output prefix set to: ' + mapping_files_pre)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
548
6f870ed59b6e Uploaded
nml
parents:
diff changeset
549 command += ['-S', sam,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
550 '-' + args.read_type, # add a dash to the front of the option
6f870ed59b6e Uploaded
nml
parents:
diff changeset
551 '--very-sensitive-local',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
552 '--no-unal',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
553 '-a', # Search for and report all alignments
6f870ed59b6e Uploaded
nml
parents:
diff changeset
554 '-x', db_full_path # The index to be aligned to
6f870ed59b6e Uploaded
nml
parents:
diff changeset
555 ]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
556
6f870ed59b6e Uploaded
nml
parents:
diff changeset
557 if args.stop_after:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
558 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
559 command += ['-u',str(int(args.stop_after))]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
560 except ValueError:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
561 print "WARNING. You asked to stop after mapping '" + args.stop_after + "' reads. I don't understand this, and will map all reads. Please speficy an integer with --stop_after or leave this as default to map 1 million reads."
6f870ed59b6e Uploaded
nml
parents:
diff changeset
562
6f870ed59b6e Uploaded
nml
parents:
diff changeset
563 if args.other:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
564 command += args.other.split()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
565
6f870ed59b6e Uploaded
nml
parents:
diff changeset
566 logging.info('Aligning reads to index {} using bowtie2...'.format(db_full_path))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
567
6f870ed59b6e Uploaded
nml
parents:
diff changeset
568 run_command(command)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
569
6f870ed59b6e Uploaded
nml
parents:
diff changeset
570 return(sam)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
571
6f870ed59b6e Uploaded
nml
parents:
diff changeset
572 def get_pileup(args,mapping_files_pre,raw_bowtie_sam,bowtie_sam_mod,fasta,pileup_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
573 # Analyse output with SAMtools
6f870ed59b6e Uploaded
nml
parents:
diff changeset
574 logging.info('Processing Bowtie2 output with SAMtools...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
575 logging.info('Generate and sort BAM file...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
576 out_file_bam = mapping_files_pre + ".unsorted.bam"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
577 run_command(['samtools', 'view', '-b', '-o', out_file_bam,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
578 '-q', str(args.mapq), '-S', bowtie_sam_mod])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
579 out_file_bam_sorted = mapping_files_pre + ".sorted"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
580 run_command(['samtools', 'sort', out_file_bam, out_file_bam_sorted])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
581
6f870ed59b6e Uploaded
nml
parents:
diff changeset
582 # Delete interim files (sam, modified sam, unsorted bam) unless otherwise specified.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
583 # Note users may also want to delete final sorted bam and pileup on completion to save space.
6f870ed59b6e Uploaded
nml
parents:
diff changeset
584 if not args.keep_interim_alignment:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
585 logging.info('Deleting sam and bam files that are not longer needed...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
586 del_filenames = [raw_bowtie_sam, bowtie_sam_mod, out_file_bam]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
587 for f in del_filenames:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
588 logging.info('Deleting ' + f)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
589 os.remove(f)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
590
6f870ed59b6e Uploaded
nml
parents:
diff changeset
591 logging.info('Generate pileup...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
592 with open(pileup_file, 'w') as sam_pileup:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
593 run_command(['samtools', 'mpileup', '-L', '1000', '-f', fasta,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
594 '-Q', str(args.baseq), '-q', str(args.mapq), out_file_bam_sorted + '.bam'],
6f870ed59b6e Uploaded
nml
parents:
diff changeset
595 stdout=sam_pileup)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
596
6f870ed59b6e Uploaded
nml
parents:
diff changeset
597 def calculate_ST(allele_scores, ST_db, gene_names, sample_name, mlst_delimiter, avg_depth_allele, mix_rates):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
598 allele_numbers = [] # clean allele calls for determing ST. order is taken from gene names, as in ST definitions file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
599 alleles_with_flags = [] # flagged alleles for printing (* if mismatches, ? if depth issues)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
600 mismatch_flags = [] # allele/diffs
6f870ed59b6e Uploaded
nml
parents:
diff changeset
601 uncertainty_flags = [] #allele/uncertainty
6f870ed59b6e Uploaded
nml
parents:
diff changeset
602 # st_flags = [] # (* if mismatches, ? if depth issues)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
603 depths = [] # depths for each typed locus
6f870ed59b6e Uploaded
nml
parents:
diff changeset
604 mafs = [] # minor allele freqencies for each typed locus
6f870ed59b6e Uploaded
nml
parents:
diff changeset
605
6f870ed59b6e Uploaded
nml
parents:
diff changeset
606 # get allele numbers & info
6f870ed59b6e Uploaded
nml
parents:
diff changeset
607 for gene in gene_names:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
608 if gene in allele_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
609 (allele,diffs,depth_problem,divergence) = allele_scores[gene]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
610 allele_number = allele.split(mlst_delimiter)[-1]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
611 depths.append(avg_depth_allele[allele])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
612 mix_rate = mix_rates[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
613 mafs.append(mix_rate)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
614 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
615 allele_number = "-"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
616 diffs = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
617 depth_problem = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
618 mix_rate = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
619 allele_numbers.append(allele_number)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
620
6f870ed59b6e Uploaded
nml
parents:
diff changeset
621 allele_with_flags = allele_number
6f870ed59b6e Uploaded
nml
parents:
diff changeset
622 if diffs != "":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
623 if diffs != "trun":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
624 allele_with_flags+="*" # trun indicates only that a truncated form had lower score, which isn't a mismatch
6f870ed59b6e Uploaded
nml
parents:
diff changeset
625 mismatch_flags.append(allele+"/"+diffs)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
626 if depth_problem != "":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
627 allele_with_flags+="?"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
628 uncertainty_flags.append(allele+"/"+depth_problem)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
629 alleles_with_flags.append(allele_with_flags)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
630
6f870ed59b6e Uploaded
nml
parents:
diff changeset
631 # calculate ST (no flags)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
632 if ST_db:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
633 allele_string = " ".join(allele_numbers) # for determining ST
6f870ed59b6e Uploaded
nml
parents:
diff changeset
634 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
635 clean_st = ST_db[allele_string]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
636 except KeyError:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
637 print "This combination of alleles was not found in the sequence type database:",
6f870ed59b6e Uploaded
nml
parents:
diff changeset
638 print sample_name,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
639 for gene in allele_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
640 (allele,diffs,depth_problems,divergence) = allele_scores[gene]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
641 print allele,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
642 print
6f870ed59b6e Uploaded
nml
parents:
diff changeset
643 clean_st = "NF"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
644 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
645 clean_st = "ND"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
646
6f870ed59b6e Uploaded
nml
parents:
diff changeset
647 # add flags for reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
648 st = clean_st
6f870ed59b6e Uploaded
nml
parents:
diff changeset
649 if len(mismatch_flags) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
650 if mismatch_flags!=["trun"]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
651 st += "*" # trun indicates only that a truncated form had lower score, which isn't a mismatch
6f870ed59b6e Uploaded
nml
parents:
diff changeset
652 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
653 mismatch_flags = ['0'] # record no mismatches
6f870ed59b6e Uploaded
nml
parents:
diff changeset
654 if len(uncertainty_flags) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
655 st += "?"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
656 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
657 uncertainty_flags = ['-']
6f870ed59b6e Uploaded
nml
parents:
diff changeset
658
6f870ed59b6e Uploaded
nml
parents:
diff changeset
659 # mean depth across loci
6f870ed59b6e Uploaded
nml
parents:
diff changeset
660 if len(depths) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
661 mean_depth = float(sum(depths))/len(depths)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
662 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
663 mean_depth = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
664
6f870ed59b6e Uploaded
nml
parents:
diff changeset
665 # maximum maf across locus
6f870ed59b6e Uploaded
nml
parents:
diff changeset
666 if len(mafs) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
667 max_maf = max(mafs)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
668 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
669 max_maf = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
670
6f870ed59b6e Uploaded
nml
parents:
diff changeset
671 return (st,clean_st,alleles_with_flags,mismatch_flags,uncertainty_flags,mean_depth,max_maf)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
672
6f870ed59b6e Uploaded
nml
parents:
diff changeset
673 def parse_ST_database(ST_filename,gene_names_from_fai):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
674 # Read ST definitions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
675 ST_db = {} # key = allele string, value = ST
6f870ed59b6e Uploaded
nml
parents:
diff changeset
676 gene_names = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
677 num_gene_cols_expected = len(gene_names_from_fai)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
678 print "Attempting to read " + str(num_gene_cols_expected) + " loci from ST database " + ST_filename
6f870ed59b6e Uploaded
nml
parents:
diff changeset
679 with open(ST_filename) as f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
680 count = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
681 for line in f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
682 count += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
683 line_split = line.rstrip().split("\t")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
684 if count == 1: # Header
6f870ed59b6e Uploaded
nml
parents:
diff changeset
685 gene_names = line_split[1:min(num_gene_cols_expected+1,len(line_split))]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
686 for g in gene_names_from_fai:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
687 if g not in gene_names:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
688 print "Warning: gene " + g + " in database file isn't among the columns in the ST definitions: " + ",".join(gene_names)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
689 print " Any sequences with this gene identifer from the database will not be included in typing."
6f870ed59b6e Uploaded
nml
parents:
diff changeset
690 if len(line_split) == num_gene_cols_expected+1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
691 gene_names.pop() # we read too many columns
6f870ed59b6e Uploaded
nml
parents:
diff changeset
692 num_gene_cols_expected -= 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
693 for g in gene_names:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
694 if g not in gene_names_from_fai:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
695 print "Warning: gene " + g + " in ST definitions file isn't among those in the database " + ",".join(gene_names_from_fai)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
696 print " This will result in all STs being called as unknown (but allele calls will be accurate for other loci)."
6f870ed59b6e Uploaded
nml
parents:
diff changeset
697 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
698 ST = line_split[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
699 if ST not in ST_db.values():
6f870ed59b6e Uploaded
nml
parents:
diff changeset
700 ST_string = " ".join(line_split[1:num_gene_cols_expected+1])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
701 ST_db[ST_string] = ST
6f870ed59b6e Uploaded
nml
parents:
diff changeset
702 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
703 print "Warning: this ST is not unique in the ST definitions file: " + ST
6f870ed59b6e Uploaded
nml
parents:
diff changeset
704 print "Read ST database " + ST_filename + " successfully"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
705 return (ST_db, gene_names)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
706
6f870ed59b6e Uploaded
nml
parents:
diff changeset
707 def get_allele_name_from_db(allele,unique_allele_symbols,unique_cluster_symbols,run_type,args):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
708
6f870ed59b6e Uploaded
nml
parents:
diff changeset
709 if run_type != "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
710 # header format: >[cluster]___[gene]___[allele]___[uniqueID] [info]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
711 allele_parts = allele.split()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
712 allele_detail = allele_parts.pop(0)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
713 allele_info = allele_detail.split("__")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
714
6f870ed59b6e Uploaded
nml
parents:
diff changeset
715 if len(allele_info)>2:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
716 cluster_id = allele_info[0] # ID number for the cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
717 gene_name = allele_info[1] # gene name/symbol for the cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
718 allele_name = allele_info[2] # specific allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
719 seqid = allele_info[3] # unique identifier for this seq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
720 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
721 cluster_id = gene_name = allele_name = seqid = allele_parts[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
722
6f870ed59b6e Uploaded
nml
parents:
diff changeset
723 if not unique_allele_symbols:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
724 allele_name += "_" + seqid
6f870ed59b6e Uploaded
nml
parents:
diff changeset
725
6f870ed59b6e Uploaded
nml
parents:
diff changeset
726 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
727 gene_name = allele.split(args.mlst_delimiter)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
728 allele_name = allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
729 seqid = None
6f870ed59b6e Uploaded
nml
parents:
diff changeset
730 cluster_id = None
6f870ed59b6e Uploaded
nml
parents:
diff changeset
731
6f870ed59b6e Uploaded
nml
parents:
diff changeset
732 return gene_name, allele_name, cluster_id, seqid
6f870ed59b6e Uploaded
nml
parents:
diff changeset
733
6f870ed59b6e Uploaded
nml
parents:
diff changeset
734 def create_allele_pileup(allele_name, all_pileup_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
735 outpileup = allele_name + "." + all_pileup_file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
736 with open(outpileup, 'w') as allele_pileup:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
737 with open(all_pileup_file) as all_pileup:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
738 for line in all_pileup:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
739 if line.split()[0] == allele_name:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
740 allele_pileup.write(line)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
741 return outpileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
742
6f870ed59b6e Uploaded
nml
parents:
diff changeset
743 def parse_scores(run_type,args,scores, hash_edge_depth,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
744 avg_depth_allele, coverage_allele, mismatch_allele, indel_allele,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
745 missing_allele, size_allele, next_to_del_depth_allele,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
746 unique_cluster_symbols,unique_allele_symbols, pileup_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
747
6f870ed59b6e Uploaded
nml
parents:
diff changeset
748 # sort into hash for each gene locus
6f870ed59b6e Uploaded
nml
parents:
diff changeset
749 scores_by_gene = collections.defaultdict(dict) # key1 = gene, key2 = allele, value = score
6f870ed59b6e Uploaded
nml
parents:
diff changeset
750
6f870ed59b6e Uploaded
nml
parents:
diff changeset
751 if run_type=="mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
752 for allele in scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
753 if coverage_allele[allele] > args.min_coverage:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
754 allele_info = allele.split(args.mlst_delimiter)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
755 scores_by_gene[allele_info[0]][allele] = scores[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
756 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
757 for allele in scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
758 if coverage_allele[allele] > args.min_coverage:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
759 gene_name = get_allele_name_from_db(allele,unique_allele_symbols,unique_cluster_symbols,run_type,args)[2] # cluster ID
6f870ed59b6e Uploaded
nml
parents:
diff changeset
760 scores_by_gene[gene_name][allele] = scores[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
761
6f870ed59b6e Uploaded
nml
parents:
diff changeset
762 # determine best allele for each gene locus/cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
763 results = {} # key = gene, value = (allele,diffs,depth)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
764
6f870ed59b6e Uploaded
nml
parents:
diff changeset
765 for gene in scores_by_gene:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
766
6f870ed59b6e Uploaded
nml
parents:
diff changeset
767 gene_hash = scores_by_gene[gene]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
768 scores_sorted = sorted(gene_hash.iteritems(),key=operator.itemgetter(1)) # sort by score
6f870ed59b6e Uploaded
nml
parents:
diff changeset
769 (top_allele,top_score) = scores_sorted[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
770
6f870ed59b6e Uploaded
nml
parents:
diff changeset
771 # check if depth is adequate for confident call
6f870ed59b6e Uploaded
nml
parents:
diff changeset
772 adequate_depth = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
773 depth_problem = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
774 if hash_edge_depth[top_allele][0] > args.min_edge_depth and hash_edge_depth[top_allele][1] > args.min_edge_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
775 if next_to_del_depth_allele[top_allele] != "NA":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
776 if float(next_to_del_depth_allele[top_allele]) > args.min_edge_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
777 if avg_depth_allele[top_allele] > args.min_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
778 adequate_depth = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
779 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
780 depth_problem="depth"+str(avg_depth_allele[top_allele])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
781 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
782 depth_problem = "del"+str(next_to_del_depth_allele[top_allele])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
783 elif avg_depth_allele[top_allele] > args.min_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
784 adequate_depth = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
785 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
786 depth_problem="depth"+str(avg_depth_allele[top_allele])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
787 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
788 depth_problem = "edge"+str(min(hash_edge_depth[top_allele][0],hash_edge_depth[top_allele][1]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
789
6f870ed59b6e Uploaded
nml
parents:
diff changeset
790 # check if there are confident differences against this allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
791 differences = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
792 if mismatch_allele[top_allele] > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
793 differences += str(mismatch_allele[top_allele])+"snp"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
794 if indel_allele[top_allele] > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
795 differences += str(indel_allele[top_allele])+"indel"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
796 if missing_allele[top_allele] > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
797 differences += str(missing_allele[top_allele])+"holes"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
798
6f870ed59b6e Uploaded
nml
parents:
diff changeset
799 divergence = float(mismatch_allele[top_allele]) / float( size_allele[top_allele] - missing_allele[top_allele] )
6f870ed59b6e Uploaded
nml
parents:
diff changeset
800
6f870ed59b6e Uploaded
nml
parents:
diff changeset
801 # check for truncated
6f870ed59b6e Uploaded
nml
parents:
diff changeset
802 if differences != "" or not adequate_depth:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
803 # if there are SNPs or not enough depth to trust the result, no need to screen next best match
6f870ed59b6e Uploaded
nml
parents:
diff changeset
804 results[gene] = (top_allele, differences, depth_problem, divergence)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
805 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
806 # looks good but this could be a truncated version of the real allele; check for longer versions
6f870ed59b6e Uploaded
nml
parents:
diff changeset
807 truncation_override = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
808 if len(scores_sorted) > 1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
809 (next_best_allele,next_best_score) = scores_sorted[1]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
810 if size_allele[next_best_allele] > size_allele[top_allele]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
811 # next best is longer, top allele could be a truncation?
6f870ed59b6e Uploaded
nml
parents:
diff changeset
812 if (mismatch_allele[next_best_allele] + indel_allele[next_best_allele] + missing_allele[next_best_allele]) == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
813 # next best also has no mismatches
6f870ed59b6e Uploaded
nml
parents:
diff changeset
814 if (next_best_score - top_score)/top_score < 0.1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
815 # next best has score within 10% of this one
6f870ed59b6e Uploaded
nml
parents:
diff changeset
816 truncation_override = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
817 if truncation_override:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
818 results[gene] = (next_best_allele, "trun", "", divergence) # no diffs but report this call is based on truncation test
6f870ed59b6e Uploaded
nml
parents:
diff changeset
819 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
820 results[gene] = (top_allele, "", "",divergence) # no caveats to report
6f870ed59b6e Uploaded
nml
parents:
diff changeset
821
6f870ed59b6e Uploaded
nml
parents:
diff changeset
822 # Check if there are any potential new alleles
6f870ed59b6e Uploaded
nml
parents:
diff changeset
823 #depth_problem = results[gene][2]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
824 #divergence = results[gene][3]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
825 if depth_problem == "" and divergence > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
826 new_allele = True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
827 # Get the consensus for this new allele and write it to file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
828 if args.report_new_consensus or args.report_all_consensus:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
829 new_alleles_filename = args.output + ".new_consensus_alleles.fasta"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
830 allele_pileup_file = create_allele_pileup(top_allele, pileup_file) # XXX Creates a new pileup file for that allele. Not currently cleaned up
6f870ed59b6e Uploaded
nml
parents:
diff changeset
831 read_pileup_data(allele_pileup_file, size_allele, args.prob_err, consensus_file = new_alleles_filename)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
832 if args.report_all_consensus:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
833 new_alleles_filename = args.output + ".all_consensus_alleles.fasta"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
834 allele_pileup_file = create_allele_pileup(top_allele, pileup_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
835 read_pileup_data(allele_pileup_file, size_allele, args.prob_err, consensus_file = new_alleles_filename)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
836
6f870ed59b6e Uploaded
nml
parents:
diff changeset
837 return results # (allele, diffs, depth_problem, divergence)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
838
6f870ed59b6e Uploaded
nml
parents:
diff changeset
839
6f870ed59b6e Uploaded
nml
parents:
diff changeset
840 def get_readFile_components(full_file_path):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
841 (file_path,file_name) = os.path.split(full_file_path)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
842 m1 = re.match("(.*).gz",file_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
843 ext = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
844 if m1 != None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
845 # gzipped
6f870ed59b6e Uploaded
nml
parents:
diff changeset
846 ext = ".gz"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
847 file_name = m1.groups()[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
848 (file_name_before_ext,ext2) = os.path.splitext(file_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
849 full_ext = ext2+ext
6f870ed59b6e Uploaded
nml
parents:
diff changeset
850 return(file_path,file_name_before_ext,full_ext)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
851
6f870ed59b6e Uploaded
nml
parents:
diff changeset
852 def read_file_sets(args):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
853
6f870ed59b6e Uploaded
nml
parents:
diff changeset
854 fileSets = {} # key = id, value = list of files for that sample
6f870ed59b6e Uploaded
nml
parents:
diff changeset
855 num_single_readsets = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
856 num_paired_readsets = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
857
6f870ed59b6e Uploaded
nml
parents:
diff changeset
858 if args.input_se:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
859 # single end
6f870ed59b6e Uploaded
nml
parents:
diff changeset
860 for fastq in args.input_se:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
861 (file_path,file_name_before_ext,full_ext) = get_readFile_components(fastq)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
862 m=re.match("(.*)(_S.*)(_L.*)(_R.*)(_.*)", file_name_before_ext)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
863 if m==None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
864 fileSets[file_name_before_ext] = [fastq]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
865 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
866 fileSets[m.groups()[0]] = [fastq] # Illumina names
6f870ed59b6e Uploaded
nml
parents:
diff changeset
867 num_single_readsets += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
868
6f870ed59b6e Uploaded
nml
parents:
diff changeset
869 elif args.input_pe:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
870 # paired end
6f870ed59b6e Uploaded
nml
parents:
diff changeset
871 forward_reads = {} # key = sample, value = full path to file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
872 reverse_reads = {} # key = sample, value = full path to file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
873 num_paired_readsets = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
874 num_single_readsets = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
875 for fastq in args.input_pe:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
876 (file_path,file_name_before_ext,full_ext) = get_readFile_components(fastq)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
877 # try to match to MiSeq format:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
878 m=re.match("(.*)(_S.*)(_L.*)(_R.*)(_.*)", file_name_before_ext)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
879 if m==None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
880 # not default Illumina file naming format, expect simple/ENA format
6f870ed59b6e Uploaded
nml
parents:
diff changeset
881 m=re.match("(.*)("+args.forward+")$",file_name_before_ext)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
882 if m!=None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
883 # store as forward read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
884 (baseName,read) = m.groups()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
885 forward_reads[baseName] = fastq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
886 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
887 m=re.match("(.*)("+args.reverse+")$",file_name_before_ext)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
888 if m!=None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
889 # store as reverse read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
890 (baseName,read) = m.groups()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
891 reverse_reads[baseName] = fastq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
892 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
893 logging.info("Could not determine forward/reverse read status for input file " + fastq)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
894 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
895 # matches default Illumina file naming format, e.g. m.groups() = ('samplename', '_S1', '_L001', '_R1', '_001')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
896 baseName, read = m.groups()[0], m.groups()[3]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
897 if read == "_R1":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
898 forward_reads[baseName] = fastq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
899 elif read == "_R2":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
900 reverse_reads[baseName] = fastq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
901 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
902 logging.info( "Could not determine forward/reverse read status for input file " + fastq )
6f870ed59b6e Uploaded
nml
parents:
diff changeset
903 logging.info( " this file appears to match the MiSeq file naming convention (samplename_S1_L001_[R1]_001), but we were expecting [R1] or [R2] to designate read as forward or reverse?" )
6f870ed59b6e Uploaded
nml
parents:
diff changeset
904 fileSets[file_name_before_ext] = fastq
6f870ed59b6e Uploaded
nml
parents:
diff changeset
905 num_single_readsets += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
906 # store in pairs
6f870ed59b6e Uploaded
nml
parents:
diff changeset
907 for sample in forward_reads:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
908 if sample in reverse_reads:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
909 fileSets[sample] = [forward_reads[sample],reverse_reads[sample]] # store pair
6f870ed59b6e Uploaded
nml
parents:
diff changeset
910 num_paired_readsets += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
911 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
912 fileSets[sample] = [forward_reads[sample]] # no reverse found
6f870ed59b6e Uploaded
nml
parents:
diff changeset
913 num_single_readsets += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
914 logging.info('Warning, could not find pair for read:' + forward_reads[sample])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
915 for sample in reverse_reads:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
916 if sample not in fileSets:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
917 fileSets[sample] = reverse_reads[sample] # no forward found
6f870ed59b6e Uploaded
nml
parents:
diff changeset
918 num_single_readsets += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
919 logging.info('Warning, could not find pair for read:' + reverse_reads[sample])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
920
6f870ed59b6e Uploaded
nml
parents:
diff changeset
921 if num_paired_readsets > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
922 logging.info('Total paired readsets found:' + str(num_paired_readsets))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
923 if num_single_readsets > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
924 logging.info('Total single reads found:' + str(num_single_readsets))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
925
6f870ed59b6e Uploaded
nml
parents:
diff changeset
926 return fileSets
6f870ed59b6e Uploaded
nml
parents:
diff changeset
927
6f870ed59b6e Uploaded
nml
parents:
diff changeset
928 def read_results_from_file(infile):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
929
6f870ed59b6e Uploaded
nml
parents:
diff changeset
930 if os.stat(infile).st_size == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
931 logging.info("WARNING: Results file provided is empty: " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
932 return False, False, False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
933
6f870ed59b6e Uploaded
nml
parents:
diff changeset
934 results_info = infile.split("__")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
935 if len(results_info) > 1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
936
6f870ed59b6e Uploaded
nml
parents:
diff changeset
937 if re.search("compiledResults",infile)!=None:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
938 dbtype = "compiled"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
939 dbname = results_info[0] # output identifier
6f870ed59b6e Uploaded
nml
parents:
diff changeset
940 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
941 dbtype = results_info[1] # mlst or genes
6f870ed59b6e Uploaded
nml
parents:
diff changeset
942 dbname = results_info[2] # database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
943
6f870ed59b6e Uploaded
nml
parents:
diff changeset
944 logging.info("Processing " + dbtype + " results from file " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
945
6f870ed59b6e Uploaded
nml
parents:
diff changeset
946 if dbtype == "genes":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
947 results = collections.defaultdict(dict) # key1 = sample, key2 = gene, value = allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
948 with open(infile) as f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
949 header = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
950 for line in f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
951 line_split = line.rstrip().split("\t")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
952 if len(header) == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
953 header = line_split
6f870ed59b6e Uploaded
nml
parents:
diff changeset
954 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
955 sample = line_split[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
956 for i in range(1,len(line_split)):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
957 gene = header[i] # cluster_id
6f870ed59b6e Uploaded
nml
parents:
diff changeset
958 results[sample][gene] = line_split[i]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
959
6f870ed59b6e Uploaded
nml
parents:
diff changeset
960 elif dbtype == "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
961 results = {} # key = sample, value = MLST string
6f870ed59b6e Uploaded
nml
parents:
diff changeset
962 with open(infile) as f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
963 header = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
964 for line in f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
965 if header > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
966 results[line.split("\t")[0]] = line.rstrip()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
967 if "maxMAF" not in header:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
968 results[line.split("\t")[0]] += "\tNC" # empty column for maxMAF
6f870ed59b6e Uploaded
nml
parents:
diff changeset
969 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
970 header = line.rstrip()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
971 results[line.split("\t")[0]] = line.rstrip() # store header line too (index "Sample")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
972 if "maxMAF" not in header:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
973 results[line.split("\t")[0]] += "\tmaxMAF" # add column for maxMAF
6f870ed59b6e Uploaded
nml
parents:
diff changeset
974
6f870ed59b6e Uploaded
nml
parents:
diff changeset
975 elif dbtype == "compiled":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
976 results = collections.defaultdict(dict) # key1 = sample, key2 = gene, value = allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
977 with open(infile) as f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
978 header = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
979 mlst_cols = 0 # INDEX of the last mlst column
6f870ed59b6e Uploaded
nml
parents:
diff changeset
980 n_cols = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
981 for line in f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
982 line_split = line.rstrip().split("\t")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
983 if len(header) == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
984 header = line_split
6f870ed59b6e Uploaded
nml
parents:
diff changeset
985 n_cols = len(header)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
986 if n_cols > 1:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
987 if header[1] == "ST":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
988 # there is mlst data reported
6f870ed59b6e Uploaded
nml
parents:
diff changeset
989 mlst_cols = 2 # first locus column
6f870ed59b6e Uploaded
nml
parents:
diff changeset
990 while header[mlst_cols] != "depth":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
991 mlst_cols += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
992 results["Sample"]["mlst"] = "\t".join(line_split[0:(mlst_cols+1)])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
993 results["Sample"]["mlst"] += "\tmaxMAF" # add to mlst header even if not encountered in this file, as it may be in others
6f870ed59b6e Uploaded
nml
parents:
diff changeset
994 if header[mlst_cols+1] == "maxMAF":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
995 mlst_cols += 1 # record maxMAF column within MLST data, if present
6f870ed59b6e Uploaded
nml
parents:
diff changeset
996 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
997 # no mlst data reported
6f870ed59b6e Uploaded
nml
parents:
diff changeset
998 dbtype = "genes"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
999 logging.info("No MLST data in compiled results file " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1000 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1001 # no mlst data reported
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1002 dbtype = "genes"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1003 logging.info("No MLST data in compiled results file " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1004
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1005 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1006 sample = line_split[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1007 if mlst_cols > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1008 results[sample]["mlst"] = "\t".join(line_split[0:(mlst_cols+1)])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1009 if "maxMAF" not in header:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1010 results[sample]["mlst"] += "\t" # add to mlst section even if not encountered in this file, as it may be in others
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1011 if n_cols > mlst_cols:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1012 # read genes component
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1013 for i in range(mlst_cols+1,n_cols):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1014 # note i=1 if mlst_cols==0, ie we are reading all
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1015 gene = header[i]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1016 if len(line_split) > i:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1017 results[sample][gene] = line_split[i]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1018 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1019 results[sample][gene] = "-"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1020 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1021 results = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1022 dbtype = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1023 dbname = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1024 logging.info("Couldn't decide what to do with file results file provided: " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1025
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1026 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1027 results = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1028 dbtype = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1029 dbname = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1030 logging.info("Couldn't decide what to do with file results file provided: " + infile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1031
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1032 return results, dbtype, dbname
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1033
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1034 def read_scores_file(scores_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1035 hash_edge_depth = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1036 avg_depth_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1037 coverage_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1038 mismatch_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1039 indel_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1040 missing_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1041 size_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1042 next_to_del_depth_allele = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1043 mix_rates = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1044 scores = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1045
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1046 f = file(scores_file,"r")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1047
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1048 for line in f:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1049 line_split = line.rstrip().split("\t")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1050 allele = line_split[0]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1051 if allele != "Allele": # skip header row
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1052 scores[allele] = float(line_split[1])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1053 mix_rates[allele] = float(line_split[11])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1054 avg_depth_allele[allele] = float(line_split[2])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1055 hash_edge_depth[allele] = (float(line_split[3]),float(line_split[4]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1056 coverage_allele[allele] = float(line_split[5])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1057 size_allele[allele] = int(line_split[6])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1058 mismatch_allele[allele] = int(line_split[7])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1059 indel_allele[allele] = int(line_split[8])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1060 missing_allele[allele] = int(line_split[9])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1061 next_to_del_depth = line_split[10]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1062 next_to_del_depth_allele[allele] = line_split[10]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1063
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1064 return hash_edge_depth, avg_depth_allele, coverage_allele, mismatch_allele, indel_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1065 missing_allele, size_allele, next_to_del_depth_allele, scores, mix_rates
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1066
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1067 def run_srst2(args, fileSets, dbs, run_type):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1068
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1069 db_reports = [] # list of db-specific output files to return
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1070 db_results_list = [] # list of results hashes, one per db
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1071
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1072 for fasta in dbs:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1073 db_reports, db_results_list = process_fasta_db(args, fileSets, run_type, db_reports, db_results_list, fasta)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1074
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1075 return db_reports, db_results_list
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1076
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1077 def process_fasta_db(args, fileSets, run_type, db_reports, db_results_list, fasta):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1078
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1079 check_command_version(['samtools'],
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1080 'Version: 0.1.18',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1081 'samtools',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1082 '0.1.18')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1083
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1084 logging.info('Processing database ' + fasta)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1085
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1086 db_path, db_name = os.path.split(fasta) # database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1087 (db_name,db_ext) = os.path.splitext(db_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1088 db_results = "__".join([args.output,run_type,db_name,"results.txt"])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1089 db_report = file(db_results,"w")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1090 db_reports.append(db_results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1091
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1092 # Get sequence lengths and gene names
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1093 # lengths are needed for MLST heuristic to distinguish alleles from their truncated forms
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1094 # gene names read from here are needed for non-MLST dbs
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1095 fai_file = fasta + '.fai'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1096 if not os.path.exists(fai_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1097 run_command(['samtools', 'faidx', fasta])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1098 size, gene_names, unique_gene_symbols, unique_allele_symbols, cluster_symbols = \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1099 parse_fai(fai_file,run_type,args.mlst_delimiter)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1100
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1101 # Prepare for MLST reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1102 ST_db = False
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1103 if run_type == "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1104 results = {} # key = sample, value = ST string for printing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1105 if args.mlst_definitions:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1106 # store MLST profiles, replace gene names (we want the order as they appear in this file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1107 ST_db, gene_names = parse_ST_database(args.mlst_definitions,gene_names)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1108 db_report.write("\t".join(["Sample","ST"]+gene_names+["mismatches","uncertainty","depth","maxMAF"]) + "\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1109 results["Sample"] = "\t".join(["Sample","ST"]+gene_names+["mismatches","uncertainty","depth","maxMAF"])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1110
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1111 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1112 # store final results for later tabulation
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1113 results = collections.defaultdict(dict) #key1 = sample, key2 = gene, value = allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1114
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1115 gene_list = [] # start with empty gene list; will add genes from each genedb test
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1116
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1117 # determine maximum mismatches per read to use for pileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1118 if run_type == "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1119 max_mismatch = args.mlst_max_mismatch
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1120 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1121 max_mismatch = args.gene_max_mismatch
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1122
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1123 # Align and score each read set against this DB
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1124 for sample_name in fileSets:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1125 logging.info('Processing sample ' + sample_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1126 fastq_inputs = fileSets[sample_name] # reads
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1127
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1128 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1129 # try mapping and scoring this fileset against the current database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1130 # update the gene_list list and results dict with data from this strain
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1131 # __mlst__ will be printed during this routine if this is a mlst run
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1132 # __fullgenes__ will be printed during this routine if requested and this is a gene_db run
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1133 gene_list, results = \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1134 map_fileSet_to_db(args,sample_name,fastq_inputs,db_name,fasta,size,gene_names,\
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1135 unique_gene_symbols, unique_allele_symbols,run_type,ST_db,results,gene_list,db_report,cluster_symbols,max_mismatch)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1136 # if we get an error from one of the commands we called
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1137 # log the error message and continue onto the next fasta db
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1138 except CommandError as e:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1139 logging.error(e.message)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1140 # record results as unknown, so we know that we did attempt to analyse this readset
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1141 if run_type == "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1142 st_result_string = "\t".join( [sample_name,"-"] + ["-"] * (len(gene_names)+3)) # record missing results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1143 db_report.write( st_result_string + "\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1144 logging.info(" " + st_result_string)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1145 results[sample_name] = st_result_string
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1146 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1147 logging.info(" failed gene detection")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1148 results[sample_name]["failed"] = True # so we know that we tried this strain
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1149
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1150 if run_type != "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1151 # tabulate results across samples for this gene db (i.e. __genes__ file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1152 logging.info('Tabulating results for database {} ...'.format(fasta))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1153 gene_list.sort()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1154 db_report.write("\t".join(["Sample"]+gene_list)+"\n") # report header row
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1155 for sample_name in fileSets:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1156 db_report.write(sample_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1157 if sample_name in results:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1158 # print results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1159 if "failed" not in results[sample_name]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1160 for cluster_id in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1161 if cluster_id in results[sample_name]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1162 db_report.write("\t"+results[sample_name][cluster_id]) # print full allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1163 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1164 db_report.write("\t-") # no hits for this gene cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1165 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1166 # no data on this, as the sample failed mapping
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1167 for cluster_id in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1168 db_report.write("\t?") #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1169 results[sample_name][cluster_id] = "?" # record as unknown
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1170 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1171 # no data on this because genes were not found (but no mapping errors)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1172 for cluster_id in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1173 db_report.write("\t?") #
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1174 results[sample_name][cluster_id] = "-" # record as absent
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1175 db_report.write("\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1176
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1177 # Finished with this database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1178 logging.info('Finished processing for database {} ...'.format(fasta))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1179 db_report.close()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1180 db_results_list.append(results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1181
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1182 return db_reports, db_results_list
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1183
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1184 def map_fileSet_to_db(args,sample_name,fastq_inputs,db_name,fasta,size,gene_names,\
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1185 unique_gene_symbols, unique_allele_symbols,run_type,ST_db,results,gene_list,db_report,cluster_symbols,max_mismatch):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1186
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1187 mapping_files_pre = args.output + '__' + sample_name + '.' + db_name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1188 pileup_file = mapping_files_pre + '.pileup'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1189 scores_file = mapping_files_pre + '.scores'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1190
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1191 # Get or read scores
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1192
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1193 if args.use_existing_scores and os.path.exists(scores_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1194
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1195 logging.info(' Using existing scores in ' + scores_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1196
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1197 # read in scores and info from existing scores file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1198 hash_edge_depth, avg_depth_allele, coverage_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1199 mismatch_allele, indel_allele, missing_allele, size_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1200 next_to_del_depth_allele, scores, mix_rates = read_scores_file(scores_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1201
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1202 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1203
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1204 # Get or read pileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1205
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1206 if args.use_existing_pileup and os.path.exists(pileup_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1207 logging.info(' Using existing pileup in ' + pileup_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1208
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1209 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1210
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1211 # run bowtie against this db
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1212 bowtie_sam = run_bowtie(mapping_files_pre,sample_name,fastq_inputs,args,db_name,fasta)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1213
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1214 # Modify Bowtie's SAM formatted output so that we get secondary
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1215 # alignments in downstream pileup
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1216 (raw_bowtie_sam,bowtie_sam_mod) = modify_bowtie_sam(bowtie_sam,max_mismatch)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1217
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1218 # generate pileup from sam (via sorted bam)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1219 get_pileup(args,mapping_files_pre,raw_bowtie_sam,bowtie_sam_mod,fasta,pileup_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1220
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1221 # Get scores
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1222
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1223 # Process the pileup and extract info for scoring and reporting on each allele
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1224 logging.info(' Processing SAMtools pileup...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1225 hash_alignment, hash_max_depth, hash_edge_depth, avg_depth_allele, coverage_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1226 mismatch_allele, indel_allele, missing_allele, size_allele, next_to_del_depth_allele= \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1227 read_pileup_data(pileup_file, size, args.prob_err)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1228
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1229 # Generate scores for all alleles (prints these and associated info if verbose)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1230 # result = dict, with key=allele, value=score
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1231 logging.info(' Scoring alleles...')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1232 scores, mix_rates = score_alleles(args, mapping_files_pre, hash_alignment, hash_max_depth, hash_edge_depth, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1233 avg_depth_allele, coverage_allele, mismatch_allele, indel_allele, missing_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1234 size_allele, next_to_del_depth_allele, run_type)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1235
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1236 # GET BEST SCORE for each gene/cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1237 # result = dict, with key = gene, value = (allele,diffs,depth_problem)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1238 # for MLST DBs, key = gene = locus, allele = gene-number
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1239 # for gene DBs, key = gene = cluster ID, allele = cluster__gene__allele__id
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1240 # for gene DBs, only those alleles passing the coverage cutoff are returned
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1241
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1242 allele_scores = parse_scores(run_type, args, scores, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1243 hash_edge_depth, avg_depth_allele, coverage_allele, mismatch_allele, \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1244 indel_allele, missing_allele, size_allele, next_to_del_depth_allele,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1245 unique_gene_symbols, unique_allele_symbols, pileup_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1246
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1247 # REPORT/RECORD RESULTS
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1248
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1249 # Report MLST results to __mlst__ file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1250 if run_type == "mlst" and len(allele_scores) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1251
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1252 # Calculate ST and get info for reporting
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1253 (st,clean_st,alleles_with_flags,mismatch_flags,uncertainty_flags,mean_depth,max_maf) = \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1254 calculate_ST(allele_scores, ST_db, gene_names, sample_name, args.mlst_delimiter, avg_depth_allele, mix_rates)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1255
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1256 # Print to MLST report, log and save the result
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1257 st_result_string = "\t".join([sample_name,st]+alleles_with_flags+[";".join(mismatch_flags),";".join(uncertainty_flags),str(mean_depth),str(max_maf)])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1258 db_report.write( st_result_string + "\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1259 logging.info(" " + st_result_string)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1260 results[sample_name] = st_result_string
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1261
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1262 # Make sure scores are printed if there was uncertainty in the call
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1263 scores_output_file = mapping_files_pre + '.scores'
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1264 if uncertainty_flags != ["-"] and not args.save_scores and not os.path.exists(scores_output_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1265 # print full score set
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1266 logging.info("Printing all MLST scores to " + scores_output_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1267 scores_output = file(scores_output_file, 'w')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1268 scores_output.write("Allele\tScore\tAvg_depth\tEdge1_depth\tEdge2_depth\tPercent_coverage\tSize\tMismatches\tIndels\tTruncated_bases\tDepthNeighbouringTruncation\tMmaxMAF\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1269 for allele in scores.keys():
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1270 score = scores[allele]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1271 scores_output.write('\t'.join([allele, str(score), str(avg_depth_allele[allele]), \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1272 str(hash_edge_depth[allele][0]), str(hash_edge_depth[allele][1]), \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1273 str(coverage_allele[allele]), str(size_allele[allele]), str(mismatch_allele[allele]), \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1274 str(indel_allele[allele]), str(missing_allele[allele]), str(next_to_del_depth_allele[allele]), str(round(mix_rates[allele],3))]) + '\n')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1275 scores_output.close()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1276
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1277 # Record gene results for later processing and optionally print detailed gene results to __fullgenes__ file
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1278 elif run_type == "genes" and len(allele_scores) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1279 if args.no_gene_details:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1280 full_results = "__".join([args.output,"full"+run_type,db_name,"results.txt"])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1281 logging.info("Printing verbose gene detection results to " + full_results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1282 f = file(full_results,"w")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1283 f.write("\t".join(["Sample","DB","gene","allele","coverage","depth","diffs","uncertainty","divergence","length", "maxMAF","clusterid","seqid","annotation"])+"\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1284 for gene in allele_scores:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1285 (allele,diffs,depth_problem,divergence) = allele_scores[gene] # gene = top scoring alleles for each cluster
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1286 gene_name, allele_name, cluster_id, seqid = \
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1287 get_allele_name_from_db(allele,unique_allele_symbols,unique_gene_symbols,run_type,args)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1288
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1289 # store for gene result table only if divergence passes minimum threshold:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1290 if divergence*100 <= float(args.max_divergence):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1291 column_header = cluster_symbols[cluster_id]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1292 results[sample_name][column_header] = allele_name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1293 if diffs != "":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1294 results[sample_name][column_header] += "*"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1295 if depth_problem != "":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1296 results[sample_name][column_header] += "?"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1297 if gene not in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1298 gene_list.append(column_header)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1299
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1300 # write details to full genes report
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1301 if args.no_gene_details:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1302
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1303 # get annotation info
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1304 header_string = os.popen(" ".join(["grep",allele,fasta]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1305 try:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1306 header = header_string.read().rstrip().split()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1307 header.pop(0) # remove allele name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1308 if len(header) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1309 annotation = " ".join(header) # put back the spaces
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1310 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1311 annotation = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1312
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1313 except:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1314 annotation = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1315
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1316 f.write("\t".join([sample_name,db_name,gene_name,allele_name,str(round(coverage_allele[allele],3)),str(avg_depth_allele[allele]),diffs,depth_problem,str(round(divergence*100,3)),str(size_allele[allele]),str(round(mix_rates[allele],3)),cluster_id,seqid,annotation])+"\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1317
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1318 # log the gene detection result
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1319 logging.info(" " + str(len(allele_scores)) + " genes identified in " + sample_name)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1320
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1321 # Finished with this read set
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1322 logging.info(' Finished processing for read set {} ...'.format(sample_name))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1323
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1324 return gene_list, results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1325
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1326 def compile_results(args,mlst_results,db_results,compiled_output_file):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1327
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1328 o = file(compiled_output_file,"w")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1329
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1330 # get list of all samples and genes present in these datasets
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1331 sample_list = [] # each entry is a sample present in at least one db
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1332 gene_list = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1333 mlst_cols = 0
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1334 mlst_header_string = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1335 blank_mlst_section = ""
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1336
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1337 mlst_results_master = {} # compilation of all MLST results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1338 db_results_master = collections.defaultdict(dict) # compilation of all gene results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1339 st_counts = {} # key = ST, value = count
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1340
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1341 if len(mlst_results) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1342
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1343 for mlst_result in mlst_results:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1344
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1345 # check length of the mlst string
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1346 if "Sample" in mlst_result:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1347 test_string = mlst_result["Sample"]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1348 if mlst_cols == 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1349 mlst_header_string = test_string
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1350 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1351 test_string = mlst_result[mlst_result.keys()[0]] # no header line?
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1352 test_string_split = test_string.split("\t")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1353 this_mlst_cols = len(test_string)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1354
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1355 if (mlst_cols == 0) or (mlst_cols == this_mlst_cols):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1356 mlst_cols = this_mlst_cols
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1357 blank_mlst_section = "\t" * (mlst_cols-1) # blank MLST string in case some samples missing
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1358 # use this data
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1359 for sample in mlst_result:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1360 mlst_results_master[sample] = mlst_result[sample]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1361 if sample not in sample_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1362 sample_list.append(sample)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1363 elif mlst_cols != this_mlst_cols:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1364 # don't process this data further
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1365 logging.info("Problem reconciling MLST data from two files, first MLST results encountered had " + str(mlst_cols) + " columns, this one has " + str(this_mlst_cols) + " columns?")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1366 if args.mlst_db:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1367 logging.info("Compiled report will contain only the MLST data from this run, not previous outputs")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1368 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1369 logging.info("Compiled report will contain only the data from the first MLST result set provided")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1370
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1371 if len(db_results) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1372 for results in db_results:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1373 for sample in results:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1374 if sample not in sample_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1375 sample_list.append(sample)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1376 for gene in results[sample]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1377 if gene != "failed":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1378 db_results_master[sample][gene] = results[sample][gene]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1379 if gene not in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1380 gene_list.append(gene)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1381
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1382 if "Sample" in sample_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1383 sample_list.remove("Sample")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1384 sample_list.sort()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1385 gene_list.sort()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1386
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1387 # print header
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1388 header_elements = []
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1389 if len(mlst_results) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1390 header_elements.append(mlst_header_string)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1391 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1392 header_elements.append("Sample")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1393 if (gene_list) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1394 header_elements += gene_list
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1395 o.write("\t".join(header_elements)+"\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1396
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1397 # print results for all samples
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1398 for sample in sample_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1399
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1400 sample_info = [] # first entry is mlst string OR sample name, rest are genes
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1401
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1402 # print mlst if provided, otherwise just print sample name
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1403 if len(mlst_results_master) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1404 if sample in mlst_results_master:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1405 sample_info.append(mlst_results_master[sample])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1406 this_st = mlst_results_master[sample].split("\t")[1]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1407 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1408 sample_info.append(sample+blank_mlst_section)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1409 this_st = "unknown"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1410 # record the MLST result
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1411 if this_st in st_counts:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1412 st_counts[this_st] += 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1413 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1414 st_counts[this_st] = 1
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1415 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1416 sample_info.append(sample)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1417
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1418 # get gene info if provided
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1419 if sample in db_results_master:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1420 for gene in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1421 if gene in db_results_master[sample]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1422 sample_info.append(db_results_master[sample][gene])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1423 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1424 sample_info.append("-")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1425 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1426 for gene in gene_list:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1427 sample_info.append("?") # record no gene data on this strain
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1428
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1429 o.write("\t".join(sample_info)+"\n")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1430
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1431 o.close()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1432
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1433 logging.info("Compiled data on " + str(len(sample_list)) + " samples printed to: " + compiled_output_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1434
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1435 # log ST counts
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1436 if len(mlst_results_master) > 0:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1437 logging.info("Detected " + str(len(st_counts.keys())) + " STs: ")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1438 sts = st_counts.keys()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1439 sts.sort()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1440 for st in sts:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1441 logging.info("ST" + st + "\t" + str(st_counts[st]))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1442
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1443 return True
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1444
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1445
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1446 def main():
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1447 args = parse_args()
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1448 if args.log is True:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1449 logfile = args.output + ".log"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1450 else:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1451 logfile = None
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1452 logging.basicConfig(
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1453 filename=logfile,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1454 level=logging.DEBUG,
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1455 filemode='w',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1456 format='%(asctime)s %(message)s',
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1457 datefmt='%m/%d/%Y %H:%M:%S')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1458 logging.info('program started')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1459 logging.info('command line: {0}'.format(' '.join(sys.argv)))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1460
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1461 # Delete consensus file if it already exists (so can use append file in funtions)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1462 if args.report_new_consensus or args.report_all_consensus:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1463 new_alleles_filename = args.output + ".consensus_alleles.fasta"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1464 if os.path.exists(new_alleles_filename):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1465 os.remove(new_alleles_filename)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1466
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1467 # vars to store results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1468 mlst_results_hashes = [] # dict (sample->MLST result string) for each MLST output files created/read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1469 gene_result_hashes = [] # dict (sample->gene->result) for each gene typing output files created/read
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1470
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1471 # parse list of file sets to analyse
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1472 fileSets = read_file_sets(args) # get list of files to process
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1473
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1474 # run MLST scoring
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1475 if fileSets and args.mlst_db:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1476
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1477 if not args.mlst_definitions:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1478
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1479 # print warning to screen to alert user, may want to stop and restart
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1480 print "Warning, MLST allele sequences were provided without ST definitions:"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1481 print " allele sequences: " + str(args.mlst_db)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1482 print " these will be mapped and scored, but STs can not be calculated"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1483
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1484 # log
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1485 logging.info("Warning, MLST allele sequences were provided without ST definitions:")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1486 logging.info(" allele sequences: " + str(args.mlst_db))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1487 logging.info(" these will be mapped and scored, but STs can not be calculated")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1488
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1489 bowtie_index(args.mlst_db) # index the MLST database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1490
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1491 # score file sets against MLST database
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1492 mlst_report, mlst_results = run_srst2(args,fileSets,args.mlst_db,"mlst")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1493
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1494 logging.info('MLST output printed to ' + mlst_report[0])
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1495
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1496 #mlst_reports_files += mlst_report
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1497 mlst_results_hashes += mlst_results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1498
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1499 # run gene detection
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1500 if fileSets and args.gene_db:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1501
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1502 bowtie_index(args.gene_db) # index the gene databases
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1503
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1504 db_reports, db_results = run_srst2(args,fileSets,args.gene_db,"genes")
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1505
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1506 for outfile in db_reports:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1507 logging.info('Gene detection output printed to ' + outfile)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1508
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1509 gene_result_hashes += db_results
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1510
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1511 # process prior results files
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1512 if args.prev_output:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1513
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1514 unique_results_files = list(OrderedDict.fromkeys(args.prev_output))
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1515
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1516 for results_file in unique_results_files:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1517
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1518 results, dbtype, dbname = read_results_from_file(results_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1519
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1520 if dbtype == "mlst":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1521 mlst_results_hashes.append(results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1522
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1523 elif dbtype == "genes":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1524 gene_result_hashes.append(results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1525
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1526 elif dbtype == "compiled":
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1527 # store mlst in its own db
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1528 mlst_results = {}
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1529 for sample in results:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1530 if "mlst" in results[sample]:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1531 mlst_results[sample] = results[sample]["mlst"]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1532 del results[sample]["mlst"]
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1533 mlst_results_hashes.append(mlst_results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1534 gene_result_hashes.append(results)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1535
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1536 # compile results if multiple databases or datasets provided
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1537 if ( (len(gene_result_hashes) + len(mlst_results_hashes)) > 1 ):
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1538 compiled_output_file = args.output + "__compiledResults.txt"
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1539 compile_results(args,mlst_results_hashes,gene_result_hashes,compiled_output_file)
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1540
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1541 elif args.prev_output:
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1542 logging.info('One previous output file was provided, but there is no other data to compile with.')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1543
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1544 logging.info('SRST2 has finished.')
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1545
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1546
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1547 if __name__ == '__main__':
6f870ed59b6e Uploaded
nml
parents:
diff changeset
1548 main()