mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Extract genomic DNA will now result in warnings rather than error conditions in most cases.
This commit is contained in:
Greg Von Kuster
@@ -96,7 +96,7 @@ def main():
|
||||
# No overlap with any partition? For now throw this since the
|
||||
# partitions tile the encode regions completely, indicate an interval
|
||||
# that does not even overlap an encode region
|
||||
warning = "warning: Interval (%s, %d, %d) does not overlap any partition" % ( chr, start, end ) + ", line[" + str( line_count ) + "]"
|
||||
warning = "warning: Interval (%s, %d, %d) does not overlap any partition" % ( chr, start, end ) + ", line[" + str( line_count ) + "]. "
|
||||
warnings.append( warning )
|
||||
name = "no_overlap"
|
||||
score = 0
|
||||
@@ -112,14 +112,8 @@ def main():
|
||||
in_file.close()
|
||||
|
||||
if warnings:
|
||||
warn_msg = "This tool is useful on ENCODE regions only."
|
||||
if len( warnings ) > 2:
|
||||
warn_msg += "More than 2 warnings: "
|
||||
for warning in warnings[0:2]:
|
||||
warn_msg += warning + ", "
|
||||
else:
|
||||
for warning in warnings:
|
||||
warn_msg += warning + ", "
|
||||
warn_msg = "Total of %d warnings, 1st is: " % len( warnings )
|
||||
warn_msg += warnings[0]
|
||||
print warn_msg
|
||||
if skipped_lines:
|
||||
print "Skipped %d invalid lines starting at line # %d: %s" % ( skipped_lines, first_invalid_line, invalid_line )
|
||||
|
||||
@@ -80,26 +80,20 @@ def __main__():
|
||||
pass
|
||||
dbkey = sys.argv[7]
|
||||
output_format = sys.argv[8]
|
||||
|
||||
# TODO: is this still necessary? If so, let's get it fixed!
|
||||
if (re.search("^mm\d$", dbkey)):
|
||||
dbkey = "musMus" + dbkey[-1]
|
||||
if (re.search("^rn\d$", dbkey)):
|
||||
dbkey = "ratNor" + dbkey[-1]
|
||||
|
||||
nibs = {}
|
||||
twobits = {}
|
||||
nib_path = check_nib_file( dbkey )
|
||||
twobit_path = check_twobit_file( dbkey )
|
||||
if not os.path.exists( nib_path ) and not os.path.exists( twobit_path ):
|
||||
# If this occurs, we need to fix the metadata validator.
|
||||
stop_err( "No sequences are available for %s, request them by reporting this error." % dbkey )
|
||||
stop_err( "No sequences are available for '%s', request them by reporting this error." % dbkey )
|
||||
|
||||
skipped_lines = 0
|
||||
first_invalid_line = 0
|
||||
invalid_line = ''
|
||||
fout = open( output_filename, "w" )
|
||||
err_msg = ''
|
||||
warnings = []
|
||||
warning = ''
|
||||
|
||||
for i, line in enumerate( open( input_filename ) ):
|
||||
line = line.rstrip( '\r\n' )
|
||||
@@ -109,64 +103,91 @@ def __main__():
|
||||
chrom = fields[chrom_col]
|
||||
start = int( fields[start_col] )
|
||||
end = int( fields[end_col] )
|
||||
if includes_strand_col:
|
||||
strand = fields[strand_col]
|
||||
if strand not in ['+', '-']:
|
||||
strand = "+"
|
||||
|
||||
sequence = ''
|
||||
if os.path.exists( "%s/%s.nib" % ( nib_path, chrom) ):
|
||||
if chrom in nibs:
|
||||
nib = nibs[chrom]
|
||||
else:
|
||||
nibs[chrom] = nib = bx.seq.nib.NibFile( file( "%s/%s.nib" % ( nib_path, chrom ) ) )
|
||||
try:
|
||||
sequence = nib.get( start, end-start )
|
||||
except:
|
||||
err_msg = "Unable to fetch the sequence from %d to %d from %s." %( start, end-start, nib_path )
|
||||
break
|
||||
elif os.path.exists( twobit_path ):
|
||||
if chrom in twobits:
|
||||
t = twobits[chrom]
|
||||
else:
|
||||
twobits[chrom] = t = bx.seq.twobit.TwoBitFile( file( twobit_path ) )
|
||||
try:
|
||||
sequence = t[chrom][start:end]
|
||||
except:
|
||||
err_msg = "Unable to fetch the sequence from %d to %d from %s." %( start, end-start, twobit_path )
|
||||
break
|
||||
else:
|
||||
err_msg = "Sequence %s was not found for build %s. Most likely your data lists the wrong chromosome number for this organism. Check your build selection." % ( chrom, dbkey )
|
||||
break
|
||||
|
||||
if not sequence:
|
||||
err_msg = "%s_%s_%s is either invalid or not present in the specified build." %( chrom, start, end )
|
||||
break
|
||||
if includes_strand_col and strand == "-":
|
||||
sequence = reverse_complement( sequence )
|
||||
|
||||
if output_format == "fasta" :
|
||||
l = len( sequence )
|
||||
c = 0
|
||||
fields = [dbkey, str( chrom ), str( start ), str( end ), strand]
|
||||
meta_data = "_".join( fields )
|
||||
print >> fout, ">%s" %meta_data
|
||||
while c < l:
|
||||
b = min( c + 50, l )
|
||||
print >> fout, sequence[c:b]
|
||||
c = b
|
||||
else: # output_format == "interval"
|
||||
meta_data = "\t".join( fields )
|
||||
print >> fout, meta_data, "\t", sequence
|
||||
except:
|
||||
warning = "Chrom: '%s', start: '%s', end: '%s' is either invalid or not present in build '%s'." %( chrom, start, end, dbkey )
|
||||
warnings.append( warning )
|
||||
skipped_lines += 1
|
||||
if not invalid_line:
|
||||
first_invalid_line = i + 1
|
||||
invalid_line = line
|
||||
continue
|
||||
|
||||
if includes_strand_col:
|
||||
strand = fields[strand_col]
|
||||
if strand not in ['+', '-']:
|
||||
strand = "+"
|
||||
|
||||
sequence = ''
|
||||
if os.path.exists( "%s/%s.nib" % ( nib_path, chrom) ):
|
||||
if chrom in nibs:
|
||||
nib = nibs[chrom]
|
||||
else:
|
||||
nibs[chrom] = nib = bx.seq.nib.NibFile( file( "%s/%s.nib" % ( nib_path, chrom ) ) )
|
||||
try:
|
||||
sequence = nib.get( start, end-start )
|
||||
except:
|
||||
warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'." %( start, end-start, dbkey )
|
||||
warnings.append( warning )
|
||||
skipped_lines += 1
|
||||
if not invalid_line:
|
||||
first_invalid_line = i + 1
|
||||
invalid_line = line
|
||||
continue
|
||||
elif os.path.exists( twobit_path ):
|
||||
if chrom in twobits:
|
||||
t = twobits[chrom]
|
||||
else:
|
||||
twobits[chrom] = t = bx.seq.twobit.TwoBitFile( file( twobit_path ) )
|
||||
try:
|
||||
sequence = t[chrom][start:end]
|
||||
except:
|
||||
warning = "Unable to fetch the sequence from '%d' to '%d' for build '%s'." %( start, end-start, dbkey )
|
||||
warnings.append( warning )
|
||||
skipped_lines += 1
|
||||
if not invalid_line:
|
||||
first_invalid_line = i + 1
|
||||
invalid_line = line
|
||||
continue
|
||||
else:
|
||||
warning = "Chrom '%s' was not found for build '%s'." % ( chrom, dbkey )
|
||||
warnings.append( warning )
|
||||
skipped_lines += 1
|
||||
if not invalid_line:
|
||||
first_invalid_line = i + 1
|
||||
invalid_line = line
|
||||
continue
|
||||
if not sequence:
|
||||
warning = "Chrom: '%s', start: '%s', end: '%s' is either invalid or not present in build '%s'." %( chrom, start, end, dbkey )
|
||||
warnings.append( warning )
|
||||
skipped_lines += 1
|
||||
if not invalid_line:
|
||||
first_invalid_line = i + 1
|
||||
invalid_line = line
|
||||
continue
|
||||
if includes_strand_col and strand == "-":
|
||||
sequence = reverse_complement( sequence )
|
||||
|
||||
if output_format == "fasta" :
|
||||
l = len( sequence )
|
||||
c = 0
|
||||
fields = [dbkey, str( chrom ), str( start ), str( end ), strand]
|
||||
meta_data = "_".join( fields )
|
||||
fout.write( ">%s\n" % meta_data )
|
||||
while c < l:
|
||||
b = min( c + 50, l )
|
||||
fout.write( "%s\n" % str( sequence[c:b] ) )
|
||||
c = b
|
||||
else: # output_format == "interval"
|
||||
meta_data = "\t".join( fields )
|
||||
fout.write( "%s\t%s\n" % ( meta_data, str( sequence ) ) )
|
||||
|
||||
fout.close()
|
||||
if err_msg:
|
||||
stop_err( err_msg )
|
||||
|
||||
if warnings:
|
||||
warn_msg = "Total of %d warnings, 1st is: " % len( warnings )
|
||||
warn_msg += warnings[0]
|
||||
print warn_msg
|
||||
if skipped_lines:
|
||||
print 'Data issue: skipped %d invalid lines starting at line #%d, "%s"' % ( skipped_lines, first_invalid_line, invalid_line )
|
||||
print 'Skipped %d invalid lines starting at line #%d, "%s"' % ( skipped_lines, first_invalid_line, invalid_line )
|
||||
|
||||
if __name__ == "__main__": __main__()
|
||||
@@ -5,7 +5,7 @@
|
||||
<param format="interval" name="input" type="data" label="Fetch sequences corresponding to Query">
|
||||
<validator type="dataset_metadata_in_file" filename="/depot/data2/galaxy/alignseq.loc" metadata_name="dbkey" metadata_column="1" message="Unspecified build (click the pencil icon in the history item) or sequences are not available for the build specified (you may request them via email)." split=" " line_startswith="seq" />
|
||||
</param>
|
||||
<param name="out_format" type="select" label="Output Type">
|
||||
<param name="out_format" type="select" label="Output data type">
|
||||
<option value="fasta">FASTA</option>
|
||||
<option value="interval">Interval</option>
|
||||
</param>
|
||||
@@ -35,11 +35,18 @@
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
Make sure that the genome build is specified for the interval dataset you are extracting sequences for (click the pencil icon if it is not specified).
|
||||
Make sure that the genome build is specified for the dataset from which you are extracting sequences (click the pencil icon in the history item if it is not specified).
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
All of the following will cause a line from the input dataset to be skipped and a warning generated. The number of warnings and skipped lines is documented in the resulting history item.
|
||||
- Any lines that do not contain at least 3 columns, a chromosome and numerical start and end coordinates.
|
||||
- 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.
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**Extract genomic DNA using coordinates from ASSEMBLED genomes and UNassembled genomes** - previously were achieved by two seperate tools.
|
||||
**Extract genomic DNA using coordinates from ASSEMBLED genomes and UNassembled genomes** previously were achieved by two seperate tools.
|
||||
|
||||
-----
|
||||
|
||||
@@ -53,13 +60,13 @@ If strand is not defined, the default value is "+".
|
||||
|
||||
**Example**
|
||||
|
||||
Input dataset::
|
||||
If the input dataset is::
|
||||
|
||||
chr7 127475281 127475310 NM_000230 0 +
|
||||
chr7 127485994 127486166 NM_000230 0 +
|
||||
chr7 127486011 127486166 D49487 0 +
|
||||
|
||||
Fetch genomic DNAs of the above data::
|
||||
Extracting sequences with **FASTA** output data type returns::
|
||||
|
||||
>hg17_chr7_127475281_127475310_+
|
||||
GTAGGAATCGCAGCGCCAGCGGTTGCAAG
|
||||
@@ -74,11 +81,11 @@ Fetch genomic DNAs of the above data::
|
||||
CACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCAC
|
||||
ACACG
|
||||
|
||||
Fetch genomic DNAs of the above data and return as interval format::
|
||||
Extrracting sequences with **Interval** output data type returns::
|
||||
|
||||
chr7 127475281 127475310 NM_000230 0 + GTAGGAATCGCAGCGCCAGCGGTTGCAAG
|
||||
chr7 127485994 127486166 NM_000230 0 + GCCCAAGAAGCCCATCCTGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG
|
||||
chr7 127486011 127486166 D49487 0 + TGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG
|
||||
chr7 127475281 127475310 NM_000230 0 + GTAGGAATCGCAGCGCCAGCGGTTGCAAG
|
||||
chr7 127485994 127486166 NM_000230 0 + GCCCAAGAAGCCCATCCTGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG
|
||||
chr7 127486011 127486166 D49487 0 + TGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
|
||||
Reference in New Issue
Block a user