Mercurial > repos > greg > extract_genomic_dna
changeset 0:cff5b7c9be55 draft
Uploaded
| author | greg |
|---|---|
| date | Thu, 14 Jan 2016 07:55:22 -0500 |
| parents | |
| children | 311febbd33d6 |
| files | extract_genomic_dna.py extract_genomic_dna.xml tool_dependencies.xml |
| diffstat | 3 files changed, 477 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extract_genomic_dna.py Thu Jan 14 07:55:22 2016 -0500 @@ -0,0 +1,264 @@ +#!/usr/bin/env python +import argparse +import os +import subprocess +import sys +import tempfile + +import bx.seq.nib +import bx.seq.twobit +from bx.intervals.io import Header, Comment +from galaxy.datatypes.util import gff_util +from galaxy.tools.util.galaxyops import parse_cols_arg + + +def convert_to_twobit(reference_genome): + """ + Create 2bit file history fasta dataset. + """ + try: + seq_path = tempfile.NamedTemporaryFile(dir=".").name + cmd = "faToTwoBit %s %s" % (reference_genome, seq_path) + tmp_name = tempfile.NamedTemporaryFile(dir=".").name + tmp_stderr = open(tmp_name, 'wb') + proc = subprocess.Popen(args=cmd, shell=True, stderr=tmp_stderr.fileno()) + returncode = proc.wait() + tmp_stderr.close() + if returncode != 0: + # Get stderr, allowing for case where it's very large. + tmp_stderr = open(tmp_name, 'rb') + stderr = '' + buffsize = 1048576 + try: + while True: + stderr += tmp_stderr.read(buffsize) + if not stderr or len(stderr) % buffsize != 0: + break + except OverflowError: + pass + tmp_stderr.close() + os.remove(tmp_name) + stop_err(stderr) + return seq_path + except Exception, e: + stop_err('Error running faToTwoBit. ' + str(e)) + + +def get_lines(feature): + # Get feature's line(s). + if isinstance(feature, gff_util.GFFFeature): + return feature.lines() + else: + return [feature.rstrip('\r\n')] + + +def reverse_complement(s): + complement_dna = {"A": "T", "T": "A", "C": "G", "G": "C", "a": "t", "t": "a", "c": "g", "g": "c", "N": "N", "n": "n"} + reversed_s = [] + for i in s: + reversed_s.append(complement_dna[i]) + reversed_s.reverse() + return "".join(reversed_s) + + +def stop_err(msg): + sys.stderr.write(msg) + sys.exit(1) + +parser = argparse.ArgumentParser() +parser.add_option('--input_format', dest='input_format', help="Input dataset format") +parser.add_option('--input', dest='input', help="Input dataset") +parser.add_option('--dbkey', dest='dbkey', help="Input dataset genome build") +parser.add_option('--interpret_features', dest='interpret_features', default=None, help="Interpret features if input format is gff") +parser.add_option('--columns', dest='columns', help="Columns to use in input file") +parser.add_option('--reference_genome_source', dest='reference_genome_source', help="Source of reference genome file") +parser.add_option('--reference_genome', dest='reference_genome', help="Reference genome file") +parser.add_option('--output_format', dest='output_format', help="Output format") +parser.add_option('--output', dest='output', help="Output dataset") +args = parser.parse_args() + +input_is_gff = args.input_format == 'gff' +interpret_features = args.interpret_features == "yes" +if len(args.cols.split(',')) == 5: + # Bed file. + chrom_col, start_col, end_col, strand_col, name_col = parse_cols_arg(args.cols) +else: + # Gff file. + chrom_col, start_col, end_col, strand_col = parse_cols_arg(args.cols) + name_col = False + +if args.reference_genome_source == "history": + seq_path = convert_to_twobit(args.reference_genome) +else: + seq_path = args.reference_genome +seq_dir = os.path.split(seq_path)[0] + +includes_strand_col = strand_col >= 0 +strand = None +nibs = {} +skipped_lines = 0 +first_invalid_line = 0 +invalid_lines = [] +warnings = [] +warning = '' +twobitfile = None +line_count = 1 +file_iterator = open(args.input) +if interpret_features: + file_iterator = gff_util.GFFReaderWrapper(file_iterator, fix_strand=False) +out = open(args.output, 'wt') + +for feature in file_iterator: + # Ignore comments, headers. + if isinstance(feature, (Header, Comment)): + line_count += 1 + continue + name = "" + if interpret_features: + # Processing features. + gff_util.convert_gff_coords_to_bed(feature) + chrom = feature.chrom + start = feature.start + end = feature.end + strand = feature.strand + else: + # Processing lines, either interval or GFF format. + line = feature.rstrip('\r\n') + if line and not line.startswith("#"): + fields = line.split('\t') + try: + chrom = fields[chrom_col] + start = int(fields[start_col]) + end = int(fields[end_col]) + if name_col: + name = fields[name_col] + if input_is_gff: + start, end = gff_util.convert_gff_coords_to_bed([start, end]) + if includes_strand_col: + strand = fields[strand_col] + except: + warning = "Invalid chrom, start or end column values. " + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + if start > end: + warning = "Invalid interval, start '%d' > end '%d'. " % (start, end) + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + if strand not in ['+', '-']: + strand = '+' + sequence = '' + else: + continue + # Open sequence file and get sequence for feature/interval. + if os.path.exists("%s/%s.nib" % (seq_dir, chrom)): + if chrom in nibs: + nib = nibs[chrom] + else: + nibs[chrom] = nib = bx.seq.nib.NibFile(file("%s/%s.nib" % (seq_path, chrom))) + try: + sequence = nib.get(start, end - start) + except Exception, e: + warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'. " % (start, end - start, args.dbkey) + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + elif os.path.isfile(os.path.join(seq_dir, '%s.2bit' % args.dbkey)): + if not(twobitfile): + twobitfile = bx.seq.twobit.TwoBitFile(file(seq_path)) + try: + if interpret_features: + # Create sequence from intervals within a feature. + sequence = '' + for interval in feature.intervals: + sequence += twobitfile[interval.chrom][interval.start:interval.end] + else: + sequence = twobitfile[chrom][start:end] + except: + warning = "Unable to fetch the sequence from '%d' to '%d' for chrom '%s'. " % (start, end - start, chrom) + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + else: + warning = "Chromosome by name '%s' was not found for build '%s'. " % (chrom, args.dbkey) + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + if sequence == '': + warning = "Chrom: '%s', start: '%s', end: '%s' is either invalid or not present in build '%s'. " % (chrom, start, end, args.dbkey) + warnings.append(warning) + if not invalid_lines: + invalid_lines = get_lines(feature) + first_invalid_line = line_count + skipped_lines += len(invalid_lines) + continue + if includes_strand_col and strand == "-": + sequence = reverse_complement(sequence) + if args.output_format == "fasta": + l = len(sequence) + c = 0 + if input_is_gff: + start, end = gff_util.convert_bed_coords_to_gff([start, end]) + fields = [args.dbkey, str(chrom), str(start), str(end), strand] + meta_data = "_".join(fields) + if name.strip(): + out.write(">%s %s\n" % (meta_data, name)) + else: + out.write(">%s\n" % meta_data) + while c < l: + b = min(c + 50, l) + out.write("%s\n" % str(sequence[c:b])) + c = b + else: + # output_format == "interval". + if interpret_features: + meta_data = "\t".join([feature.chrom, + "galaxy_extract_genomic_dna", + "interval", + str(feature.start), + str(feature.end), + feature.score, + feature.strand, + ".", + gff_util.gff_attributes_to_str(feature.attributes, "GTF")]) + else: + # Where is fields being set here? + meta_data = "\t".join(fields) + if input_is_gff: + format_str = "%s seq \"%s\";\n" + else: + format_str = "%s\t%s\n" + out.write(format_str % (meta_data, str(sequence))) + # Update line count. + if isinstance(feature, gff_util.GFFFeature): + line_count += len(feature.intervals) + else: + line_count += 1 +out.close() + +if warnings: + warn_msg = "%d warnings, 1st is: " % len(warnings) + warn_msg += warnings[0] + print warn_msg +if skipped_lines: + # Error message includes up to the first 10 skipped lines. + print 'Skipped %d invalid lines, 1st is #%d, "%s"' % (skipped_lines, first_invalid_line, '\n'.join(invalid_lines[:10])) + +if args.reference_genome_source == "history": + os.remove(seq_path)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extract_genomic_dna.xml Thu Jan 14 07:55:22 2016 -0500 @@ -0,0 +1,207 @@ +<tool id="Extract genomic DNA 1" name="Extract Genomic DNA" version="3.0.0"> + <description>using coordinates from assembled/unassembled genomes</description> + <command> + <![CDATA[ + #set input_format $input_format_cond.input_format + #set input $input_format_cond.input + #set dbkey = $input.metadata.dbkey + #set datatype = $input.datatype + mkdir -p output_dir && + python $__tool_directory__/extract_genomic_dna.py + --input_format $input_format + --input "$input" + --dbkey $dbkey + #if str($input_format) == "gff": + --interpret_features $input_format_cond.interpret_features + #end if + #if isinstance($datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__): + --columns "1,4,5,7" + #else: + --columns "${input.metadata.chromCol},${input.metadata.startCol},${input.metadata.endCol},${input.metadata.strandCol},${input.metadata.nameCol}" + #end if + --reference_genome_source $reference_genome_cond.reference_genome_source + #if str($reference_genome_cond.reference_genome_source) == "cached" + --reference_genome $reference_genome_cond.reference_genome.fields.path + #else: + --reference_genome $reference_genome_cond.reference_genome + #end if + --output_format $output_format + --output $output + ]]> + </command> + <inputs> + <conditional name="input_format_cond"> + <param name="input_format" type="select" label="Input file format"> + <option value="gff" selected="True">Gff</option> + <option value="interval">Interval</option> + </param> + <when value="gff"> + <param name="input" type="data" format="gff" label="Fetch sequences for intervals in"> + <validator type="unspecified_build" /> + </param> + <param name="interpret_features" type="select" label="Interpret features when possible"> + <option value="yes">Yes</option> + <option value="no">No</option> + </param> + </when> + <when value="interval"> + <param name="input" type="data" format="interval" label="Fetch sequences for intervals in"> + <validator type="unspecified_build" /> + </param> + </when> + </conditional> + <conditional name="reference_genome_cond"> + <param name="reference_genome_source" type="select" label="Choose the source for the reference genome"> + <option value="cached">locally cached</option> + <option value="history">from history</option> + </param> + <when value="cached"> + <param name="reference_genome" type="select" label="Using reference genome"> + <options from_data_table="alignseq_seq"> + <filter type="data_meta" key="dbkey" ref="input" column="dbkey"/> + </options> + <validator type="no_options" message="A built-in reference genome is not available for the build associated with the selected input file"/> + </param> + </when> + <when value="history"> + <param name="reference_genome" type="data" format="fasta" label="Using reference genome"> + <options> + <filter type="data_meta" key="dbkey" ref="input_bam" /> + </options> + <validator type="no_options" message="The current history does not include a fasta dataset with the build associated with the selected input file"/> + </param> + </when> + </conditional> + <param name="output_format" type="select" label="Select output format"> + <option value="fasta" selected="True">fasta</option> + <option value="interval">interval</option> + </param> + </inputs> + <outputs> + <data name="output" format="gff"> + <change_format> + <when output_format="interval" format="interval" /> + </change_format> + </data> + </outputs> + <tests> + <test> + <param name="input" value="1.bed" dbkey="hg17" ftype="bed" /> + <param name="interpret_features" value="yes"/> + <param name="index_source" value="cached"/> + <param name="out_format" value="fasta"/> + <output name="out_file1"> + <assert_contents> + <!-- First few lines... --> + <has_text text=">hg17_chr1_147962192_147962580_- CCDS989.1_cds_0_0_chr1_147962193_r" /> + <has_text text="ACTTGATCCTGCTCCCTCGGTGTCTGCATTGACTCCTCATGCTGGGACTG" /> + <has_text text="GACCCGTCAACCCCCCTGCTCGCTGCTCACGTACCTTCATCACTTTTAGT" /> + <has_text text="GATGATGCAACTTTCGAGGAATGGTTCCCCCAAGGGCGGCCCCCAAAAGT" /> + <!-- Last few lines... --> + <has_text text="GCTGTGGCACAGAACATGGACTCTGTGTTTAAGGAGCTCTTGGGAAAGAC" /> + <has_text text="CTCTGTCCGCCAGGGCCTTGGGCCAGCATCTACCACCTCTCCCAGTCCTG" /> + <has_text text="GGCCCCGAAGCCCAAAGGCCCCGCCCAGCAGCCGCCTGGGCAGGAACAAA" /> + <has_text text="GGCTTCTCCCGGGGCCCTGGGGCCCCAGCCTCACCCTCAGCTTCCCACCC" /> + <has_text text="CCAGGGCCTAGACACGACCCCCAAGCCACACTGA" /> + </assert_contents> + </output> + </test> + <test> + <param name="input" value="droPer1.bed" dbkey="droPer1" ftype="bed" /> + <param name="interpret_features" value="yes"/> + <param name="index_source" value="cached"/> + <param name="out_format" value="fasta"/> + <output name="out_file1" file="extract_genomic_dna_out2.fasta" /> + </test> + <test> + <param name="input" value="1.bed" dbkey="hg17" ftype="bed" /> + <param name="interpret_features" value="yes"/> + <param name="index_source" value="cached"/> + <param name="out_format" value="interval"/> + <output name="out_file1" file="extract_genomic_dna_out3.interval" /> + </test> + <!-- Test GFF file support. --> + <test> + <param name="input" value="gff_filter_by_attribute_out1.gff" dbkey="mm9" ftype="gff" /> + <param name="interpret_features" value="no"/> + <param name="index_source" value="cached"/> + <param name="out_format" value="interval"/> + <output name="out_file1" file="extract_genomic_dna_out4.gff" /> + </test> + <test> + <param name="input" value="gff_filter_by_attribute_out1.gff" dbkey="mm9" ftype="gff" /> + <param name="interpret_features" value="no"/> + <param name="out_format" value="fasta"/> + <param name="index_source" value="cached"/> + <output name="out_file1" file="extract_genomic_dna_out5.fasta" /> + </test> + <!-- Test custom sequences support and GFF feature interpretation. --> + <test> + <param name="input" value="cufflinks_out1.gtf" dbkey="mm9" ftype="gff" /> + <param name="interpret_features" value="no"/> + <param name="index_source" value="history"/> + <param name="ref_file" value="tophat_in1.fasta"/> + <param name="out_format" value="fasta"/> + <output name="out_file1" file="extract_genomic_dna_out6.fasta" /> + </test> + <test> + <param name="input" value="cufflinks_out1.gtf" dbkey="mm9" ftype="gff" /> + <param name="interpret_features" value="yes"/> + <param name="index_source" value="history"/> + <param name="ref_file" value="tophat_in1.fasta"/> + <param name="out_format" value="fasta"/> + <output name="out_file1" file="extract_genomic_dna_out7.fasta" /> + </test> + </tests> + <help> + +.. class:: warningmark + +The following will cause a line from the input dataset to be skipped and a warning generated. + + - Sequences that fall outside of the range of a line's start and end coordinates. + - Chromosome start or end coordinates that are invalid for the specified build. + +----- + +**What it does** + +This tool uses coordinate, strand, and build information to fetch genomic DNA from gff data, producing fasta data. + +----- + +**Example** + +If the input dataset is:: + + chr7 127475281 127475310 NM_000230 0 + + chr7 127485994 127486166 NM_000230 0 + + chr7 127486011 127486166 D49487 0 + + +Extracting sequences returns:: + + >hg17_chr7_127475281_127475310_+ NM_000230 + GTAGGAATCGCAGCGCCAGCGGTTGCAAG + >hg17_chr7_127485994_127486166_+ NM_000230 + GCCCAAGAAGCCCATCCTGGGAAGGAAAATGCATTGGGGAACCCTGTGCG + GATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATC + CAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAG + GATCAATGACATTTCACACACG + >hg17_chr7_127486011_127486166_+ D49487 + TGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGG + CCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGA + CACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCAC + ACACG + + </help> + <citations> + <citation type="bibtex"> + @unpublished{None, + author = {}, + title = {None}, + year = {None}, + eprint = {None}, + url = {http://www.bx.psu.edu/~anton/labSite/} + }</citation> + </citations> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Thu Jan 14 07:55:22 2016 -0500 @@ -0,0 +1,6 @@ +<?xml version="1.0"?> +<tool_dependency> + <package name="faToTwoBit" version="35x1"> + <repository changeset_revision="d2839a511c5a" name="package_fatotwobit_35x1" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" /> + </package> +</tool_dependency>