annotate allele-counts.py @ 11:cf2af5c3118c draft default tip

"planemo upload for repository https://github.com/galaxyproject/dunovo commit ac9577f0047ad984b307fe2c4b40e2eb45a0e6e2"
author nick
date Fri, 03 Apr 2020 16:05:04 -0400
parents 6cc488e11544
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
1 #!/usr/bin/python3
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
2 """
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
3 Run with -h option or see DESCRIPTION for description.
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
4 This script's functionality is being obsoleted by the new, and much more sanely
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
5 written, nvc-filter.py.
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
7 New in this version:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
8 - Calculate strand bias
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
9 - Rename MINOR.FREQ.PERC to MAF
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
10
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
11 Naive Variant Caller variant count parsing one-liner:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
12 $ cat variants.vcf | grep -v '^#' | cut -f 10 | cut -d ':' -f 4 | tr ',=' '\t:'
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
13 """
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
14 import os
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
15 import sys
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
16 import errno
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
17 import random
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
18 from optparse import OptionParser
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
19
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
20 COLUMNS = ['sample', 'chr', 'pos', 'A', 'C', 'G', 'T', 'coverage', 'alleles',
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
21 'major', 'minor', 'freq', 'bias']
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
22 COLUMN_LABELS = ['SAMPLE', 'CHR', 'POS', 'A', 'C', 'G', 'T', 'CVRG', 'ALLELES',
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
23 'MAJOR', 'MINOR', 'MAF', 'BIAS']
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
24 CANONICAL_VARIANTS = ['A', 'C', 'G', 'T']
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
25 USAGE = """Usage: %prog [options] -i variants.vcf -o alleles.tsv
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
26 cat variants.vcf | %prog [options] > alleles.tsv"""
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
27 OPT_DEFAULTS = {'infile':'-', 'outfile':'-', 'freq_thres':1.0, 'covg_thres':100,
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
28 'print_header':False, 'stdin':False, 'stranded':False, 'no_filter':False,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
29 'debug_loc':'', 'seed':''}
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
30 DESCRIPTION = """This will parse the VCF output of the "Naive Variant Caller"
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
31 (aka "BAM Coverage") Galaxy tool. For each position reported, it counts the
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
32 number of reads of each base, determines the major allele, minor allele (second
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
33 most frequent variant), and number of alleles above a threshold. So currently
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
34 it only considers SNVs (ACGT), including in the coverage figure. By default it
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
35 reads from stdin and prints to stdout.
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
36 Prints a tab-delimited set of statistics to stdout.
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
37 To print output column labels, run "$ echo -n | ./allele-counts.py -H".
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
38 The columns are: 1:SAMPLE 2:CHR 3:POS 4:A 5:C 6:G 7:T 8:CVRG 9:ALLELES 10:MAJOR
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
39 11:MINOR 12:MAF 13:BIAS,
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
40 unless the --stranded option is used, in which case they are:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
41 1:SAMPLE 2:CHR 3:POS 4:+A 5:+C 6:+G 7:+T 8:-A 9:-C 10:-G 11:-T 12:CVRG
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
42 13:ALLELES 14:MAJOR 15:MINOR 16:MAF 17:BIAS.
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
43 """
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
44 EPILOG = """Requirements:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
45 The input VCF must report the variants for each strand.
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
46 The variants should be case-sensitive (e.g. all capital base letters).
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
47 Strand bias: Both strands must show the same bases passing the frequency
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
48 threshold (but not necessarily in the same order). If the site fails this test,
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
49 the number of alleles is reported as 0."""
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
50
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
51
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
52 def get_options(defaults, usage, description='', epilog=''):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
53 """Get options, print usage text."""
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
54
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
55 parser = OptionParser(usage=usage, description=description, epilog=epilog)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
56
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
57 parser.add_option('-i', '--infile', dest='infile',
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
58 default=defaults.get('infile'),
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
59 help='Read input VCF data from this file instead of stdin.')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
60 parser.add_option('-o', '--outfile', dest='outfile',
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
61 default=defaults.get('outfile'),
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
62 help='Print output data to this file instead of stdout.')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
63 parser.add_option('-f', '--freq-thres', dest='freq_thres', type='float',
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
64 default=defaults.get('freq_thres'),
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
65 help=('Frequency threshold for counting alleles, given in percentage: -f 1 '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
66 +'= 1% frequency. Default is %default%.'))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
67 parser.add_option('-c', '--covg-thres', dest='covg_thres', type='int',
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
68 default=defaults.get('covg_thres'),
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
69 help=('Coverage threshold. Each site must be supported by at least this '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
70 +'many reads on each strand. Otherwise the site will not be printed in '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
71 +'the output. The default is %default reads per strand.'))
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
72 parser.add_option('-n', '--no-filter', dest='no_filter', action='store_const',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
73 const=not(defaults.get('no_filter')), default=defaults.get('no_filter'),
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
74 help=('Operate without a frequency threshold or coverage threshold. '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
75 +'Equivalent to "-c 0 -f 0".'))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
76 parser.add_option('-H', '--header', dest='print_header', action='store_const',
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
77 const=not(defaults.get('print_header')), default=defaults.get('print_header'),
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
78 help=('Print header line. This is a #-commented line with the column '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
79 +'labels. Off by default.'))
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
80 parser.add_option('-s', '--stranded', dest='stranded', action='store_const',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
81 const=not(defaults.get('stranded')), default=defaults.get('stranded'),
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
82 help='Report variant counts by strand, in separate columns. Off by default.')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
83 parser.add_option('-r', '--rand-seed', dest='seed',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
84 default=defaults.get('seed'), help=('Seed for random number generator.'))
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
85 parser.add_option('-d', '--debug', dest='debug_loc',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
86 default=defaults.get('debug_loc'),
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
87 help=('Turn on debug mode and specify a single site to process using UCSC '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
88 +'coordinate format. You can also append a sample ID after another ":" '
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
89 +'to restrict it further.'))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
90
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
91 (options, args) = parser.parse_args()
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
92
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
93 return (options, args)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
94
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
95
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
96 def main():
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
97
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
98 (options, args) = get_options(OPT_DEFAULTS, USAGE, DESCRIPTION, EPILOG)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
99
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
100 infile = options.infile
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
101 outfile = options.outfile
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
102 print_header = options.print_header
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
103 freq_thres = options.freq_thres / 100
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
104 covg_thres = options.covg_thres
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
105 stranded = options.stranded
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
106 debug_loc = options.debug_loc
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
107 seed = options.seed
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
108
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
109 if options.no_filter:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
110 freq_thres = 0
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
111 covg_thres = 0
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
112
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
113 if seed:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
114 random.seed(seed)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
115
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
116 debug = False
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
117 print_sample = ''
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
118 if debug_loc:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
119 debug = True
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
120 coords = debug_loc.split(':')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
121 print_chr = coords[0]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
122 print_pos = ''
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
123 if len(coords) > 1: print_pos = coords[1]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
124 if len(coords) > 2: print_sample = coords[2]
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
125
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
126 # set infile_handle to either stdin or the input file
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
127 if infile == OPT_DEFAULTS.get('infile'):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
128 infile_handle = sys.stdin
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
129 sys.stderr.write("Reading from standard input..\n")
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
130 else:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
131 if os.path.exists(infile):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
132 infile_handle = open(infile, 'r')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
133 else:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
134 fail('Error: Input VCF file '+infile+' not found.')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
135
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
136 # set outfile_handle to either stdout or the output file
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
137 if outfile == OPT_DEFAULTS.get('outfile'):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
138 outfile_handle = sys.stdout
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
139 else:
2
318fdf77aa54 Current version, from another toolshed
nick
parents: 1
diff changeset
140 try:
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
141 outfile_handle = open(outfile, 'w')
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
142 except IOError:
2
318fdf77aa54 Current version, from another toolshed
nick
parents: 1
diff changeset
143 fail('Error: The given output filename '+outfile+' could not be opened.')
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
144
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
145 # Take care of column names, print header
2
318fdf77aa54 Current version, from another toolshed
nick
parents: 1
diff changeset
146 if len(COLUMNS) != len(COLUMN_LABELS):
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
147 fail('Error: Internal column names list do not match column labels list.')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
148 if stranded:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
149 COLUMNS[3:7] = ['+A', '+C', '+G', '+T', '-A', '-C', '-G', '-T']
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
150 COLUMN_LABELS[3:7] = ['+A', '+C', '+G', '+T', '-A', '-C', '-G', '-T']
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
151 if print_header:
2
318fdf77aa54 Current version, from another toolshed
nick
parents: 1
diff changeset
152 outfile_handle.write('#'+'\t'.join(COLUMN_LABELS)+"\n")
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
153
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
154 # main loop
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
155 # each iteration processes one VCF line and prints one or more output lines
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
156 # one VCF line = one site, one or more samples
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
157 # one output line = one site, one sample
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
158 sample_names = []
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
159 for line in infile_handle:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
160 line = line.rstrip('\r\n')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
161
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
162 # header lines
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
163 if line[0] == '#':
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
164 if line[0:6].upper() == '#CHROM':
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
165 sample_names = line.split('\t')[9:]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
166 continue
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
167
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
168 if not sample_names:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
169 fail("Error in input VCF: Data line encountered before header line. "
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
170 +"Failed on line:\n"+line)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
171
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
172 site_data = read_site(line, sample_names, CANONICAL_VARIANTS)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
173
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
174 if debug:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
175 if print_pos != site_data['pos']:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
176 continue
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
177 if print_chr != site_data['chr'] and print_chr != '':
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
178 continue
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
179 if print_sample != '':
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
180 for sample in site_data['samples'].keys():
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
181 if sample.lower() != print_sample.lower():
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
182 site_data['samples'].pop(sample, None)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
183 if len(site_data['samples']) == 0:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
184 sys.stderr.write("Error: Sample '"+print_sample+"' not found.\n")
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
185 sys.exit(1)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
186
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
187 site_summary = summarize_site(site_data, sample_names, CANONICAL_VARIANTS,
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
188 freq_thres, covg_thres, stranded, debug=debug)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
189
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
190 if debug and site_summary[0]['print']:
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
191 print(line.split('\t')[9].split(':')[-1])
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
192
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
193 try:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
194 print_site(outfile_handle, site_summary, COLUMNS)
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
195 except IOError as ioe:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
196 if ioe.errno == errno.EPIPE:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
197 sys.exit(0)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
198
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
199 # keeps Galaxy from giving an error if there were messages on stderr
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
200 sys.exit(0)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
201
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
202
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
203
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
204 def read_site(line, sample_names, canonical):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
205 """Read in a line, parse the variants into a data structure, and return it.
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
206 The line should be actual site data, not a header line, so check beforehand.
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
207 Only the variants in 'canonical' will be read; all others are ignored.
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
208 Note: the line is assumed to have been chomped.
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
209 The returned data is stored in a dict, with the following structure:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
210 {
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
211 'chr': 'chr1',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
212 'pos': '2617',
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
213 'samples': {
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
214 'THYROID': {
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
215 '+A': 32,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
216 '-A': 45,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
217 '-G': 1,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
218 },
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
219 'BLOOD': {
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
220 '+A': 2,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
221 '-C': 1,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
222 '+G': 37,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
223 '-G': 42,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
224 },
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
225 },
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
226 }
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
227 """
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
228
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
229 site = {}
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
230 fields = line.split('\t')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
231
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
232 if len(fields) < 9:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
233 fail("Error in input VCF: wrong number of fields in data line. "
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
234 "Failed on line:\n"+line)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
235
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
236 site['chr'] = fields[0]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
237 site['pos'] = fields[1]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
238 samples = fields[9:]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
239
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
240 if len(samples) < len(sample_names):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
241 fail("Error in input VCF: missing sample fields in data line. "
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
242 "Failed on line:\n"+line)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
243 elif len(samples) > len(sample_names):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
244 fail("Error in input VCF: more sample fields in data line than in header. "
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
245 "Failed on line:\n"+line)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
246
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
247 sample_counts = {}
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
248 for i in range(len(samples)):
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
249
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
250 variant_counts = {}
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
251 counts = samples[i].split(':')[-1]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
252 counts = counts.split(',')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
253
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
254 for count in counts:
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
255 if not count or count == '.':
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
256 continue
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
257 fields = count.split('=')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
258 if len(fields) != 2:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
259 fail("Error in input VCF: Incorrect variant data format (must contain "
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
260 "a single '='). Failed on data \"{}\" in line:\n{}"
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
261 .format(count, line))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
262 (variant, reads) = fields
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
263 if variant[1:] not in canonical:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
264 continue
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
265 if not variant.startswith('-') and not variant.startswith('+'):
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
266 fail("Error in input VCF: variant data not strand-specific. Failed on "
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
267 "data \"{}\" on line:\n{}".format(variant, line))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
268 try:
2
318fdf77aa54 Current version, from another toolshed
nick
parents: 1
diff changeset
269 variant_counts[variant] = int(float(reads))
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
270 except ValueError:
8
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
271 fail("Error in input VCF: Variant count not a valid number. Failed on "
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
272 "variant count string \"{}\"\nIn the following line:\n{}"
411adeff1eec Handle &#34;.&#34; sample columns, update tests to work with BIAS column.
nick
parents: 6
diff changeset
273 .format(reads, line))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
274
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
275 sample_counts[sample_names[i]] = variant_counts
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
276
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
277 site['samples'] = sample_counts
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
278
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
279 return site
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
280
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
281
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
282 def summarize_site(site, sample_names, canonical, freq_thres, covg_thres,
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
283 stranded=False, debug=False):
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
284 """Take the raw data from the VCF line and transform it into the summary data
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
285 to be printed in the output format."""
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
286
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
287 site_summary = []
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
288 for sample_name in sample_names:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
289
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
290 sample = {'print':False}
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
291 variants = site['samples'].get(sample_name)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
292
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
293 sample['sample'] = sample_name
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
294 sample['chr'] = site['chr']
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
295 sample['pos'] = site['pos']
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
296
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
297 coverage = sum(variants.values())
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
298
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
299 # get stranded coverage
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
300 covg_plus = 0
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
301 covg_minus = 0
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
302 for variant in variants:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
303 if variant[0] == '+':
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
304 covg_plus += variants[variant]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
305 elif variant[0] == '-':
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
306 covg_minus += variants[variant]
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
307 # stranded coverage threshold
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
308 if covg_plus < covg_thres or covg_minus < covg_thres:
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
309 site_summary.append(sample)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
310 continue
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
311 else:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
312 sample['print'] = True
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
313
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
314 # get an ordered list of read counts for all variants (both strands)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
315 bases = get_read_counts(variants, '+-')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
316 ranked_bases = process_read_counts(bases, sort=True, debug=debug)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
317
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
318 # prepare stranded or unstranded lists of base counts
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
319 base_count_lists = []
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
320 if stranded:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
321 strands = ['+', '-']
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
322 base_count_lists.append(get_read_counts(variants, '+'))
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
323 base_count_lists.append(get_read_counts(variants, '-'))
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
324 else:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
325 strands = ['']
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
326 base_count_lists.append(ranked_bases)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
327
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
328 # record read counts into output dict
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
329 # If stranded, this will loop twice, once for each strand, and prepend '+'
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
330 # or '-' to the base name. If not stranded, it will loop once, and prepend
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
331 # nothing ('').
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
332 for (strand, base_count_list) in zip(strands, base_count_lists):
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
333 for base_count in base_count_list:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
334 sample[strand+base_count[0]] = base_count[1]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
335 # fill in any zeros
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
336 for base in canonical:
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
337 if strand+base not in sample:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
338 sample[strand+base] = 0
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
339
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
340 sample['alleles'] = count_alleles(variants, freq_thres, debug=debug)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
341
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
342 # If there's a tie for 2nd, randomly choose one to be 2nd
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
343 if len(ranked_bases) >= 3 and ranked_bases[1][1] == ranked_bases[2][1]:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
344 swap = random.choice([True,False])
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
345 if swap:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
346 tmp_base = ranked_bases[1]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
347 ranked_bases[1] = ranked_bases[2]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
348 ranked_bases[2] = tmp_base
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
349
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
350 if debug: print("ranked +-: "+str(ranked_bases))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
351
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
352 sample['coverage'] = coverage
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
353 try:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
354 sample['major'] = ranked_bases[0][0]
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
355 except IndexError:
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
356 sample['major'] = '.'
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
357 try:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
358 sample['minor'] = ranked_bases[1][0]
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
359 sample['freq'] = round(ranked_bases[1][1]/coverage, 5)
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
360 except IndexError:
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
361 sample['minor'] = '.'
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
362 sample['freq'] = 0.0
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
363
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
364 # Now that we've decided major and minor, we can calculate strand bias
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
365 bias = get_strand_bias(sample, variants)
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
366 if bias is None:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
367 sample['bias'] = '.'
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
368 else:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
369 sample['bias'] = round(bias, 5)
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
370
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
371 site_summary.append(sample)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
372
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
373 return site_summary
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
374
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
375
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
376 def get_read_counts(stranded_counts, strands='+-'):
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
377 """Do a simple sum of the read counts per variant, on the specified strands.
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
378 Arguments:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
379 stranded_counts: Dict of the stranded variants (keys) and their read counts
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
380 (values).
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
381 strands: Which strand(s) to count. Can be '+', '-', or '+-' for both (default).
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
382 Return value:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
383 summed_counts: A list of the alleles and their read counts. The elements are
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
384 tuples (variant, read count)."""
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
385
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
386 variants = stranded_counts.keys()
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
387
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
388 summed_counts = {}
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
389 for variant in variants:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
390 strand = variant[0]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
391 base = variant[1:]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
392 if strand in strands:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
393 summed_counts[base] = stranded_counts[variant] + summed_counts.get(base, 0)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
394
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
395 return list(summed_counts.items())
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
396
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
397
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
398 def process_read_counts(variant_counts, freq_thres=0, sort=False, debug=False):
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
399 """Process a list of read counts by frequency filtering and/or sorting.
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
400 Arguments:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
401 variant_counts: List of the non-stranded variants and their read counts. The
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
402 elements are tuples (variant, read count).
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
403 freq_thres: The frequency threshold each allele needs to pass to be included.
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
404 sort: Whether to sort the list in descending order of read counts.
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
405 Return value:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
406 variant_counts: A list of the processed alleles and their read counts. The
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
407 elements are tuples (variant, read count)."""
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
408
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
409 # get coverage for the specified strands
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
410 coverage = 0
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
411 for variant in variant_counts:
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
412 coverage += variant[1]
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
413
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
414 if coverage <= 0:
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
415 return []
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
416
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
417 # sort the list of alleles by read count
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
418 if sort:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
419 variant_counts.sort(reverse=True, key=lambda variant: variant[1])
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
420
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
421 if debug:
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
422 print('coverage: '+str(coverage)+', freq_thres: '+str(freq_thres))
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
423 for variant in variant_counts:
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
424 print((variant[0]+': '+str(variant[1])+'/'+str(float(coverage))+' = '+
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
425 str(variant[1]/coverage)))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
426
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
427 # remove bases below the frequency threshold
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
428 if freq_thres > 0:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
429 variant_counts = [variant for variant in variant_counts
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
430 if variant[1]/coverage >= freq_thres]
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
431
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
432 return variant_counts
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
433
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
434
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
435 def count_alleles(variant_counts, freq_thres, debug=False):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
436 """Determine how many alleles to report, based on filtering rules.
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
437 The current rule determines which bases pass the frequency threshold on each
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
438 strand individually, then compares the two sets of bases. If they are the same
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
439 (regardless of order), the allele count is the number of bases. Otherwise it
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
440 is zero."""
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
441 allele_count = 0
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
442
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
443 alleles_plus = get_read_counts(variant_counts, '+')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
444 alleles_plus = process_read_counts(alleles_plus, freq_thres=freq_thres,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
445 sort=False, debug=debug)
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
446 alleles_minus = get_read_counts(variant_counts, '-')
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
447 alleles_minus = process_read_counts(alleles_minus, freq_thres=freq_thres,
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
448 sort=False, debug=debug)
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
449
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
450 if debug:
9
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
451 print('+ '+str(alleles_plus))
6cc488e11544 "planemo upload for repository https://github.com/galaxyproject/dunovo commit 5a2e08bc1213b0437d0adcb45f7f431bd3c735f4"
nick
parents: 8
diff changeset
452 print('- '+str(alleles_minus))
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
453
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
454 # Check if each strand reports the same set of alleles.
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
455 # Sorting by base is to compare lists without regard to order (as sets).
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
456 alleles_plus_sorted = sorted([base[0] for base in alleles_plus if base[1]])
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
457 alleles_minus_sorted = sorted([base[0] for base in alleles_minus if base[1]])
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
458 if alleles_plus_sorted == alleles_minus_sorted:
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
459 allele_count = len(alleles_plus)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
460
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
461 return allele_count
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
462
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
463
6
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
464 def get_strand_bias(sample, variants):
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
465 """Based on method 1 (SB) of Guo et al., 2012
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
466 If there a denominator would be 0, there is no valid result and this will
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
467 return None. This occurs when there are no reads on one of the strands, or
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
468 when there are no minor allele reads."""
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
469 # using same variable names as in paper
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
470 a = variants.get('+'+sample['major'], 0) # forward major allele
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
471 b = variants.get('+'+sample['minor'], 0) # forward minor allele
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
472 c = variants.get('-'+sample['major'], 0) # reverse major allele
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
473 d = variants.get('-'+sample['minor'], 0) # reverse minor allele
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
474 # no minor allele reads
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
475 try:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
476 return abs(b/(a+b) - d/(c+d)) / ((b+d) / (a+b+c+d))
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
477 except ZeroDivisionError:
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
478 return None
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
479
df3b28364cd2 allele-counts.{py,xml}: Add strand bias, documentation updates.
nicksto <nmapsy@gmail.com>
parents: 5
diff changeset
480
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
481 def print_site(filehandle, site, columns):
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
482 """Print the output lines for one site (one per sample).
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
483 filehandle must be open."""
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
484 for sample in site:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
485 if sample['print']:
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
486 fields = [str(sample.get(column)) for column in columns]
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
487 filehandle.write('\t'.join(fields)+"\n")
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
488
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
489
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
490 def fail(message):
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
491 sys.stderr.write(message+'\n')
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
492 sys.exit(1)
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
493
5
31361191d2d2 Uploaded tarball.
nick
parents: 2
diff changeset
494
1
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
495 if __name__ == "__main__":
49bb46c3a1af Uploaded script
nick
parents:
diff changeset
496 main()