GTF to BEDGraph converter.

This commit is contained in:
Jeremy Goecks
2010-04-22 16:05:08 -04:00
parent 3d9c8baaad
commit f8eabb60fd
3 changed files with 153 additions and 0 deletions
+1
View File
@@ -79,6 +79,7 @@
<tool file="fastx_toolkit/fastq_to_fasta.xml" />
<tool file="filters/wiggle_to_simple.xml" />
<tool file="filters/sff_extractor.xml" />
<tool file="filters/gtf2bedgraph.xml" />
</section>
<section name="Extract Features" id="features">
<tool file="filters/ucsc_gene_bed_to_exon_bed.xml" />
+79
View File
@@ -0,0 +1,79 @@
<tool id="gtf2bedgraph" name="GTF-to-BEDGraph">
<description>converter</description>
<command interpreter="python">gtf_to_bedgraph_converter.py $input $out_file1 $attribute_name</command>
<inputs>
<param format="gtf" name="input" type="data" label="Convert this query"/>
<param name="attribute_name" type="text" label="Attribute to Use for Value"/>
</inputs>
<outputs>
<data format="bedgraph" name="out_file1" />
</outputs>
<tests>
<test>
<param name="input" value="gtf2bedgraph_in.gtf" ftype="gtf"/>
<param name="attribute_name" value="FPKM"/>
<output name="out_file1" file="gtf2bedgraph_out.bedgraph" ftype="bedgraph"/>
</test>
</tests>
<help>
**What it does**
This tool converts data from GTF format to BEDGraph format (scroll down for format description).
--------
**Example**
The following data in GFF format::
chr22 GeneA enhancer 10000000 10001000 500 + . gene_id "GeneA"; transcript_id "TranscriptAlpha"; FPKM "2.75"; frac "1.000000";
chr22 GeneA promoter 10010000 10010100 900 + . gene_id "GeneA"; transcript_id "TranscriptsAlpha"; FPKM "2.25"; frac "1.000000";
using the attribute name 'FPKM' will be converted to BEDGraph (**note** that 1 is subtracted from the start coordinate)::
chr22 9999999 10001000 2.75
chr22 10009999 10010100 2.25
------
.. class:: infomark
**About formats**
**GTF format** Gene Transfer Format is a format for describing genes and other features associated with DNA, RNA and Protein sequences. GTF lines have nine tab-separated fields::
1. seqname - Must be a chromosome or scaffold.
2. source - The program that generated this feature.
3. feature - The name of this type of feature. Some examples of standard feature types are "CDS", "start_codon", "stop_codon", and "exon".
4. start - The starting position of the feature in the sequence. The first base is numbered 1.
5. end - The ending position of the feature (inclusive).
6. score - A score between 0 and 1000. If there is no score value, enter ".".
7. strand - Valid entries include '+', '-', or '.' (for don't know/care).
8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'.
9. group - The group field is a list of attributes. Each attribute consists of a type/value pair. Attributes must end in a semi-colon, and be separated from any following attribute by exactly one space. The attribute list must begin with the two mandatory attributes: (i) gene_id value - A globally unique identifier for the genomic source of the sequence and (ii) transcript_id value - A globally unique identifier for the predicted transcript.
**BEDGraph format**
The bedGraph format is line-oriented. Bedgraph data are preceeded by a track definition line, which adds a number of options for controlling the default display of this track.
For the track definition line, all options are placed in a single line separated by spaces:
track type=bedGraph name=track_label description=center_label
visibility=display_mode color=r,g,b altColor=r,g,b
priority=priority autoScale=on|off alwaysZero=on|off
gridDefault=on|off maxHeightPixels=max:default:min
graphType=bar|points viewLimits=lower:upper
yLineMark=real-value yLineOnOff=on|off
windowingFunction=maximum|mean|minimum smoothingWindow=off|2-16
The track type is REQUIRED, and must be bedGraph:
type=bedGraph
Following the track definition line are the track data in four column BED format::
chromA chromStartA chromEndA dataValueA
chromB chromStartB chromEndB dataValueB
</help>
</tool>
@@ -0,0 +1,73 @@
#!/usr/bin/env python
import os, sys, tempfile
assert sys.version_info[:2] >= ( 2, 4 )
def __main__():
# Read parms.
input_name = sys.argv[1]
output_name = sys.argv[2]
attribute_name = sys.argv[3]
# Create temp file.
tmp_name = tempfile.NamedTemporaryFile().name
# Do conversion.
skipped_lines = 0
first_skipped_line = 0
out = open( tmp_name, 'w' )
# Write track definition line.
out.write( "track type=bedGraph\n")
# Write track data to temporary file.
i = 0
for i, line in enumerate( file( input_name ) ):
line = line.rstrip( '\r\n' )
if line and not line.startswith( '#' ):
try:
elems = line.split( '\t' )
start = str( int( elems[3] ) - 1 ) # GTF coordinates are 1-based, BedGraph are 0-based.
strand = elems[6]
if strand not in ['+', '-']:
strand = '+'
attributes_list = elems[8].split(";")
attributes = {}
for name_value_pair in attributes_list:
pair = name_value_pair.strip().split(" ")
name = pair[0].strip()
if name == '':
continue
# Need to strip double quote from values
value = pair[1].strip(" \"")
attributes[name] = value
value = attributes[ attribute_name ]
# GTF format: chrom source, name, chromStart, chromEnd, score, strand, frame, attributes.
# BedGraph format: chrom, chromStart, chromEnd, value
out.write( "%s\t%s\t%s\t%s\n" %( elems[0], start, elems[4], value ) )
except:
skipped_lines += 1
if not first_skipped_line:
first_skipped_line = i + 1
else:
skipped_lines += 1
if not first_skipped_line:
first_skipped_line = i + 1
out.close()
# Sort tmp file to create bedgraph file; sort by chromosome name and chromosome start.
cmd = "sort -k1,1 -k2,2n < %s > %s" % ( tmp_name, output_name )
try:
os.system(cmd)
os.remove(tmp_name)
except Exception, ex:
sys.stderr.write( "%s\n" % ex )
sys.exit(1)
info_msg = "%i lines converted to BEDGraph. " % ( i + 1 - skipped_lines )
if skipped_lines > 0:
info_msg += "Skipped %d blank/comment/invalid lines starting with line #%d." %( skipped_lines, first_skipped_line )
print info_msg
if __name__ == "__main__": __main__()