Mercurial > repos > nick > allele_counts
changeset 5:31361191d2d2
Uploaded tarball.
Version 1.1: Stranded output, slightly different handling of minor allele ties and 0 coverage sites, revised help text, added test datasets.
| author | nick |
|---|---|
| date | Thu, 12 Sep 2013 11:34:23 -0400 |
| parents | 898eb3daab43 |
| children | df3b28364cd2 |
| files | allele-counts.py allele-counts.xml tests/artificial-nofilt.csv.out tests/artificial-nofilt.vcf.in tests/artificial-samples.csv.out tests/artificial-samples.vcf.in tests/artificial.csv.out tests/artificial.vcf.in tests/compute-site.py tests/errors.vcf.in tests/real-mit-s.csv.out tests/real-mit-s.vcf.in tests/real-mit.csv.out tests/real-mit.vcf.in tests/real-nofilt.csv.out tests/real-nofilt.vcf.in tests/real.csv.out tests/real.vcf.in tests/run-tests.py |
| diffstat | 19 files changed, 837 insertions(+), 136 deletions(-) [+] |
line wrap: on
line diff
--- a/allele-counts.py Tue Jun 04 00:16:29 2013 -0400 +++ b/allele-counts.py Thu Sep 12 11:34:23 2013 -0400 @@ -1,25 +1,25 @@ #!/usr/bin/python # This parses the output of Dan's "Naive Variant Detector" (previously, # "BAM Coverage"). It was forked from the code of "bam-coverage.py". +# Or, to put it briefly, +# cat variants.vcf | grep -v '^#' | cut -f 10 | cut -d ':' -f 4 | tr ',=' '\t:' # # New in this version: # Made header line customizable # - separate from internal column labels, which are used as dict keys -# -# TODO: -# - test handling of -c 0 (and -f 0?) -# - should it technically handle data lines that start with a '#'? import os import sys +import random from optparse import OptionParser COLUMNS = ['sample', 'chr', 'pos', 'A', 'C', 'G', 'T', 'coverage', 'alleles', 'major', 'minor', 'freq'] #, 'bias'] COLUMN_LABELS = ['SAMPLE', 'CHR', 'POS', 'A', 'C', 'G', 'T', 'CVRG', 'ALLELES', 'MAJOR', 'MINOR', 'MINOR.FREQ.PERC.'] #, 'STRAND.BIAS'] CANONICAL_VARIANTS = ['A', 'C', 'G', 'T'] -USAGE = """Usage: cat variants.vcf | %prog [options] > alleles.csv - %prog [options] -i variants.vcf -o alleles.csv""" +USAGE = """Usage: %prog [options] -i variants.vcf -o alleles.csv + cat variants.vcf | %prog [options] > alleles.csv""" OPT_DEFAULTS = {'infile':'-', 'outfile':'-', 'freq_thres':1.0, 'covg_thres':100, - 'print_header':False, 'stdin':False} + 'print_header':False, 'stdin':False, 'stranded':False, 'no_filter':False, + 'debug_loc':'', 'seed':''} DESCRIPTION = """This will parse the VCF output of Dan's "Naive Variant Caller" (aka "BAM Coverage") Galaxy tool. For each position reported, it counts the number of reads of each base, determines the major allele, minor allele (second most frequent variant), and number of alleles above a threshold. So currently it only considers SNVs (ACGT), including in the coverage figure. By default it reads from stdin and prints to stdout.""" EPILOG = """Requirements: The input VCF must report the variants for each strand. @@ -40,27 +40,35 @@ help='Print output data to this file instead of stdout.') parser.add_option('-f', '--freq-thres', dest='freq_thres', type='float', default=defaults.get('freq_thres'), - help='Frequency threshold for counting alleles, given in percentage: -f 1 = 1% frequency. Default is %default%.') + help=('Frequency threshold for counting alleles, given in percentage: -f 1 ' + +'= 1% frequency. Default is %default%.')) parser.add_option('-c', '--covg-thres', dest='covg_thres', type='int', default=defaults.get('covg_thres'), - help='Coverage threshold. Each site must be supported by at least this many reads on each strand. Otherwise the site will not be printed in the output. The default is %default reads per strand.') + help=('Coverage threshold. Each site must be supported by at least this ' + +'many reads on each strand. Otherwise the site will not be printed in ' + +'the output. The default is %default reads per strand.')) + parser.add_option('-n', '--no-filter', dest='no_filter', action='store_const', + const=not(defaults.get('no_filter')), default=defaults.get('no_filter'), + help=('Operate without a frequency threshold or coverage threshold. ' + +'Equivalent to "-c 0 -f 0".')) parser.add_option('-H', '--header', dest='print_header', action='store_const', const=not(defaults.get('print_header')), default=defaults.get('print_header'), - help='Print header line. This is a #-commented line with the column labels. Off by default.') - parser.add_option('-d', '--debug', dest='debug', action='store_true', - default=False, - help='Turn on debug mode. You must also specify a single site to process in a final argument using UCSC coordinate format.') + help=('Print header line. This is a #-commented line with the column ' + +'labels. Off by default.')) + parser.add_option('-s', '--stranded', dest='stranded', action='store_const', + const=not(defaults.get('stranded')), default=defaults.get('stranded'), + help='Report variant counts by strand, in separate columns. Off by default.') + parser.add_option('-r', '--rand-seed', dest='seed', + default=defaults.get('seed'), help=('Seed for random number generator.')) + parser.add_option('-d', '--debug', dest='debug_loc', + default=defaults.get('debug_loc'), + help=('Turn on debug mode and specify a single site to process using UCSC ' + +'coordinate format. You can also append a sample ID after another ":" ' + +'to restrict it further.')) (options, args) = parser.parse_args() - # read in positional arguments - arguments = {} - if options.debug: - if len(args) >= 1: - arguments['print_loc'] = args[0] - args.remove(args[0]) - - return (options, arguments) + return (options, args) def main(): @@ -72,19 +80,26 @@ print_header = options.print_header freq_thres = options.freq_thres / 100.0 covg_thres = options.covg_thres - debug = options.debug + stranded = options.stranded + debug_loc = options.debug_loc + seed = options.seed + + if options.no_filter: + freq_thres = 0 + covg_thres = 0 - if debug: - print_loc = args.get('print_loc') - if print_loc: - if ':' in print_loc: - (print_chr, print_pos) = print_loc.split(':') - else: - print_pos = print_loc - else: - sys.stderr.write("Warning: No site coordinate found in arguments. " - +"Turning off debug mode.\n") - debug = False + if seed: + random.seed(seed) + + debug = False + print_sample = '' + if debug_loc: + debug = True + coords = debug_loc.split(':') + print_chr = coords[0] + print_pos = '' + if len(coords) > 1: print_pos = coords[1] + if len(coords) > 2: print_sample = coords[2] # set infile_handle to either stdin or the input file if infile == OPT_DEFAULTS.get('infile'): @@ -105,12 +120,19 @@ except IOError, e: fail('Error: The given output filename '+outfile+' could not be opened.') + # Take care of column names, print header if len(COLUMNS) != len(COLUMN_LABELS): - fail('Error: Internal column names do not match column labels.') + fail('Error: Internal column names list do not match column labels list.') + if stranded: + COLUMNS[3:7] = ['+A', '+C', '+G', '+T', '-A', '-C', '-G', '-T'] + COLUMN_LABELS[3:7] = ['+A', '+C', '+G', '+T', '-A', '-C', '-G', '-T'] if print_header: outfile_handle.write('#'+'\t'.join(COLUMN_LABELS)+"\n") - # main loop: process and print one line at a time + # main loop + # each iteration processes one VCF line and prints one or more output lines + # one VCF line = one site, one or more samples + # one output line = one site, one sample sample_names = [] for line in infile_handle: line = line.rstrip('\r\n') @@ -128,19 +150,24 @@ site_data = read_site(line, sample_names, CANONICAL_VARIANTS) if debug: - if site_data['pos'] != print_pos: + if print_pos != site_data['pos']: + continue + if print_chr != site_data['chr'] and print_chr != '': continue - try: - if site_data['chr'] != print_chr: - continue - except NameError, e: - pass # No chr specified. Just go ahead and print the line. + if print_sample != '': + for sample in site_data['samples'].keys(): + if sample.lower() != print_sample.lower(): + site_data['samples'].pop(sample, None) + if len(site_data['samples']) == 0: + sys.stderr.write("Error: Sample '"+print_sample+"' not found.\n") + sys.exit(1) + site_summary = summarize_site(site_data, sample_names, CANONICAL_VARIANTS, - freq_thres, covg_thres, debug=debug) + freq_thres, covg_thres, stranded, debug=debug) if debug and site_summary[0]['print']: - print line.split('\t')[9].split(':')[-1] + print line.split('\t')[9].split(':')[-1] print_site(outfile_handle, site_summary, COLUMNS) @@ -158,8 +185,27 @@ def read_site(line, sample_names, canonical): """Read in a line, parse the variants into a data structure, and return it. The line should be actual site data, not a header line, so check beforehand. - Notes: - - The line is assumed to have been chomped.""" + Only the variants in 'canonical' will be read; all others are ignored. + Note: the line is assumed to have been chomped. + The returned data is stored in a dict, with the following structure: + { + 'chr': 'chr1', + 'pos': '2617', + 'samples': { + 'THYROID': { + '+A': 32, + '-A': 45, + '-G': 1, + }, + 'BLOOD': { + '+A': 2, + '-C': 1, + '+G': 37, + '-G': 42, + }, + }, + } + """ site = {} fields = line.split('\t') @@ -212,7 +258,7 @@ def summarize_site(site, sample_names, canonical, freq_thres, covg_thres, - debug=False): + stranded=False, debug=False): """Take the raw data from the VCF line and transform it into the summary data to be printed in the output format.""" @@ -221,9 +267,6 @@ sample = {'print':False} variants = site['samples'].get(sample_name) - if not variants: - site_summary.append(sample) - continue sample['sample'] = sample_name sample['chr'] = site['chr'] @@ -240,31 +283,49 @@ elif variant[0] == '-': covg_minus += variants[variant] # stranded coverage threshold - if coverage <= 0 or covg_plus < covg_thres or covg_minus < covg_thres: + if covg_plus < covg_thres or covg_minus < covg_thres: site_summary.append(sample) continue else: sample['print'] = True - # get an ordered list of read counts for all variants (either strand) - ranked_bases = get_read_counts(variants, 0, strands='+-', debug=debug) + # get an ordered list of read counts for all variants (both strands) + bases = get_read_counts(variants, '+-') + ranked_bases = process_read_counts(bases, sort=True, debug=debug) + + # prepare stranded or unstranded lists of base counts + base_count_lists = [] + if stranded: + strands = ['+', '-'] + base_count_lists.append(get_read_counts(variants, '+')) + base_count_lists.append(get_read_counts(variants, '-')) + else: + strands = [''] + base_count_lists.append(ranked_bases) - # record read counts into dict for this sample - for base in ranked_bases: - sample[base[0]] = base[1] - # fill in any zeros - for variant in canonical: - if not sample.has_key(variant): - sample[variant] = 0 + # record read counts into output dict + # If stranded, this will loop twice, once for each strand, and prepend '+' + # or '-' to the base name. If not stranded, it will loop once, and prepend + # nothing (''). + for (strand, base_count_list) in zip(strands, base_count_lists): + for base_count in base_count_list: + sample[strand+base_count[0]] = base_count[1] + # fill in any zeros + for base in canonical: + if not sample.has_key(strand+base): + sample[strand+base] = 0 - sample['alleles'] = count_alleles(variants, freq_thres, debug=debug) + sample['alleles'] = count_alleles(variants, freq_thres, debug=debug) - # set minor allele to N if there's a tie for 2nd + # If there's a tie for 2nd, randomly choose one to be 2nd if len(ranked_bases) >= 3 and ranked_bases[1][1] == ranked_bases[2][1]: - ranked_bases[1] = ('N', 0) - sample['alleles'] = 1 if sample['alleles'] else 0 + swap = random.choice([True,False]) + if swap: + tmp_base = ranked_bases[1] + ranked_bases[1] = ranked_bases[2] + ranked_bases[2] = tmp_base - if debug: print ranked_bases + if debug: print "ranked +-: "+str(ranked_bases) sample['coverage'] = coverage try: @@ -283,67 +344,63 @@ return site_summary -def print_site(filehandle, site, columns): - """Print the output lines for one site (one per sample). - filehandle must be open.""" - for sample in site: - if sample['print']: - fields = [str(sample.get(column)) for column in columns] - filehandle.write('\t'.join(fields)+"\n") +def get_read_counts(stranded_counts, strands='+-'): + """Do a simple sum of the read counts per variant, on the specified strands. + Arguments: + stranded_counts: Dict of the stranded variants (keys) and their read counts + (values). + strands: Which strand(s) to count. Can be '+', '-', or '+-' for both (default). + Return value: + summed_counts: A list of the alleles and their read counts. The elements are + tuples (variant, read count).""" + + variants = stranded_counts.keys() + + summed_counts = {} + for variant in variants: + strand = variant[0] + base = variant[1:] + if strand in strands: + summed_counts[base] = stranded_counts[variant] + summed_counts.get(base, 0) + + return summed_counts.items() -def get_read_counts(variant_counts, freq_thres, strands='+-', debug=False): - """Count the number of reads for each base, and create a ranked list of - alleles passing the frequency threshold. +def process_read_counts(variant_counts, freq_thres=0, sort=False, debug=False): + """Process a list of read counts by frequency filtering and/or sorting. Arguments: - variant_counts: Dict of the stranded variants (keys) and their read counts (values). + variant_counts: List of the non-stranded variants and their read counts. The + elements are tuples (variant, read count). freq_thres: The frequency threshold each allele needs to pass to be included. - strands: Which strand(s) to count. Can be '+', '-', or '+-' for both (default). - variants: A list of the variants of interest. Other types of variants will not - be included in the returned list. If no list is given, all variants found in - the variant_counts will be used. + sort: Whether to sort the list in descending order of read counts. Return value: - ranked_bases: A list of the alleles and their read counts. The elements are - tuples (base, read count). The alleles are listed in descending order of - frequency, and only those passing the threshold are included.""" - - # Get list of all variants from variant_counts list - variants = [variant[1:] for variant in variant_counts] - # deduplicate via a dict - variant_dict = dict((variant, 1) for variant in variants) - variants = variant_dict.keys() - - ranked_bases = [] - for variant in variants: - reads = 0 - for strand in strands: - reads += variant_counts.get(strand+variant, 0) - ranked_bases.append((variant, reads)) + variant_counts: A list of the processed alleles and their read counts. The + elements are tuples (variant, read count).""" # get coverage for the specified strands coverage = 0 for variant in variant_counts: - if variant[0] in strands: - coverage += variant_counts.get(variant, 0) - # if debug: print "strands: "+strands+', covg: '+str(coverage) + coverage += variant[1] - if coverage < 1: + if coverage <= 0: return [] # sort the list of alleles by read count - ranked_bases.sort(reverse=True, key=lambda base: base[1]) + if sort: + variant_counts.sort(reverse=True, key=lambda variant: variant[1]) if debug: - print strands+' coverage: '+str(coverage)+', freq_thres: '+str(freq_thres) - for base in ranked_bases: - print (base[0]+': '+str(base[1])+'/'+str(float(coverage))+' = '+ - str(base[1]/float(coverage))) + print 'coverage: '+str(coverage)+', freq_thres: '+str(freq_thres) + for variant in variant_counts: + print (variant[0]+': '+str(variant[1])+'/'+str(float(coverage))+' = '+ + str(variant[1]/float(coverage))) # remove bases below the frequency threshold - ranked_bases = [base for base in ranked_bases - if base[1]/float(coverage) >= freq_thres] + if freq_thres > 0: + variant_counts = [variant for variant in variant_counts + if variant[1]/float(coverage) >= freq_thres] - return ranked_bases + return variant_counts def count_alleles(variant_counts, freq_thres, debug=False): @@ -354,16 +411,19 @@ is zero.""" allele_count = 0 - alleles_plus = get_read_counts(variant_counts, freq_thres, debug=debug, - strands='+') - alleles_minus = get_read_counts(variant_counts, freq_thres, debug=debug, - strands='-') + alleles_plus = get_read_counts(variant_counts, '+') + alleles_plus = process_read_counts(alleles_plus, freq_thres=freq_thres, + sort=False, debug=debug) + alleles_minus = get_read_counts(variant_counts, '-') + alleles_minus = process_read_counts(alleles_minus, freq_thres=freq_thres, + sort=False, debug=debug) if debug: print '+ '+str(alleles_plus) print '- '+str(alleles_minus) - # check if each strand reports the same set of alleles + # Check if each strand reports the same set of alleles. + # Sorting by base is to compare lists without regard to order (as sets). alleles_plus_sorted = sorted([base[0] for base in alleles_plus if base[1]]) alleles_minus_sorted = sorted([base[0] for base in alleles_minus if base[1]]) if alleles_plus_sorted == alleles_minus_sorted: @@ -372,9 +432,19 @@ return allele_count +def print_site(filehandle, site, columns): + """Print the output lines for one site (one per sample). + filehandle must be open.""" + for sample in site: + if sample['print']: + fields = [str(sample.get(column)) for column in columns] + filehandle.write('\t'.join(fields)+"\n") + + def fail(message): sys.stderr.write(message+'\n') sys.exit(1) + if __name__ == "__main__": main() \ No newline at end of file
--- a/allele-counts.xml Tue Jun 04 00:16:29 2013 -0400 +++ b/allele-counts.xml Thu Sep 12 11:34:23 2013 -0400 @@ -1,10 +1,12 @@ -<tool id="allele_counts_1" version="1.0" name="Count alleles"> - <description>and minor allele frequencies</description> - <command interpreter="python">allele-counts.py -i $input -o $output -f $freq -c $covg $header</command> +<tool id="allele_counts_1" version="1.1" name="Variant Annotator"> + <description> process variant counts</description> + <command interpreter="python">allele-counts.py -i $input -o $output -f $freq -c $covg $header $stranded $nofilt</command> <inputs> <param name="input" type="data" format="vcf" label="Input variants from Naive Variants Detector"/> <param name="freq" type="float" value="1.0" min="0" max="100" label="Minor allele frequency threshold (in percent)"/> - <param name="covg" type="integer" value="10" min="0" label="Coverage threshold (per strand)"/> + <param name="covg" type="integer" value="10" min="0" label="Coverage threshold (in reads per strand)"/> + <param name="nofilt" type="boolean" truevalue="-n" falsevalue="" checked="False" label="Do not filter sites or alleles" /> + <param name="stranded" type="boolean" truevalue="-s" falsevalue="" checked="False" label="Output stranded base counts" /> <param name="header" type="boolean" truevalue="-H" falsevalue="" checked="True" label="Write header line" /> </inputs> <outputs> @@ -21,40 +23,51 @@ **What it does** -This tool parses variant counts from a special VCF file (normally the output of the **Naive Variant Detector** tool). It counts simple (ACGT) variants, calculates numbers of alleles, and calculates minor allele frequency. It applies filters based on coverage, strand bias, and minor allele frequency cutoffs. +This tool parses variant counts from a special VCF file. It counts simple variants, calculates numbers of alleles, and calculates minor allele frequency. It can apply filters based on coverage, strand bias, and minor allele frequency cutoffs. ----- +.. class:: infomark + +**Input Format** + .. class:: warningmark -**Note** +**Note:** variants that are not A/C/G/T SNVs will be ignored! -The VCF must have a certain genotype field in the sample columns, giving the read count of each type of variant. Also, the variant data **must be stranded**. The **Naive Variant Detector** tool produces this type of VCF. +The input VCF should be like the output of the **Naive Variant Detector** tool (using the stranded option). The sample column(s) must give the read count for each variant **on each strand**. Below is an example of a valid sample column entry (the important part is after the last colon):: + + 0/0:1:0.02:+T=27,+G=1,-T=22, ----- .. class:: infomark -**Output columns** +**Output** -Each row represents one site in one sample. 12 fields give information about that site:: +Each row represents one site in one sample. For unstranded output, 12 fields give information about that site:: - 1. SAMPLE - Sample names (from VCF sample column labels) + 1. SAMPLE - Sample name (from VCF sample column labels) 2. CHR - Chromosome of the site 3. POS - Chromosomal coordinate of the site 4. A - Number of reads supporting an 'A' - 5. C - ditto, for 'C' - 6. G - ditto, for 'G' - 7. T - ditto, for 'T' + 5. C - 'C' reads + 6. G - 'G' reads + 7. T - 'T' reads 8. CVRG - Total (number of reads supporting one of the four bases above) 9. ALLELES - Number of qualifying alleles - 10. MAJOR - Major allele base - 11. MINOR - Minor allele base (2nd most prevalent variant) + 10. MAJOR - Major allele + 11. MINOR - Minor allele (2nd most prevalent variant) 12. MINOR.FREQ.PERC. - Frequency of minor allele +For stranded output, instead of using 4 columns to report read counts per base, 8 are used to report the stranded counts per base:: + + 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 + SAMPLE CHR POS +A +C +G +T -A -C -G -T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. + **Example** -This is the header line, followed by some example data lines. Note that some samples and/or sites will not be included in the output, if they fall below the coverage threshold:: +Below is a header line, followed by some example data lines. Since the input contained three samples, the data for each site is reported on three consecutive lines. However, if a sample fell below the coverage threshold at that site, the line will be omitted:: #SAMPLE CHR POS A C G T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. BLOOD_1 chr20 99 0 101 1 2 104 1 C T 0.01923 @@ -69,19 +82,17 @@ **Site printing and allele tallying requirements** -Each line is printed only when the site is covered by the threshold number of reads **on each strand**. If coverage of either strand is below the threshold, the line (sample + site combination) is omitted. +Coverage threshold: -**N.B.**: This means the total coverage for each printed site will be at least twice the number you give in the "coverage threshold" option. +If a coverage threshold is used, the number of reads **on each strand** must be at or above the threshold. If either strand is below the threshold, the line will be omitted. **N.B.** this means the total coverage for each printed site will be at least twice the number you give in the "coverage threshold" option. Also, since only simple variants are counted, a site with 100 reads, all supporting a deletion variant, would not be printed. -Also, reads supporting a variant outside the canonical 4 nucleotides will not count towards the coverage requirement. For instance, a site/sample line with 100x coverage, all of which support a deletion variant, will not be printed. +Frequency threshold: -Alleles are only counted (in column 9) if they meet or exceed the minor allele frequency threshold. So a site/sample line with types of variants, 96% A, 3.3% C, and 0.7% G, will count as 2 alleles (at 1% threshold). - -Strand bias: the alleles passing the threshold on each strand have to match (though not in order). Otherwise, the allele count will be 0. So a site/sample line whose + strand shows 70% A, 27% C, and 3% G, and - strand shows 70% A and 30% C will have an allele count of 0. The minor allele and minor allele frequency, though, will always be reported\*. +If a frequency threshold is used, alleles are only counted (in the ALLELES column) if they meet or exceed this minor allele frequency threshold. -But in this version, there is no requirement that the strands show similar allele frequencies, as long as they both pass the threshold. +Strand bias: -\*One specific case will actually affect the reported minor allele identity and frequency. If there is a tie for the minor allele (between the 2nd and 3rd most common alleles), the minor allele will be reporated as 'N', and the frequency as 0.0. +The alleles passing the threshold on each strand must match (though not in order), or the allele count will be 0. So a site with A, C, G on the plus strand and A, G on the minus strand will get an allele count of zero, though the (strand-independent) major allele, minor allele, and minor allele frequency will still be reported. If there is a tie for the minor allele, one will be randomly chosen. </help>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial-nofilt.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,27 @@ +#SAMPLE CHR POS A C G T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. +THYROID chr1 0 30 0 0 0 30 1 A . 0.0 +THYROID chr1 10 30 0 2 0 32 2 A G 0.0625 +THYROID chr1 20 31 0 1 0 32 0 A G 0.03125 +THYROID chr1 30 21 0 4 0 25 2 A G 0.16 +THYROID chr1 40 22 0 3 0 25 0 A G 0.12 +THYROID chr1 50 3 0 0 0 3 1 A . 0.0 +THYROID chr1 60 2 0 2 0 4 2 A G 0.5 +THYROID chr1 70 1 0 3 0 4 0 G A 0.25 +THYROID chr1 80 104 0 3 0 107 0 A G 0.02804 +THYROID chr1 90 100 2 11 0 113 3 A G 0.09735 +THYROID chr1 100 100 1 11 0 112 0 A G 0.09821 +THYROID chr1 120 0 0 0 0 0 0 . . 0.0 +THYROID chr1 130 0 0 2 0 2 1 G . 0.0 +THYROID chr1 140 0 0 1 0 1 0 G . 0.0 +THYROID chr1 150 0 0 4 0 4 1 G . 0.0 +THYROID chr1 160 0 0 3 0 3 0 G . 0.0 +THYROID chr1 260 106 0 14 0 120 2 A G 0.11667 +THYROID chr1 300 2 0 2 76 80 3 T G 0.025 +THYROID chr1 310 12 0 12 76 100 3 T G 0.12 +THYROID chr1 320 12 0 12 56 80 3 T A 0.15 +THYROID chr1 330 7 0 7 66 80 3 T G 0.0875 +THYROID chr1 340 1 0 1 98 100 0 T G 0.01 +THYROID chr1 350 11 0 11 78 100 0 T A 0.11 +THYROID chr1 400 32 0 8 0 40 2 A G 0.2 +THYROID chr1 410 1 0 2 97 100 0 T G 0.02 +THYROID chr1 420 104 0 0 0 104 1 A . 0.0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial-nofilt.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,50 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-H -r 1 -f 0 -c 0" +##comment="This is a test set of made-up sites, each created in order to test certain functionality. It's meant to be run with -f 10 -c 10" +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT THYROID +# General note: the only data made consistent is the CHROM, POS, REF, ALT, and the variant data (after the ':'). The other stuff isn't supposed to be consistent. +# Simplest case, but POS 0 and no minor allele +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, +# Simple, normal cases of A/G variants: above/below threshold x strand bias/no strand bias (2 x 2 = 4 cases) +chr1 10 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,-A=15,-G=1, +chr1 20 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,-A=16, +chr1 30 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,+G=2,-A=11,-G=2, +chr1 40 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+G=3,-A=11, +# Really low coverage +chr1 50 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=2,-A=1, +chr1 60 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+G=1,-A=1,-G=1 +chr1 70 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+G=1,-G=2, +chr1 80 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=102,+G=3,-A=2, +# Very low frequency tertiary allele +chr1 90 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=50,+G=5,+C=1,-A=50,-G=6,-C=1, +chr1 100 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=50,+G=5,-A=50,-G=6,-C=1, +# Only non-canonical alleles +chr1 120 . d1 . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,-d1=15, +# Same 4 cases, but with non-canonical major allele (d1) +chr1 130 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,+G=1,-d1=15,-G=1, +chr1 140 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,+G=1,-d1=16, +chr1 150 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=10,+G=2,-d1=11,-G=2, +chr1 160 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=11,+G=3,-d1=11, +# Test case where minor allele is above threshold on only one strand because of different coverage. Also, a long decimal minor allele frequency. +chr1 260 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=13,+G=7,-A=93,-G=7, +# Test case where minor alleles have equal frequency: Above/below threshold, +/- strand bias +chr1 300 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+G=1,+T=38,-A=1,-G=1,-T=38, +chr1 310 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=6,-G=6,-T=38, +chr1 320 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=6,-G=6,-T=18, +chr1 330 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=1,-G=1,-T=28, +chr1 340 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+T=80,-G=1,-T=18, +chr1 350 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+T=60,-G=11,-T=18, +# Case where + and - variants are interleaved with each other. Also, a long decimal result for the minor allele frequency. +chr1 400 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=16,-A=16,+G=4,-G=4, +# Test complex data in the ALT, INFO, and sample (before the ':') columns +chr1 410 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=16,-G=2, +chr1 420 . A . . . AC=;AF= GT:AC:AF:NC 0/0:::+A=82,-A=22, \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial-samples.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,13 @@ +BRAIN chr1 0 30 0 0 0 30 1 A . 0.0 +ARTERY chr1 0 0 0 30 0 30 1 G . 0.0 +THYROID chr1 0 0 30 0 0 30 1 C . 0.0 +BRAIN chr1 10 30 0 0 0 30 1 A . 0.0 +ARTERY chr1 10 30 0 2 0 32 1 A G 0.0625 +THYROID chr1 10 31 0 1 0 32 1 A G 0.03125 +BRAIN chr1 20 30 0 2 0 32 1 A G 0.0625 +ARTERY chr1 20 34 0 6 0 40 2 A G 0.15 +THYROID chr1 20 30 0 2 0 32 0 A G 0.0625 +BRAIN chr1 30 30 0 0 0 30 1 A . 0.0 +BRAIN chr1 40 0 0 0 30 30 1 T . 0.0 +ARTERY chr1 40 1 0 2 97 100 0 T G 0.02 +THYROID chr1 40 0 69 0 31 100 0 C T 0.31
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial-samples.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,23 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-f 10 -c 10" +##comment="This is a test set of made-up sites, in this case to test the handling of multiple read groups. It's meant to be run with -f 10 -c 10" +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT BRAIN ARTERY THYROID +# naive case +chr1 0 . A G,C . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, 0/0:1:1:+G=15,-G=15, 0/0:1:1:+C=15,-C=15, +# cases all with one allele, but different minor allele numbers below threshold +chr1 10 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, 0/0:1:1:+A=15,+G=1,-A=15,-G=1, 0/0:1:1:+A=15,+G=1,-A=16, +# cases with different numbers of alleles +chr1 20 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,-A=15,-G=1, 0/0:1:1:+A=17,+G=3,-A=17,-G=3, 0/0:1:1:+A=15,+G=2,-A=15, +# cases where some shouldn't be printed +chr1 30 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, 0/0:1:1:+A=9,-A=8, 0/0:1:1:+A=15,-A=9, +# cases with some complex pre-':' data +chr1 40 . T A,G,C . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1:1:+T=15,-T=15, 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=16,-G=2, 0/0:3,1:0.0394736842105,0.0131578947368:+C=52,+T=30,-C=17,-T=1, \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,35 @@ +THYROID chr1 0 30 0 0 0 30 1 A . 0.0 +THYROID chr1 10 30 0 2 0 32 1 A G 0.0625 +THYROID chr1 20 31 0 1 0 32 1 A G 0.03125 +THYROID chr1 30 21 0 4 0 25 2 A G 0.16 +THYROID chr1 40 22 0 3 0 25 0 A G 0.12 +THYROID chr1 50 30 0 0 0 30 1 A . 0.0 +THYROID chr1 60 31 0 0 0 31 1 A . 0.0 +THYROID chr1 70 21 0 0 0 21 1 A . 0.0 +THYROID chr1 80 22 0 0 0 22 1 A . 0.0 +THYROID chr1 82 30 0 2 0 32 1 A G 0.0625 +THYROID chr1 84 31 0 1 0 32 1 A G 0.03125 +THYROID chr1 86 21 0 4 0 25 2 A G 0.16 +THYROID chr1 88 22 0 3 0 25 0 A G 0.12 +THYROID chr1 90 30 0 0 0 30 1 A . 0.0 +THYROID chr1 100 31 0 0 0 31 1 A . 0.0 +THYROID chr1 110 21 0 0 0 21 1 A . 0.0 +THYROID chr1 120 22 0 0 0 22 1 A . 0.0 +THYROID chr1 210 20 0 0 0 20 1 A . 0.0 +THYROID chr1 220 22 0 0 0 22 1 A . 0.0 +THYROID chr1 230 182 0 18 0 200 1 A G 0.09 +THYROID chr1 240 180 0 20 0 200 2 A G 0.1 +THYROID chr1 250 178 0 22 0 200 2 A G 0.11 +THYROID chr1 260 106 0 14 0 120 0 A G 0.11667 +THYROID chr1 300 2 0 2 76 80 1 T G 0.025 +THYROID chr1 310 12 0 12 76 100 3 T G 0.12 +THYROID chr1 320 12 0 12 56 80 3 T A 0.15 +THYROID chr1 330 7 0 7 66 80 0 T G 0.0875 +THYROID chr1 340 1 0 1 98 100 1 T G 0.01 +THYROID chr1 350 11 0 11 78 100 0 T A 0.11 +THYROID chr1 400 32 0 8 0 40 2 A G 0.2 +THYROID chr1 410 1 0 2 97 100 0 T G 0.02 +THYROID chr1 420 104 0 0 0 104 1 A . 0.0 +THYROID chr1 430 30 0 0 0 30 1 A . 0.0 +THYROID chr1 440 30 0 0 0 30 1 A . 0.0 +THYROID 27 1234567890 29 0 0 0 29 1 A . 0.0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/artificial.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,72 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-r 1 -f 10 -c 10" +##comment="This is a test set of made-up sites, each created in order to test certain functionality. It's meant to be run with -f 10 -c 10" +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT THYROID +# General note: the only data made consistent is the CHROM, POS, REF, ALT, and the variant data (after the ':'). The other stuff isn't supposed to be consistent. +# Simplest case, but POS 0 and no minor allele +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, +# Simple, normal cases of A/G variants: above/below threshold x strand bias/no strand bias (2 x 2 = 4 cases) +chr1 10 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,-A=15,-G=1, +chr1 20 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,-A=16, +chr1 30 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,+G=2,-A=11,-G=2, +chr1 40 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+G=3,-A=11, +# Same 4 cases, but with minor allele = N +chr1 50 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+N=1,-A=15,-N=1, +chr1 60 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+N=1,-A=16, +chr1 70 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,+N=2,-A=11,-N=2, +chr1 80 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+N=3,-A=11, +# Same 4 cases, but with an additional noncanonical minor allele d1 +chr1 82 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,+d1=1,-A=15,-G=1,-d1=1, +chr1 84 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,+d1=1,-A=16, +chr1 86 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,+G=2,+d1=1,-A=11,-G=2,-d1=1, +chr1 88 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+G=3,+d1=1,-A=11, +# Same 4 cases, but with minor allele = d1 (non-canonical) +chr1 90 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+d1=1,-A=15,-d1=1, +chr1 100 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+d1=1,-A=16, +chr1 110 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,+d1=2,-A=11,-d1=2, +chr1 120 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+d1=3,-A=11, +# Same 4 cases, but with MAJOR allele = d1 +chr1 130 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,+G=1,-d1=15,-G=1, +chr1 140 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,+G=1,-d1=16, +chr1 150 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=10,+G=2,-d1=11,-G=2, +chr1 160 . d1 G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=11,+G=3,-d1=11, +# Test edge cases where freq == freq_thres and/or covg == covg_thres +chr1 200 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=9,-A=9, +chr1 210 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=10,-A=10, +chr1 220 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,-A=11, +chr1 230 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=91,+G=9,-A=91,-G=9, +chr1 240 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=90,+G=10,-A=90,-G=10, +chr1 250 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=89,+G=11,-A=89,-G=11, +# Test case where minor allele is above threshold on only one strand because of different coverage. Also, a long decimal minor allele frequency. +chr1 260 . A G . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=13,+G=7,-A=93,-G=7, +# Test case where minor alleles have equal frequency: Above/below threshold, +/- strand bias +chr1 300 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+G=1,+T=38,-A=1,-G=1,-T=38, +chr1 310 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=6,-G=6,-T=38, +chr1 320 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=6,-G=6,-T=18, +chr1 330 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=6,+G=6,+T=38,-A=1,-G=1,-T=28, +chr1 340 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=1,+T=80,-G=1,-T=18, +chr1 350 . T G,A . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=11,+T=60,-G=11,-T=18, +# Case where + and - variants are interleaved with each other. Also, a long decimal result for the minor allele frequency. +chr1 400 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=16,-A=16,+G=4,-G=4, +# Test complex data in the ALT, INFO, and sample (before the ':') columns +chr1 410 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=16,-G=2, +chr1 420 . A . . . AC=;AF= GT:AC:AF:NC 0/0:::+A=82,-A=22, +# Test some other types of noncanonical variants (tie for 2nd and not) +chr1 430 . A N,GAA,d2 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+N=1,+d2=1,-A=15,-GAA=2 +chr1 440 . A N,GAA,d2 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+N=1,+d2=1,-A=15,-GAA=1 +# No canonical variants present +chr1 450 . A d1 . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,-d1=20, +chr1 460 . A d1,N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+d1=15,+N=2,-d1=20,-N=2, +# Catch some divide by zero errors +chr1 470 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=0, +# Test an unusual CHROM value and a long POS value +27 1234567890 . A N . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+N=1,-A=14, \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/compute-site.py Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,161 @@ +#!/usr/bin/python +import os +import sys +from optparse import OptionParser + +OPT_DEFAULTS = {'freq_thres':0, 'covg_thres':0} +BASES = ['A', 'C', 'G', 'T'] +USAGE = ("Usage: %prog (options) 'variant_str' ('variant_str2' (etc))\n" + +" cat variants.txt > %prog (options)") + +def main(): + + parser = OptionParser(usage=USAGE) + parser.add_option('-f', '--freq-thres', dest='freq_thres', type='float', + default=OPT_DEFAULTS.get('freq_thres'), + help=('Frequency threshold for counting alleles, given in percentage: -f 1 ' + +'= 1% frequency. Default is %default%.')) + parser.add_option('-c', '--covg-thres', dest='covg_thres', type='int', + default=OPT_DEFAULTS.get('covg_thres'), + help=('Coverage threshold. Each site must be supported by at least this ' + +'many reads on each strand. Otherwise the site will not be printed in ' + +'the output. The default is %default reads per strand.')) + (options, args) = parser.parse_args() + freq_thres = options.freq_thres + covg_thres = options.covg_thres + + if len(sys.argv) > 1 and '-h' in sys.argv[1][0:3]: + script_name = os.path.basename(sys.argv[0]) + print """USAGE: + $ """+script_name+""" [sample column text] + $ """+script_name+""" '+A=10,+G=2,-A=20,-G=41,' + $ """+script_name+""" '0/1:10,1:0.25,0.025:-A=29,-AG=1,-G=10,' + $ """+script_name+""" '+A=10,+G=2,-A=20,-G=41,' '0/1:10,1:0.25,0.025:-A=29,-AG=1,-G=10,' +Or invoke with no arguments to use interactively. It will read from stdin, so +just paste one sample per line.""" + sys.exit(0) + + if len(args) > 0: + stdin = False + samples = args + else: + stdin = True + samples = sys.stdin + print "Reading from standard input.." + + for sample in samples: + print '' + sample = sample.split(':')[-1] + if not sample: + continue + print sample + counts_dict = parse_counts(sample) + compute_stats(counts_dict, freq_thres, covg_thres) + + +def parse_counts(sample_str): + counts_dict = {} + counts = sample_str.split(',') + for count in counts: + if '=' not in count: + continue + (var, reads) = count.split('=') + if var[1:] in BASES: + counts_dict[var] = int(reads) + return counts_dict + + +def compute_stats(counts_dict, freq_thres, covg_thres): + # totals for A, C, G, T + counts_unstranded = {} + for base in BASES: + counts_unstranded[base] = 0 + for strand in '+-': + counts_unstranded[base] += counts_dict.get(strand+base, 0) + print '+- '+str(counts_unstranded) + + # get counts for each strand + plus_counts = get_stranded_counts(counts_dict, '+') + minus_counts = get_stranded_counts(counts_dict, '-') + counts_lists = {'+':plus_counts, '-':minus_counts} + print ' + '+str(plus_counts) + print ' - '+str(minus_counts) + + # stranded coverage threshold + coverages = {} + failing_strands = {} + for strand in '+-': + coverages[strand] = 0 + for count in counts_lists[strand].values(): + coverages[strand] += count + if coverages[strand] < covg_thres: + failing_strands[strand] = coverages[strand] + sys.stdout.write(strand+'coverage: '+str(coverages[strand])+"\t") + coverages['+-'] = coverages['+'] + coverages['-'] + sys.stdout.write("+-coverage: "+str(coverages['+-'])+"\n") + + if failing_strands: + for strand in failing_strands: + print ('coverage on '+strand+' strand too low (' + +str(failing_strands[strand])+' < '+str(covg_thres)+")") + return + + # apply frequency threshold + for strand in counts_lists: + strand_counts = counts_lists[strand] + for variant in strand_counts.keys(): + # print (variant+" freq: "+str(strand_counts[variant])+"/" + # +str(coverages[strand])+" = " + # +str(strand_counts[variant]/float(coverages[strand]))) + if strand_counts[variant]/float(coverages[strand]) < freq_thres: + strand_counts.pop(variant) + plus_variants = sorted(plus_counts.keys()) + minus_variants = sorted(minus_counts.keys()) + if plus_variants == minus_variants: + strand_bias = False + print "no strand bias: +"+str(plus_variants)+" == -"+str(minus_variants) + sys.stdout.write("alleles: "+str(len(plus_variants))+"\t") + else: + strand_bias = True + print " strand bias: +"+str(plus_variants)+" != -"+str(minus_variants) + sys.stdout.write("alleles: 0\t") + + variants_sorted = sort_variants(counts_unstranded) + if len(variants_sorted) >= 1: + sys.stdout.write("major: "+variants_sorted[0]+"\t") + minor = '.' + if len(variants_sorted) == 2: + minor = variants_sorted[1] + elif len(variants_sorted) > 2: + if (counts_unstranded.get(variants_sorted[1]) == + counts_unstranded.get(variants_sorted[2])): + minor = 'N' + else: + minor = variants_sorted[1] + + sys.stdout.write("minor: "+minor+"\tfreq: ") + print round(float(counts_unstranded.get(minor, 0))/coverages['+-'], 5) + + +def get_stranded_counts(unstranded_counts, strand): + stranded_counts = {} + for variant in unstranded_counts: + if variant[0] == strand: + stranded_counts[variant[1:]] = unstranded_counts[variant] + return stranded_counts + +def sort_variants(variant_counts): + """Sort the list of variants based on their counts. Returns a list of just + the variants, no counts.""" + variants = variant_counts.keys() + var_del = [] + for variant in variants: + if variant_counts.get(variant) == 0: + var_del.append(variant) + for variant in var_del: + variants.remove(variant) + variants.sort(reverse=True, key=lambda variant: variant_counts.get(variant,0)) + return variants + +if __name__ == "__main__": + main()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/errors.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,25 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-f 10 -c 10" +##comment="This is a test set of made-up sites, each created in order to test a certain type of error handling. It's meant to be run with -f 10 -c 10" +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT THYROID +# A correct control +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, +# List some bases more than once +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15,+A=12, +# Comma omissions +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15 +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15-A=15, +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1,+C=1,-A=15,-G=1, +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,+G=1+C=1,-A=15,-G=1, +# Wrong number of sample columns +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC 0/0:1:1:+A=15,-A=15, 0/0:1:1:+A=15,-A=15, +chr1 0 . A . . . AC=;AF= GT:AC:AF:NC \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-mit-s.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,12 @@ +#SAMPLE CHR POS +A +C +G +T -A -C -G -T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. +S1 chrM 2000 1 9095 1 0 7 5808 0 1 14913 1 C A 0.00054 +S3 chrM 2000 0 7933 0 4 10 5242 1 2 13192 1 C A 0.00076 +S1 chrM 3000 17399 7 22 8 10567 35 22 4 28064 0 A G 0.00157 +S2 chrM 3000 12535 3 24 2 7937 13 12 2 20528 1 A G 0.00175 +S3 chrM 3000 18981 7 29 6 11286 33 17 4 30363 0 A G 0.00152 +S4 chrM 3000 9254 1 15 2 6240 16 14 1 15543 0 A G 0.00187 +S1 chrM 4000 6134 2 1 3 6124 1 1 1 12267 1 A T 0.00033 +S1 chrM 7000 0 17 1 6216 4 9 2 7529 13778 0 T C 0.00189 +S2 chrM 7000 0 7 2 5104 0 9 4 6288 11414 1 T C 0.0014 +S3 chrM 7000 0 9 0 6446 4 4 10 7506 13979 1 T C 0.00093 +S3 chrM 8000 3 1 5023 1 1 0 5043 2 10074 1 G A 0.0004
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-mit-s.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,28 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-f 0.2 -c 5000 -H -s" +##fileDate=20130521 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_2536.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S1 S2 S3 S4 +chrM 1000 . ACTC TCTC,CCTC,GCTC,ATC,A . . AC=188032,682,8,8,1;AF=0.995310134556,0.00361003186568,4.23464148467e-05,4.23464148467e-05,5.29330185583e-06 GT:AC:AF:NC 1/1:12546,51,0,0,0:0.994688020297,0.00404344723698,0,0,0:+A=2,+C=20,+T=8020,-A=14,-C=31,-T=4526, 1/1:9970,38,1,1,0:0.995208624476,0.00379317228988,9.98203234178e-05,9.98203234178e-05,0:+A=1,+d1=1,+C=21,+T=6546,+G=1,-A=7,-C=17,-T=3424, 1/1:11372,35,0,2,0:0.995448179272,0.0030637254902,0,0.000175070028011,0:+d1=2,+C=10,+T=7199,-A=15,-C=25,-T=4173, 1/1:9115,39,0,0,0:0.995304651671,0.00425857174055,0,0,0:+C=13,+T=5817,-A=4,-C=26,-T=3298, +chrM 2000 . TAC CAC,AAC,GAC,T,TC,TGAC . . AC=209371,114,8,2,2,1;AF=0.999083807733,0.000543989158392,3.81746777819e-05,9.54366944547e-06,9.54366944547e-06,4.77183472273e-06 GT:AC:AF:NC 1/1:14903,8,1,0,0,0:0.999329444109,0.000536444712667,6.70555890834e-05,0,0,0:+A=1,+C=9095,+G=1,-A=7,-C=5808,-T=1, 1/1:11838,3,0,0,0,0:0.999324666554,0.000253250042208,0,0,0,0:+C=7242,+T=4,-A=3,-C=4596,-T=1, 1/1:13175,10,1,0,0,0:0.998711340206,0.000758035172832,7.58035172832e-05,0,0,0:+C=7933,+T=4,-A=10,-C=5242,-T=2,-G=1, 1/1:10921,3,2,0,0,1:0.999176578225,0.000274473924977,0.000182982616651,0,0,9.14913083257e-05:+A=1,+TG=1,+C=6732,+T=1,-A=2,-C=4189,-T=2,-G=2, +chrM 3000 . TC AC,GC,CC,T,TGC . . AC=386813,649,441,35,33;AF=0.996740346013,0.00167234421946,0.00113636949273,9.01880549786e-05,8.5034451837e-05 GT:AC:AF:NC 1/1:27966,44,42,5,6:0.996117542297,0.00156723063224,0.00149599287622,0.000178094390027,0.000213713268032:+A=17399,+d1=1,+C=7,+T=8,+G=22,-A=10567,-d1=4,-C=35,-G=22,-T=4,-TG=6, 1/1:20472,36,16,0,0:0.997272018706,0.00175370226033,0.000779423226812,0,0:+A=12535,+C=3,+T=2,+G=24,-A=7937,-C=13,-T=2,-G=12, 1/1:30267,46,40,5,3:0.996575680748,0.00151460274604,0.00131704586612,0.000164630733265,9.87784399592e-05:+A=18981,+d1=2,+C=7,+T=6,+G=29,-A=11286,-d1=3,-C=33,-G=17,-T=4,-TG=3, 1/1:15494,29,17,2,3:0.996526884487,0.00186519166452,0.00109338821713,0.000128633907898,0.000192950861847:+A=9254,+d1=2,+C=1,+T=2,+G=15,-A=6240,-TG=3,-C=16,-T=1,-G=14, +chrM 4000 . TA AA,GA,CA,T,TCA,TTA . . AC=138825,77,37,9,8,3;AF=0.998511134127,0.000553829334254,0.000266125783992,6.47332988089e-05,5.75407100524e-05,2.15777662696e-05 GT:AC:AF:NC 1/1:12258,2,3,1,1,2:0.998940591639,0.00016298590172,0.000244478852579,8.14929508598e-05,8.14929508598e-05,0.00016298590172:+A=6134,+d1=1,+C=2,+G=1,+TT=1,+T=3,+TC=1,-A=6124,-C=1,-TT=1,-G=1,-T=1, 1/1:8934,6,3,0,0,0:0.998546998994,0.000670615848888,0.000335307924444,0,0,0:+A=4382,+T=3,+G=5,-A=4552,-C=3,-T=1,-G=1, 1/1:10135,4,1,2,1,0:0.999014292755,0.000394282897979,9.85707244948e-05,0.00019714144899,9.85707244948e-05,0:+A=4966,+d1=1,+TC=1,+G=3,-A=5169,-d1=1,-C=1,-T=2,-G=1, 1/1:7342,4,1,0,0,0:0.999183451279,0.000544365813827,0.000136091453457,0,0,0:+A=3711,+T=1,-A=3631,-C=1,-G=4, +chrM 5000 . A G,T,C . . AC=28,27,10;AF=0.000251195421066,0.000242224156028,8.97126503808e-05 GT:AC:AF:NC 0/0:2,1,0:0.000190912562047,9.54562810233e-05,0:+A=3664,+T=1,-A=6809,-G=2, 0/0:7,1,2:0.000987724001693,0.000141103428813,0.000282206857627:+A=2337,+C=1,+T=1,+G=2,-A=4740,-C=1,-G=5, 0/0:2,3,0:0.00022411474675,0.000336172120126,0:+A=3194,+T=1,-A=5725,-T=2,-G=2, 0/0:1,1,0:0.000169865806013,0.000169865806013,0:+A=2022,+T=1,-A=3863,-G=1, +chrM 6000 . T C,A,TT,TCTTA,G . . AC=118279,31,7,1,1;AF=0.999459199108,0.000261950432218,5.91500975977e-05,8.45001394252e-06,8.45001394252e-06 GT:AC:AF:NC 1/1:9612,2,0,0,1:0.999168399168,0.0002079002079,0,0,0.00010395010395:+C=4845,+T=3,-A=2,-C=4767,-T=2,-G=1, 1/1:7771,0,1,0,0:0.99987133299,0,0.000128667009779,0,0:+C=3682,-C=4089,-TT=1, 1/1:9152,1,0,0,0:0.999563128003,0.000109217999126,0,0,0:+C=4640,+T=1,-A=1,-C=4512,-T=2, 1/1:5500,0,1,0,0:0.99981821487,0,0.000181785129976,0,0:+C=2707,-C=2793,-TT=1, +chrM 7000 . GTA TTA,CTA,ATA,GA,GACTACTA,G . . AC=177915,270,72,20,3,1;AF=0.997370840434,0.00151358866266,0.000403623643376,0.000112117678716,1.68176518073e-05,5.60588393578e-06 GT:AC:AF:NC 1/1:13745,26,4,1,0,0:0.997532476958,0.0018869293853,0.000290296828507,7.25742071268e-05,0,0:+d1=1,+C=17,+T=6216,+G=1,-A=4,-C=9,-T=7529,-G=2, 1/1:11392,16,0,1,0,0:0.997985107315,0.00140166447657,0,8.76040297854e-05,0,0:+d1=1,+C=7,+T=5104,+G=2,-C=9,-T=6288,-G=4, 1/1:13952,13,4,0,0,0:0.998068531368,0.000929966378139,0.000286143500966,0,0,0:+C=9,+T=6446,-A=4,-C=4,-T=7506,-G=10, 1/1:7740,12,1,0,0,0:0.998194480268,0.0015475883415,0.000128965695125,0,0,0:+C=8,+T=3471,-A=1,-C=4,-T=4269,-G=1, +chrM 8000 . TGA GGA,AGA,CGA,T . . AC=132300,33,7,2;AF=0.999350384482,0.000249271070959,5.28756817186e-05,1.51073376339e-05 GT:AC:AF:NC 1/1:10337,2,0,2:0.999226679555,0.000193330111165,0,0.000193330111165:+A=1,+d2=1,+T=1,+G=4832,-A=1,-d2=1,-T=3,-G=5505, 1/1:8347,3,0,0:0.999640718563,0.000359281437126,0,0:+A=2,+G=4258,-A=1,-G=4089, 1/1:10066,4,1,0:0.999205876514,0.000397061743101,9.92654357753e-05,0:+A=3,+C=1,+T=1,+G=5023,-A=1,-T=2,-G=5043, 1/1:4891,2,3,0:0.99897875817,0.000408496732026,0.000612745098039,0:+A=1,+C=1,+G=1746,-A=1,-C=2,-G=3145, +chrM 9000 . TACGC AACGC,GACGC,CACGC,T,TCGACGC . . AC=125089,77,46,8,2;AF=0.998842167463,0.000614849002667,0.000367312391204,6.38804158615e-05,1.59701039654e-05 GT:AC:AF:NC 1/1:10775,9,6,1,0:0.998517282921,0.000834028356964,0.000556018904643,9.26698174405e-05,0:+A=6306,+G=1,-A=4469,-C=6,-d4=1,-G=8, 1/1:8081,7,2,0,0:0.998764058831,0.00086515881844,0.00024718823384,0,0:+A=4777,+T=1,+G=3,-A=3304,-C=2,-G=4, 1/1:9475,4,3,2,0:0.998735111205,0.000421629598398,0.000316222198798,0.000210814799199,0:+A=5435,+G=1,-A=4040,-C=3,-d4=2,-T=3,-G=3, 1/1:5709,3,1,0,0:0.999299842465,0.000525118151584,0.000175039383861,0,0:+A=3234,+G=1,-A=2475,-C=1,-G=2, +#chrM 10000 . AGT GGT,TGT,AAGT,ACGT,A,CGT,AAGAGT,AACAGT,AT . . AC=88942,145,51,12,7,6,5,2,1;AF=0.997051734768,0.00162546942436,0.000571716832016,0.000134521607533,7.84709377277e-05,6.72608037666e-05,5.60506698055e-05,2.24202679222e-05,1.12101339611e-05 GT:AC:AF:NC 1/1:7632,14,4,3,1,1,0,0,0:0.996735013713,0.00182839232075,0.000522397805929,0.000391798354447,0.000130599451482,0.000130599451482,0,0,0:+A=1,+C=1,+AC=3,+G=2535,+T=11,-A=1,-AA=4,-d2=1,-T=3,-G=5097, 1/1:6053,5,2,1,0,0,1,0,0:0.998186015831,0.000824538258575,0.00032981530343,0.000164907651715,0,0,0.000164907651715,0,0:+AC=1,+T=5,+G=1674,-A=2,-AA=2,-G=4379,-AAGA=1, 1/1:7012,14,6,0,2,0,0,0,0:0.996730632552,0.00199004975124,0.000852878464819,0,0.000284292821606,0,0,0,0:+d2=1,+T=11,+G=2241,-A=1,-AA=6,-d2=1,-T=3,-G=4771, 1/1:4646,1,2,0,0,0,0,0,0:0.998710232158,0.000214961306965,0.000429922613929,0,0,0,0,0,0:+A=1,+T=1,+G=1337,-A=2,-AA=2,-G=3309, +#chrM 11000 . CCA ACA,GCA,TCA,CA,C . . AC=54,44,36,3,2;AF=0.000362343152385,0.00029524256861,0.00024156210159,2.01301751325e-05,1.3420116755e-05 GT:AC:AF:NC 0/0:7,3,1,0,0:0.000521920668058,0.000223680286311,7.45600954369e-05,0,0:+C=10229,+T=1,-A=7,-C=3172,-G=3, 0/0:3,4,2,0,0:0.000368233705659,0.000490978274211,0.000245489137106,0,0:+A=2,+C=6575,+T=2,+G=1,-A=1,-C=1563,-G=3, 0/0:9,3,2,0,0:0.000811615114077,0.000270538371359,0.000180358914239,0,0:+A=1,+C=7906,+T=2,-A=8,-C=3169,-G=3, 0/0:0,2,3,1,0:0,0.000220507166483,0.000330760749724,0.000110253583241,0:+d1=1,+C=7066,+T=2,-C=1998,-T=1,-G=2, +#chrM 12000 . AC CC,TC,GC,A . . AC=182586,84,23,1;AF=0.999146337459,0.000459664444955,0.000125860502785,5.47219577328e-06 GT:AC:AF:NC 1/1:13592,8,0,0:0.999264813998,0.000588148801647,0,0:+A=2,+C=7113,+T=5,-C=6479,-T=3, 1/1:9868,1,0,0:0.999797365755,0.000101317122594,0,0:+C=5182,-A=1,-C=4686,-T=1, 1/1:11289,6,3,0:0.998938147067,0.000530926466684,0.000265463233342,0:+A=2,+C=6123,+T=3,+G=1,-A=1,-C=5166,-T=3,-G=2, 1/1:9708,4,0,0:0.999382334775,0.000411776816965,0,0:+A=1,+C=4892,+T=2,-A=1,-C=4816,-T=2, +#chrM 13000 . A G,C,T . . AC=56125,9,6;AF=0.999358985773,0.000160253556739,0.000106835704492 GT:AC:AF:NC 1/1:3816,0,0:0.999214454046,0,0:+A=1,+G=1319,-A=2,-G=2497, 1/1:2814,1,0:0.999289772727,0.000355113636364,0:+G=986,-A=1,-C=1,-G=1828, 1/1:3479,0,0:0.999712643678,0,0:+A=1,+G=1285,-G=2194, 1/1:2679,0,0:1.0,0,0:+G=979,-G=1700, +#chrM 14000 . CTAA TTAA,GTAA,CAA,ATAA,C . . AC=147853,32,22,21,1;AF=0.997766290558,0.000215947740647,0.000148464071695,0.000141715704799,6.74836689521e-06 GT:AC:AF:NC 1/1:10992,3,2,4,0:0.996825972613,0.000272059490342,0.000181372993561,0.000362745987123,0:+A=1,+d1=1,+C=19,+T=7038,-A=3,-d1=1,-C=7,-T=3954,-G=3, 1/1:6853,2,2,0,0:0.997525473071,0.000291120815138,0.000291120815138,0,0:+d1=2,+C=9,+T=4609,-C=4,-T=2244,-G=2, 1/1:9277,2,0,5,0:0.998063474987,0.000215169445939,0,0.000537923614847,0:+A=2,+C=9,+T=5786,-A=3,-C=2,-T=3491,-G=2, 1/1:8125,1,3,2,0:0.997422047631,0.000122759636631,0.000368278909894,0.000245519273263,0:+d1=2,+C=9,+T=5241,-A=2,-d1=1,-C=6,-T=2884,-G=1, +#chrM 15000 . AA GA,TA,CA,A . . AC=132,29,29,4;AF=0.000611470633196,0.000134338245172,0.000134338245172,1.85294131272e-05 GT:AC:AF:NC 0/0:15,2,0,0:0.000956510649152,0.00012753475322,0,0:+A=8715,+T=1,+G=5,-A=6950,-T=1,-G=10, 0/0:13,1,0,0:0.00121985549404,9.38350380032e-05,0,0:+A=6102,+T=1,+G=6,-A=4541,-G=7, 0/0:1,2,5,1:7.86101721563e-05,0.000157220344313,0.000393050860781,7.86101721563e-05:+A=7378,+C=5,+T=2,-A=5334,-d1=1,-G=1, 0/0:6,0,0,0:0.000532765050613,0,0,0:+A=6215,+G=4,-A=5041,-G=2, +#chrM 16000 . AG GG,TG,CG,A . . AC=224750,90,3,1;AF=0.999244175707,0.000400142272808,1.33380757603e-05,4.44602525342e-06 GT:AC:AF:NC 1/1:15911,2,0,0:0.99949745587,0.000125636032414,0,0:+A=3,+T=1,+G=7156,-A=3,-T=1,-G=8755, 1/1:11311,2,0,0:0.999381516169,0.000176709666019,0,0:+A=3,+G=4973,-A=2,-T=2,-G=6338, 1/1:14157,2,0,0:0.999012066897,0.000141133300402,0,0:+A=5,+T=1,+G=6766,-A=7,-T=1,-G=7391, 1/1:11039,3,0,0:0.999185372918,0.000271542360608,0,0:+A=4,+T=1,+G=4959,-A=2,-T=2,-G=6080,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-mit.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,12 @@ +#SAMPLE CHR POS A C G T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. +S1 chrM 2000 8 14903 1 1 14913 1 C A 0.00054 +S3 chrM 2000 10 13175 1 6 13192 1 C A 0.00076 +S1 chrM 3000 27966 42 44 12 28064 0 A G 0.00157 +S2 chrM 3000 20472 16 36 4 20528 1 A G 0.00175 +S3 chrM 3000 30267 40 46 10 30363 0 A G 0.00152 +S4 chrM 3000 15494 17 29 3 15543 0 A G 0.00187 +S1 chrM 4000 12258 3 2 4 12267 1 A T 0.00033 +S1 chrM 7000 4 26 3 13745 13778 0 T C 0.00189 +S2 chrM 7000 0 16 6 11392 11414 1 T C 0.0014 +S3 chrM 7000 4 13 10 13952 13979 1 T C 0.00093 +S3 chrM 8000 4 1 10066 3 10074 1 G A 0.0004
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-mit.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,28 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-f 0.2 -c 5000 -H" +##fileDate=20130521 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_2536.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S1 S2 S3 S4 +chrM 1000 . ACTC TCTC,CCTC,GCTC,ATC,A . . AC=188032,682,8,8,1;AF=0.995310134556,0.00361003186568,4.23464148467e-05,4.23464148467e-05,5.29330185583e-06 GT:AC:AF:NC 1/1:12546,51,0,0,0:0.994688020297,0.00404344723698,0,0,0:+A=2,+C=20,+T=8020,-A=14,-C=31,-T=4526, 1/1:9970,38,1,1,0:0.995208624476,0.00379317228988,9.98203234178e-05,9.98203234178e-05,0:+A=1,+d1=1,+C=21,+T=6546,+G=1,-A=7,-C=17,-T=3424, 1/1:11372,35,0,2,0:0.995448179272,0.0030637254902,0,0.000175070028011,0:+d1=2,+C=10,+T=7199,-A=15,-C=25,-T=4173, 1/1:9115,39,0,0,0:0.995304651671,0.00425857174055,0,0,0:+C=13,+T=5817,-A=4,-C=26,-T=3298, +chrM 2000 . TAC CAC,AAC,GAC,T,TC,TGAC . . AC=209371,114,8,2,2,1;AF=0.999083807733,0.000543989158392,3.81746777819e-05,9.54366944547e-06,9.54366944547e-06,4.77183472273e-06 GT:AC:AF:NC 1/1:14903,8,1,0,0,0:0.999329444109,0.000536444712667,6.70555890834e-05,0,0,0:+A=1,+C=9095,+G=1,-A=7,-C=5808,-T=1, 1/1:11838,3,0,0,0,0:0.999324666554,0.000253250042208,0,0,0,0:+C=7242,+T=4,-A=3,-C=4596,-T=1, 1/1:13175,10,1,0,0,0:0.998711340206,0.000758035172832,7.58035172832e-05,0,0,0:+C=7933,+T=4,-A=10,-C=5242,-T=2,-G=1, 1/1:10921,3,2,0,0,1:0.999176578225,0.000274473924977,0.000182982616651,0,0,9.14913083257e-05:+A=1,+TG=1,+C=6732,+T=1,-A=2,-C=4189,-T=2,-G=2, +chrM 3000 . TC AC,GC,CC,T,TGC . . AC=386813,649,441,35,33;AF=0.996740346013,0.00167234421946,0.00113636949273,9.01880549786e-05,8.5034451837e-05 GT:AC:AF:NC 1/1:27966,44,42,5,6:0.996117542297,0.00156723063224,0.00149599287622,0.000178094390027,0.000213713268032:+A=17399,+d1=1,+C=7,+T=8,+G=22,-A=10567,-d1=4,-C=35,-G=22,-T=4,-TG=6, 1/1:20472,36,16,0,0:0.997272018706,0.00175370226033,0.000779423226812,0,0:+A=12535,+C=3,+T=2,+G=24,-A=7937,-C=13,-T=2,-G=12, 1/1:30267,46,40,5,3:0.996575680748,0.00151460274604,0.00131704586612,0.000164630733265,9.87784399592e-05:+A=18981,+d1=2,+C=7,+T=6,+G=29,-A=11286,-d1=3,-C=33,-G=17,-T=4,-TG=3, 1/1:15494,29,17,2,3:0.996526884487,0.00186519166452,0.00109338821713,0.000128633907898,0.000192950861847:+A=9254,+d1=2,+C=1,+T=2,+G=15,-A=6240,-TG=3,-C=16,-T=1,-G=14, +chrM 4000 . TA AA,GA,CA,T,TCA,TTA . . AC=138825,77,37,9,8,3;AF=0.998511134127,0.000553829334254,0.000266125783992,6.47332988089e-05,5.75407100524e-05,2.15777662696e-05 GT:AC:AF:NC 1/1:12258,2,3,1,1,2:0.998940591639,0.00016298590172,0.000244478852579,8.14929508598e-05,8.14929508598e-05,0.00016298590172:+A=6134,+d1=1,+C=2,+G=1,+TT=1,+T=3,+TC=1,-A=6124,-C=1,-TT=1,-G=1,-T=1, 1/1:8934,6,3,0,0,0:0.998546998994,0.000670615848888,0.000335307924444,0,0,0:+A=4382,+T=3,+G=5,-A=4552,-C=3,-T=1,-G=1, 1/1:10135,4,1,2,1,0:0.999014292755,0.000394282897979,9.85707244948e-05,0.00019714144899,9.85707244948e-05,0:+A=4966,+d1=1,+TC=1,+G=3,-A=5169,-d1=1,-C=1,-T=2,-G=1, 1/1:7342,4,1,0,0,0:0.999183451279,0.000544365813827,0.000136091453457,0,0,0:+A=3711,+T=1,-A=3631,-C=1,-G=4, +chrM 5000 . A G,T,C . . AC=28,27,10;AF=0.000251195421066,0.000242224156028,8.97126503808e-05 GT:AC:AF:NC 0/0:2,1,0:0.000190912562047,9.54562810233e-05,0:+A=3664,+T=1,-A=6809,-G=2, 0/0:7,1,2:0.000987724001693,0.000141103428813,0.000282206857627:+A=2337,+C=1,+T=1,+G=2,-A=4740,-C=1,-G=5, 0/0:2,3,0:0.00022411474675,0.000336172120126,0:+A=3194,+T=1,-A=5725,-T=2,-G=2, 0/0:1,1,0:0.000169865806013,0.000169865806013,0:+A=2022,+T=1,-A=3863,-G=1, +chrM 6000 . T C,A,TT,TCTTA,G . . AC=118279,31,7,1,1;AF=0.999459199108,0.000261950432218,5.91500975977e-05,8.45001394252e-06,8.45001394252e-06 GT:AC:AF:NC 1/1:9612,2,0,0,1:0.999168399168,0.0002079002079,0,0,0.00010395010395:+C=4845,+T=3,-A=2,-C=4767,-T=2,-G=1, 1/1:7771,0,1,0,0:0.99987133299,0,0.000128667009779,0,0:+C=3682,-C=4089,-TT=1, 1/1:9152,1,0,0,0:0.999563128003,0.000109217999126,0,0,0:+C=4640,+T=1,-A=1,-C=4512,-T=2, 1/1:5500,0,1,0,0:0.99981821487,0,0.000181785129976,0,0:+C=2707,-C=2793,-TT=1, +chrM 7000 . GTA TTA,CTA,ATA,GA,GACTACTA,G . . AC=177915,270,72,20,3,1;AF=0.997370840434,0.00151358866266,0.000403623643376,0.000112117678716,1.68176518073e-05,5.60588393578e-06 GT:AC:AF:NC 1/1:13745,26,4,1,0,0:0.997532476958,0.0018869293853,0.000290296828507,7.25742071268e-05,0,0:+d1=1,+C=17,+T=6216,+G=1,-A=4,-C=9,-T=7529,-G=2, 1/1:11392,16,0,1,0,0:0.997985107315,0.00140166447657,0,8.76040297854e-05,0,0:+d1=1,+C=7,+T=5104,+G=2,-C=9,-T=6288,-G=4, 1/1:13952,13,4,0,0,0:0.998068531368,0.000929966378139,0.000286143500966,0,0,0:+C=9,+T=6446,-A=4,-C=4,-T=7506,-G=10, 1/1:7740,12,1,0,0,0:0.998194480268,0.0015475883415,0.000128965695125,0,0,0:+C=8,+T=3471,-A=1,-C=4,-T=4269,-G=1, +chrM 8000 . TGA GGA,AGA,CGA,T . . AC=132300,33,7,2;AF=0.999350384482,0.000249271070959,5.28756817186e-05,1.51073376339e-05 GT:AC:AF:NC 1/1:10337,2,0,2:0.999226679555,0.000193330111165,0,0.000193330111165:+A=1,+d2=1,+T=1,+G=4832,-A=1,-d2=1,-T=3,-G=5505, 1/1:8347,3,0,0:0.999640718563,0.000359281437126,0,0:+A=2,+G=4258,-A=1,-G=4089, 1/1:10066,4,1,0:0.999205876514,0.000397061743101,9.92654357753e-05,0:+A=3,+C=1,+T=1,+G=5023,-A=1,-T=2,-G=5043, 1/1:4891,2,3,0:0.99897875817,0.000408496732026,0.000612745098039,0:+A=1,+C=1,+G=1746,-A=1,-C=2,-G=3145, +chrM 9000 . TACGC AACGC,GACGC,CACGC,T,TCGACGC . . AC=125089,77,46,8,2;AF=0.998842167463,0.000614849002667,0.000367312391204,6.38804158615e-05,1.59701039654e-05 GT:AC:AF:NC 1/1:10775,9,6,1,0:0.998517282921,0.000834028356964,0.000556018904643,9.26698174405e-05,0:+A=6306,+G=1,-A=4469,-C=6,-d4=1,-G=8, 1/1:8081,7,2,0,0:0.998764058831,0.00086515881844,0.00024718823384,0,0:+A=4777,+T=1,+G=3,-A=3304,-C=2,-G=4, 1/1:9475,4,3,2,0:0.998735111205,0.000421629598398,0.000316222198798,0.000210814799199,0:+A=5435,+G=1,-A=4040,-C=3,-d4=2,-T=3,-G=3, 1/1:5709,3,1,0,0:0.999299842465,0.000525118151584,0.000175039383861,0,0:+A=3234,+G=1,-A=2475,-C=1,-G=2, +#chrM 10000 . AGT GGT,TGT,AAGT,ACGT,A,CGT,AAGAGT,AACAGT,AT . . AC=88942,145,51,12,7,6,5,2,1;AF=0.997051734768,0.00162546942436,0.000571716832016,0.000134521607533,7.84709377277e-05,6.72608037666e-05,5.60506698055e-05,2.24202679222e-05,1.12101339611e-05 GT:AC:AF:NC 1/1:7632,14,4,3,1,1,0,0,0:0.996735013713,0.00182839232075,0.000522397805929,0.000391798354447,0.000130599451482,0.000130599451482,0,0,0:+A=1,+C=1,+AC=3,+G=2535,+T=11,-A=1,-AA=4,-d2=1,-T=3,-G=5097, 1/1:6053,5,2,1,0,0,1,0,0:0.998186015831,0.000824538258575,0.00032981530343,0.000164907651715,0,0,0.000164907651715,0,0:+AC=1,+T=5,+G=1674,-A=2,-AA=2,-G=4379,-AAGA=1, 1/1:7012,14,6,0,2,0,0,0,0:0.996730632552,0.00199004975124,0.000852878464819,0,0.000284292821606,0,0,0,0:+d2=1,+T=11,+G=2241,-A=1,-AA=6,-d2=1,-T=3,-G=4771, 1/1:4646,1,2,0,0,0,0,0,0:0.998710232158,0.000214961306965,0.000429922613929,0,0,0,0,0,0:+A=1,+T=1,+G=1337,-A=2,-AA=2,-G=3309, +#chrM 11000 . CCA ACA,GCA,TCA,CA,C . . AC=54,44,36,3,2;AF=0.000362343152385,0.00029524256861,0.00024156210159,2.01301751325e-05,1.3420116755e-05 GT:AC:AF:NC 0/0:7,3,1,0,0:0.000521920668058,0.000223680286311,7.45600954369e-05,0,0:+C=10229,+T=1,-A=7,-C=3172,-G=3, 0/0:3,4,2,0,0:0.000368233705659,0.000490978274211,0.000245489137106,0,0:+A=2,+C=6575,+T=2,+G=1,-A=1,-C=1563,-G=3, 0/0:9,3,2,0,0:0.000811615114077,0.000270538371359,0.000180358914239,0,0:+A=1,+C=7906,+T=2,-A=8,-C=3169,-G=3, 0/0:0,2,3,1,0:0,0.000220507166483,0.000330760749724,0.000110253583241,0:+d1=1,+C=7066,+T=2,-C=1998,-T=1,-G=2, +#chrM 12000 . AC CC,TC,GC,A . . AC=182586,84,23,1;AF=0.999146337459,0.000459664444955,0.000125860502785,5.47219577328e-06 GT:AC:AF:NC 1/1:13592,8,0,0:0.999264813998,0.000588148801647,0,0:+A=2,+C=7113,+T=5,-C=6479,-T=3, 1/1:9868,1,0,0:0.999797365755,0.000101317122594,0,0:+C=5182,-A=1,-C=4686,-T=1, 1/1:11289,6,3,0:0.998938147067,0.000530926466684,0.000265463233342,0:+A=2,+C=6123,+T=3,+G=1,-A=1,-C=5166,-T=3,-G=2, 1/1:9708,4,0,0:0.999382334775,0.000411776816965,0,0:+A=1,+C=4892,+T=2,-A=1,-C=4816,-T=2, +#chrM 13000 . A G,C,T . . AC=56125,9,6;AF=0.999358985773,0.000160253556739,0.000106835704492 GT:AC:AF:NC 1/1:3816,0,0:0.999214454046,0,0:+A=1,+G=1319,-A=2,-G=2497, 1/1:2814,1,0:0.999289772727,0.000355113636364,0:+G=986,-A=1,-C=1,-G=1828, 1/1:3479,0,0:0.999712643678,0,0:+A=1,+G=1285,-G=2194, 1/1:2679,0,0:1.0,0,0:+G=979,-G=1700, +#chrM 14000 . CTAA TTAA,GTAA,CAA,ATAA,C . . AC=147853,32,22,21,1;AF=0.997766290558,0.000215947740647,0.000148464071695,0.000141715704799,6.74836689521e-06 GT:AC:AF:NC 1/1:10992,3,2,4,0:0.996825972613,0.000272059490342,0.000181372993561,0.000362745987123,0:+A=1,+d1=1,+C=19,+T=7038,-A=3,-d1=1,-C=7,-T=3954,-G=3, 1/1:6853,2,2,0,0:0.997525473071,0.000291120815138,0.000291120815138,0,0:+d1=2,+C=9,+T=4609,-C=4,-T=2244,-G=2, 1/1:9277,2,0,5,0:0.998063474987,0.000215169445939,0,0.000537923614847,0:+A=2,+C=9,+T=5786,-A=3,-C=2,-T=3491,-G=2, 1/1:8125,1,3,2,0:0.997422047631,0.000122759636631,0.000368278909894,0.000245519273263,0:+d1=2,+C=9,+T=5241,-A=2,-d1=1,-C=6,-T=2884,-G=1, +#chrM 15000 . AA GA,TA,CA,A . . AC=132,29,29,4;AF=0.000611470633196,0.000134338245172,0.000134338245172,1.85294131272e-05 GT:AC:AF:NC 0/0:15,2,0,0:0.000956510649152,0.00012753475322,0,0:+A=8715,+T=1,+G=5,-A=6950,-T=1,-G=10, 0/0:13,1,0,0:0.00121985549404,9.38350380032e-05,0,0:+A=6102,+T=1,+G=6,-A=4541,-G=7, 0/0:1,2,5,1:7.86101721563e-05,0.000157220344313,0.000393050860781,7.86101721563e-05:+A=7378,+C=5,+T=2,-A=5334,-d1=1,-G=1, 0/0:6,0,0,0:0.000532765050613,0,0,0:+A=6215,+G=4,-A=5041,-G=2, +#chrM 16000 . AG GG,TG,CG,A . . AC=224750,90,3,1;AF=0.999244175707,0.000400142272808,1.33380757603e-05,4.44602525342e-06 GT:AC:AF:NC 1/1:15911,2,0,0:0.99949745587,0.000125636032414,0,0:+A=3,+T=1,+G=7156,-A=3,-T=1,-G=8755, 1/1:11311,2,0,0:0.999381516169,0.000176709666019,0,0:+A=3,+G=4973,-A=2,-T=2,-G=6338, 1/1:14157,2,0,0:0.999012066897,0.000141133300402,0,0:+A=5,+T=1,+G=6766,-A=7,-T=1,-G=7391, 1/1:11039,3,0,0:0.999185372918,0.000271542360608,0,0:+A=4,+T=1,+G=4959,-A=2,-T=2,-G=6080,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-nofilt.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,15 @@ +#SAMPLE CHR POS A C G T CVRG ALLELES MAJOR MINOR MINOR.FREQ.PERC. +THYROID chr1 246704250 29 0 0 0 29 1 A . 0.0 +THYROID chr1 246704257 0 0 0 71 71 1 T . 0.0 +THYROID chr1 246704268 104 0 0 0 104 1 A . 0.0 +THYROID chr1 246704269 0 0 0 105 105 1 T . 0.0 +THYROID chr1 246704363 0 72 3 0 75 0 C G 0.04 +THYROID chr1 246704437 5 130 0 0 135 0 C A 0.03704 +THYROID chr1 246707878 0 0 131 0 131 1 G . 0.0 +THYROID chr1 246714587 30 0 43 0 73 2 G A 0.41096 +THYROID chr1 246729215 1 0 1 88 90 0 T G 0.01111 +THYROID chr1 246729216 1 0 1 90 92 0 T G 0.01087 +THYROID chr1 246729378 16 7 0 0 23 0 A C 0.30435 +THYROID chr1 246729392 29 0 10 0 39 0 A G 0.25641 +THYROID chr7 91502881 0 0 0 26 26 1 T . 0.0 +THYROID chr7 91502897 7 36 0 0 43 0 C A 0.16279
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real-nofilt.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,27 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-H -r 1 -f 0 -c 0" +##comment="This is a test set of sites copied from real data." +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT THYROID +chr1 246704250 . A N . . AC=1;AF=0.0333333333333 GT:AC:AF:NC 0/0:1:0.0333333333333:+A=15,+N=1,-A=14, +chr1 246704257 . T N . . AC=2;AF=0.027397260274 GT:AC:AF:NC 0/0:2:0.027397260274:+T=52,+N=2,-T=19, +chr1 246704268 . A . . . AC=;AF= GT:AC:AF:NC 0/0:::+A=82,-A=22, +chr1 246704269 . T . . . AC=;AF= GT:AC:AF:NC 0/0:::+T=83,-T=22, +chr1 246704363 . C G,N . . AC=3,1;AF=0.0394736842105,0.0131578947368 GT:AC:AF:NC 0/0:3,1:0.0394736842105,0.0131578947368:+C=52,+G=3,+N=1,-C=20, +chr1 246704437 . C A,N . . AC=5,1;AF=0.0367647058824,0.00735294117647 GT:AC:AF:NC 0/0:5,1:0.0367647058824,0.00735294117647:+C=72,-A=5,-C=58,-N=1, +chr1 246707878 . G GA . . AC=6;AF=0.043795620438 GT:AC:AF:NC 0/0:6:0.043795620438:+G=90,-G=41,-GA=6, +chr1 246714587 . G A . . AC=30;AF=0.41095890411 GT:AC:AF:NC 0/1:30:0.41095890411:+A=10,+G=2,-A=20,-G=41, +chr1 246729215 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=7,-G=1, +chr1 246729216 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=9,-G=1, +chr1 246729378 . A C . . AC=7;AF=0.304347826087 GT:AC:AF:NC 0/1:7:0.304347826087:-A=16,-C=7, +chr1 246729392 . A G,AG . . AC=10,1;AF=0.25,0.025 GT:AC:AF:NC 0/1:10,1:0.25,0.025:-A=29,-AG=1,-G=10, +chr7 91502881 . TC T . . AC=3;AF=0.103448275862 GT:AC:AF:NC 0/0:3:0.103448275862:+T=14,-d1=3,-T=12, +chr7 91502897 . C A . . AC=7;AF=0.162790697674 GT:AC:AF:NC 0/0:7:0.162790697674:+A=7,+C=17,-C=19, \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real.csv.out Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,11 @@ +THYROID chr1 246704250 29 0 0 0 29 1 A . 0.0 +THYROID chr1 246704257 0 0 0 71 71 1 T . 0.0 +THYROID chr1 246704268 104 0 0 0 104 1 A . 0.0 +THYROID chr1 246704269 0 0 0 105 105 1 T . 0.0 +THYROID chr1 246704363 0 72 3 0 75 0 C G 0.04 +THYROID chr1 246704437 5 130 0 0 135 0 C A 0.03704 +THYROID chr1 246707878 0 0 131 0 131 1 G . 0.0 +THYROID chr1 246714587 30 0 43 0 73 2 G A 0.41096 +THYROID chr1 246729216 1 0 1 90 92 0 T G 0.01087 +THYROID chr7 91502881 0 0 0 26 26 1 T . 0.0 +THYROID chr7 91502897 7 36 0 0 43 0 C A 0.16279
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/real.vcf.in Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,27 @@ +##fileformat=VCFv4.1 +##comment="ARGS=-r 1 -f 1 -c 10" +##comment="This is a test set of sites copied from real data. It's meant to be run with -f 1 -c 10" +##fileDate=19700101 +##source=Dan +##reference=file:///scratch/dan/galaxy/galaxy-central/database/files/002/dataset_0000.dat +##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> +##FORMAT=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=AF,Number=.,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> +##FORMAT=<ID=NC,Number=.,Type=String,Description="Nucleotide and indel counts"> +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT THYROID +chr1 246704250 . A N . . AC=1;AF=0.0333333333333 GT:AC:AF:NC 0/0:1:0.0333333333333:+A=15,+N=1,-A=14, +chr1 246704257 . T N . . AC=2;AF=0.027397260274 GT:AC:AF:NC 0/0:2:0.027397260274:+T=52,+N=2,-T=19, +chr1 246704268 . A . . . AC=;AF= GT:AC:AF:NC 0/0:::+A=82,-A=22, +chr1 246704269 . T . . . AC=;AF= GT:AC:AF:NC 0/0:::+T=83,-T=22, +chr1 246704363 . C G,N . . AC=3,1;AF=0.0394736842105,0.0131578947368 GT:AC:AF:NC 0/0:3,1:0.0394736842105,0.0131578947368:+C=52,+G=3,+N=1,-C=20, +chr1 246704437 . C A,N . . AC=5,1;AF=0.0367647058824,0.00735294117647 GT:AC:AF:NC 0/0:5,1:0.0367647058824,0.00735294117647:+C=72,-A=5,-C=58,-N=1, +chr1 246707878 . G GA . . AC=6;AF=0.043795620438 GT:AC:AF:NC 0/0:6:0.043795620438:+G=90,-G=41,-GA=6, +chr1 246714587 . G A . . AC=30;AF=0.41095890411 GT:AC:AF:NC 0/1:30:0.41095890411:+A=10,+G=2,-A=20,-G=41, +chr1 246729215 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=7,-G=1, +chr1 246729216 . T G,A . . AC=1,1;AF=0.0111111111111,0.0111111111111 GT:AC:AF:NC 0/0:1,1:0.0111111111111,0.0111111111111:+A=1,+T=81,-T=9,-G=1, +chr1 246729378 . A C . . AC=7;AF=0.304347826087 GT:AC:AF:NC 0/1:7:0.304347826087:-A=16,-C=7, +chr1 246729392 . A G,AG . . AC=10,1;AF=0.25,0.025 GT:AC:AF:NC 0/1:10,1:0.25,0.025:-A=29,-AG=1,-G=10, +chr7 91502881 . TC T . . AC=3;AF=0.103448275862 GT:AC:AF:NC 0/0:3:0.103448275862:+T=14,-d1=3,-T=12, +chr7 91502897 . C A . . AC=7;AF=0.162790697674 GT:AC:AF:NC 0/0:7:0.162790697674:+A=7,+C=17,-C=19, \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/run-tests.py Thu Sep 12 11:34:23 2013 -0400 @@ -0,0 +1,54 @@ +#!/usr/bin/env python +import os +import sys +import subprocess + +DATASETS = [ + 'artificial', + 'artificial-samples', + 'artificial-nofilt', + 'real', + 'real-mit', + 'real-mit-s', + 'real-nofilt', +] +IN_EXT = '.vcf.in' +OUT_EXT = '.csv.out' +ARGS_KEY = '##comment="ARGS=' + +def main(): + + test_dir = os.path.dirname(os.path.relpath(sys.argv[0])) + if test_dir: + test_dir += os.sep + + for dataset in DATASETS: + infile = test_dir+dataset+IN_EXT + outfile = test_dir+dataset+OUT_EXT + + if not os.path.exists(infile): + sys.stderr.write("Error: file not found: "+infile+"\n") + continue + if not os.path.exists(outfile): + sys.stderr.write("Error: file not found: "+outfile+"\n") + continue + + options = read_options(infile) + script_cmd = 'allele-counts.py '+options+' -i '+infile + bash_cmd = 'diff '+outfile+' <('+script_cmd+')' + # print infile+":" + print script_cmd + subprocess.call(['bash', '-c', bash_cmd]) + + +def read_options(infile): + with open(infile, 'r') as infilehandle: + for line in infilehandle: + line.strip() + if ARGS_KEY == line[:len(ARGS_KEY)]: + return line[len(ARGS_KEY):-2] + return '' + + +if __name__ == '__main__': + main() \ No newline at end of file