opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

update short read tools.

Browse Source
This commit is contained in:
Wen-Yu Chung
2008-03-25 22:07:46 +00:00
parent a11018e0d0
commit c8a7bae1a2
8 changed files with 236 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/metag_tools/blat_coverage_report.xml
+14 -9
View File
@@ -17,7 +17,7 @@
.. class:: warningmark
Only Accept **BLAT output type pslx**.
**Note**. Only works for BLAT output type **pslx** (hint: add -out=pslx in the command).
-----
@@ -31,23 +31,22 @@
- 3rd column: the nucleotide from reference genome at the chromosome location (2nd column).
- 4th column: total coverage of the reads (number of reads that mapped to the chromosome location).
- 4th column: total coverage of the reads (number of reads that were mapped to the chromosome location).
- 5th column: percentage of reads that support a nucleotide **A** at this location.
- 5th column: percentage of reads that support nucleotide **A** at this location.
- 6th column: percentage of reads that support a nucleotide **T** at this location.
- 6th column: percentage of reads that support nucleotide **T** at this location.
- 7th column: percentage of reads that support a nucleotide **C** at this location.
- 7th column: percentage of reads that support nucleotide **C** at this location.
- 8th column: percentage of reads that support a nucleotide **G** at this location.
If there is no read to support a particular nucleotide, the column is left empty.
- 8th column: percentage of reads that support nucleotide **G** at this location.
-----
**Example**
- The BLAT results look like the following::
- The BLAT pslx results look like the following (tab separated with sequence at the end)::
30 0 0 0 0 0 0 0 + seq0 30 0 30 chr 4639675 4549207 4549237 1 30, 0, 4549207, cggacagcgccgccaccaacaaagccacca, cggacagcgccgccaccaacaaagccacca,
30 0 0 0 0 0 0 0 + seq1 30 0 30 chr 4639675 614777 614807 1 30, 0, 614777, aaaacaccggatgctccggcgctggcagat, aaaacaccggatgctccggcgctggcagat,
@@ -70,5 +69,11 @@
Only show part of the result.
-----
**Reference**
BLAT: Kent, W James, BLAT--the BLAST-like alignment tool. (2002) Genome Research:12(4) 656-664.
</help>
</tool>
tools/metag_tools/blat_wrapper.py
+27 -7
View File
@@ -2,6 +2,12 @@
import os, sys, tempfile
def stop_err(msg):
sys.stderr.write(msg)
sys.stderr.write("\n")
sys.exit()
def check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR ):
nib_file = "%s/alignseq.loc" % GALAXY_DATA_INDEX_DIR
nib_path = ''
@@ -39,13 +45,28 @@ def __main__():
target_file = sys.argv[2]
query_file = sys.argv[3]
output_file = sys.argv[4]
min_iden = sys.argv[5]
tile_size = sys.argv[6]
one_off = sys.argv[7]
try:
min_iden = float(sys.argv[5])
except:
stop_err('Invalid value for minimal identity')
try:
tile_size = int(sys.argv[6])
assert tile_size >= 6 and tile_size <= 18
except:
stop_err('Invalid value for tile size. DNA word size must be between 6 and 18.')
try:
one_off = int(sys.argv[7])
except:
stop_err('Invalid value for mismatch numbers in the word')
GALAXY_DATA_INDEX_DIR = sys.argv[8]
all_files = []
if (source_format == '0'):
# check target genome
dbkey = target_file
if dbkey == '?':
@@ -54,8 +75,7 @@ def __main__():
nib_path = check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR )
twobit_path = check_twobit_file( dbkey, GALAXY_DATA_INDEX_DIR )
if not os.path.exists( nib_path ) and not os.path.exists( twobit_path ):
print >> sys.stdout, "No sequences are available for %s, request them by reporting this error." % dbkey
sys.exit()
stop_err("No sequences are available for %s, request them by reporting this error." % dbkey)
# check the query file, see whether all of them are legitimate sequence
if (nib_path and os.path.isdir(nib_path)):
@@ -65,8 +85,8 @@ def __main__():
compress_files = [twobit_path]
target_path = ""
else:
print >> sys.stdout, "Requested genome build has no available sequence."
sys.exit()
stop_err("Requested genome build has no available sequence.")
for file in compress_files:
file = target_path + '/' + file
file = os.path.normpath(file)
tools/metag_tools/blat_wrapper.xml
+11 -7
View File
@@ -21,7 +21,7 @@
</conditional>
<param name="input_query" type="data" format="fasta" label="Query Sequence"/>
<param name="iden" type="float" size="15" value="90.0" label="Minimal Identity (-minIdentity)" />
<param name="tile_size" type="integer" size="15" value="11" label="Minimal Size of Match (-tileSize)" help="for shorter reads, use size = 8"/>
<param name="tile_size" type="integer" size="15" value="11" label="Minimal Size of Exact Match (-tileSize)" help="for shorter reads, use size = 8. Must be between 6 and 18."/>
<param name="one_off" type="integer" size="15" value="0" label="Number of Mismatch in the Word (-oneOff)" />
</inputs>
<outputs>
@@ -40,29 +40,33 @@
</tests>
<help>
.. class:: warningmark
Use smaller word size (*Minimal Size of Exact Match*) will increase the computational time.
-----
**What it does**
This tool runs sequence alignment program on your short reads dataset against a genome build.
This tool runs alignment program **BLAT**. Your short reads file is searched against a genome build (select from table) or another uploaded file.
-----
**Example**
Input as multiple fasta file::
- Input as multiple fasta file::
&gt;seq1
TGGTAATGGTGGTTTTTTTTTTTTTTTTTTATTTTT
Search against ce2 (C. elegans March 2004)
(This is output type = pslx)::
- Search against ce2 (C. elegans March 2004)::
25 1 0 0 0 0 0 0 + seq1 36 10 36 chrI 15080483 9704438 9704464 1 26, 10, 9704438, ggttttttttttttttttttattttt, ggtttttttttttttttttttttttt,
27 0 0 0 0 0 1 32 + seq1 36 9 36 chrI 15080483 1302536 1302595 2 21,6, 9,30, 1302536,1302589, tggtttttttttttttttttt,attttt, tggtttttttttttttttttt,attttt,
-----
.. class:: infomark
**Reference**
BLAT: Kent, W James, BLAT--the BLAST-like alignment tool. (2002) Genome Research:12(4) 656-664.
tools/metag_tools/megablast_wrapper.py
+28 -9
View File
@@ -6,20 +6,42 @@ run megablast for metagenomics data
import sys, os, tempfile, subprocess
#from megablast_xml_parser import *
def stop_err(msg):
sys.stderr.write(msg)
sys.stderr.write("\n")
sys.exit()
def __main__():
# file I/O
db_build = sys.argv[1]
query_filename = sys.argv[2]
output_filename = sys.argv[3]
query_filename = sys.argv[2].strip()
output_filename = sys.argv[3].strip()
# megablast parameters
mega_word_size = sys.argv[4] # -W
mega_iden_cutoff = sys.argv[5] # -p
mega_disc_word = sys.argv[6] # -t
try:
mega_word_size = int(sys.argv[4]) # -W
except:
stop_err('Invalid value for word size')
try:
mega_iden_cutoff = float(sys.argv[5]) # -p
except:
stop_err('Invalid value for identity cut-off')
try:
mega_disc_word = sys.argv[6] # -t
except:
stop_err('Invalid value for discontiguous word template')
mega_disc_type = sys.argv[7] # -N
mega_filter = sys.argv[8] # -F
GALAXY_DATA_INDEX_DIR = sys.argv[9]
DB_LOC = "%s/blastdb.loc" % GALAXY_DATA_INDEX_DIR
output_file = open(output_filename, 'w')
# prepare the database
@@ -35,8 +57,7 @@ def __main__():
# prepare to run megablast
retcode = subprocess.call('which megablast 2>&1', shell='True')
if retcode < 0:
print >> sys.stderr, "Cannot locate megablast."
sys.exit()
stop_err("Cannot locate megablast.")
for chunk in db[(db_build)]:
megablast_arguments = ["megablast", "-d", chunk, "-i", query_filename]
@@ -44,9 +65,7 @@ def __main__():
megablast_user_inputs = ["-W", mega_word_size, "-p", mega_iden_cutoff, "-t", mega_disc_word, "-N", mega_disc_type, "-F", mega_filter]
megablast_command = " ".join(megablast_arguments) + " " + " ".join(megablast_parameters) + " " + " ".join(megablast_user_inputs) + " 2>&1"
# use Anton's parser
megablast_output = os.popen(megablast_command)
#parse_megablast_xml_output(megablast_output,output_file)
# to avoid reading whole file into memory
for i, line in enumerate(megablast_output):
line = line.rstrip('\r\n')
tools/metag_tools/megablast_xml_parser.py
+53 -32
View File
@@ -1,10 +1,13 @@
import cElementTree #python 2.5 xml.etree.cElementTree
#! /usr/bin/python
import cElementTree
import sys, os
def parse_megablast_xml_output(infile,outfile = sys.stdout):
source = infile #sys.argv[1]
def parse_megablast_xml_output(infile_name,outfile_name):
source = infile_name
outfile = open(outfile_name, 'w')
hspTags = [
"Hsp_bit-score",
"Hsp_evalue",
@@ -24,41 +27,59 @@ def parse_megablast_xml_output(infile,outfile = sys.stdout):
hspData = []
# get an iterable
context = cElementTree.iterparse(source, events=("start", "end"))
try:
context = cElementTree.iterparse(source, events=("start", "end"))
except:
print >> sys.stderr, "Please note: this tool is for megablast output option -m 7 only."
print >> sys.stderr, "The file has inappropriate format for this tool."
sys.exit()
# turn it into an iterator
context = iter(context)
# get the root element
event, root = context.next()
try:
event, root = context.next()
except:
print >> sys.stderr, "Please note: this tool is for megablast output option -m 7 only."
print >> sys.stderr, "The file has inappropriate format for this tool."
sys.exit()
try:
for event, elem in context:
# for every <Iteration> tag
if event == "end" and elem.tag == "Iteration":
query = elem.findtext("Iteration_query-def")
qLen = elem.findtext("Iteration_query-len")
# for every <Hit> within <Iteration>
for hit in elem.findall("Iteration_hits/Hit/"):
subject = hit.findtext("Hit_id")
sLen = hit.findtext("Hit_len")
# for every <Hsp> within <Hit>
for hsp in hit.findall("Hit_hsps/Hsp"):
for tag in hspTags:
hspData.append(hsp.findtext(tag))
print >> outfile, query, '\t', qLen, '\t', subject, '\t', sLen, '\t', hspData
hspData = []
for event, elem in context:
# for every <Iteration> tag
if event == "end" and elem.tag == "Iteration":
query = elem.findtext("Iteration_query-def")
qLen = elem.findtext("Iteration_query-len")
# for every <Hit> within <Iteration>
for hit in elem.findall("Iteration_hits/Hit/"):
subject = hit.findtext("Hit_id")
sLen = hit.findtext("Hit_len")
# for every <Hsp> within <Hit>
for hsp in hit.findall("Hit_hsps/Hsp"):
for tag in hspTags:
hspData.append(hsp.findtext(tag))
print >> outfile, query, '\t', qLen, '\t', subject, '\t', sLen, '\t', hspData
hspData = []
# prevents ElementTree from growing large datastructure
root.clear()
elem.clear()
# prevents ElementTree from growing large datastructure
root.clear()
elem.clear()
except:
print >> sys.stderr, "Your file may contain tags that are not recognized by the parser."
print >> sys.stderr, "Please use megablast output option -m 7 only."
sys.exit()
outfile.close()
return
def __main__():
megablast_command = "megablast -d Ecoli.fa -i Ecoli.30bp.random.fa -m 7 "
infile = os.popen(megablast_command)
outfile = None
if len(sys.argv) > 1:
outfile = open(sys.argv[1],'w')
parse_megablast_xml_output(infile, outfile)
infile_name = sys.argv[1]
outfile_name = sys.argv[2]
parse_megablast_xml_output(infile_name, outfile_name)
if __name__ == "__main__": __main__()
tools/metag_tools/megablast_xml_parser.xml
+54
View File
@@ -0,0 +1,54 @@
<tool id="megablast_xml_parser" name="Parse megablast xml output">
<description> </description>
<command interpreter="python">megablast_xml_parser.py $input1 $output1</command>
<inputs>
<param name="input1" type="data" format="txt" label="Megablast XML Output" />
</inputs>
<outputs>
<data name="output1" format="tabular"/>
</outputs>
<tests>
<test>
<param name="input1" value="megablast_xml_parser_test1.txt" />
<output name="output1" file="megablast_xml_parser_test1.out" ftype="tabular" />
</test>
</tests>
<help>
.. class:: warningmark
**TIP**. To get xml output from megablast, add option **-m 7** in the command line.
.. class:: warningmark
**Important**. Please **zip** your xml file and upload the zipped file. This will help speeding up the process. (hint: under shell command, try: gzip -c your_xml_file > your_xml_file.gz)
-----
**What it does**
This tool is for users to upload their own megablast xml output and converts to tabular format, showing all information for retrieving alignment blocks.
-----
**Example**
- Query sequence::
&gt;seq1
CGGACAGCGCCGCCACCAACAAAGCCACCA
- Parsed output::
seq1 30 gnl|BL_ORD_ID|0 5528445 ['59.96', '8.38112e-11', '1', '30', '5436010', '5436039', '1', '1', '30', '30', 'CGGACAGCGCCGCCACCAACAAAGCCACCA', 'CGGACAGCGCCGCCACCAACAAAGCCACCA', '||||||||||||||||||||||||||||||']
-----
**Reference**
**megablast**: Zhang et al. A Greedy Algorithm for Aligning DNA Sequences. 2000. JCB: 203-214.
</help>
</tool>
tools/metag_tools/short_reads_trim_seq.py
+16 -4
View File
@@ -28,11 +28,15 @@ def unzip(zip_file):
def show_output(outfile_seq, seq_title, segments):
if (len(segments) > 1):
for i in range(len(segments)):
print >> outfile_seq, "%s_%d\n%s" % (seq_title, i, segments[i])
elif (len(segments[0]) > 0):
print >> outfile_seq, "%s\n%s" % (seq_title, segments[0])
else:
return
return
@@ -105,6 +109,7 @@ def __main__():
infile_seq_name = sys.argv[5].strip()
infile_score_name = sys.argv[6].strip()
special_argument = sys.argv[7].strip()
if seq_method == '454':
keep_homopolymers = special_argument
else:
@@ -119,10 +124,8 @@ def __main__():
unzip_score_infile = unzip(infile_score_name)
else: unzip_score_infile = infile_score_name
outfile_seq = open(outfile_seq_name,'w')
if (os.path.exists(unzip_seq_infile) and os.path.exists(unzip_score_infile)):
# read one sequence
to_find_score = True
@@ -156,7 +159,7 @@ def __main__():
try:
each_score = int(each_score)
except:
stop_err('Score file contains non-numerical values at line %d' %(i))
stop_err('Score file contains non-numerical values: %s' %(each_score))
if not score: score = score_line
else:
score = score + ' ' + score_line
@@ -177,7 +180,7 @@ def __main__():
try:
each_score = int(each_score)
except:
stop_err('Score file contains non-numerical values at line %d' %(i))
stop_err('Score file contains non-numerical values: %s' %(each_score))
if not score: score = score_line
else:
score = score + ' ' + score_line
@@ -193,7 +196,16 @@ def __main__():
if line.startswith('#'):
continue
seq_title = '>' + str(i)
# the last column of solexa file is the sequence column
seq = line.split()[-1]
seq.replace('.','N')
seq.replace('-','N')
try:
assert seq.isalpha() is True
except:
stop_err('Invalid characters in the read at line %d' %(i))
score = score_fh.readline()
each_loc = score.split('\t')
score = []
tools/metag_tools/short_reads_trim_seq.xml
+33 -15
View File
@@ -28,7 +28,7 @@
</when>
</conditional>
<param name="trim" type="integer" size="5" value="20" label="Minimal quality score" />
<param name="length" type="integer" size="5" value="0" label="Minimal length of trimmed reads to report" help="To output the read if its length is longer than this threshold. Use 0 to return the longest substring" />
<param name="length" type="integer" size="5" value="0" label="Minimal length of trimmed reads to report" help="To output the read if its length is longer than this threshold. Use 0 to return the longest substring." />
</page>
</inputs>
@@ -61,25 +61,21 @@
.. class:: warningmark
Please select correct sequencing method.
**Note**. Please select sequencing method.
.. class:: warningmark
.. class:: warningmark
Your **Quality Score** file needs to be a specific Galaxy **qualityscore** datatype (unless it is in zip format). Please click on the pencil symbol in the history panel to change the file format.
**TIP**. To use this tool your quality score dataset needs to be in *Quality Score* format. Click pencil icon next to your dataset to set datatype to *Quality Score*.
-----
**What it does**
This tool takes raw data from sequencing facility and outputs a multi-fasta file.
This tool takes Sequencing Files and Quality Files generated by Roche (454), Illumina (Solexa), or ABI SOLiD machines and trims the sequences based on the quality scores.
Accept the following two file formats:
- fasta: for **454** and **SOLiD** data
- tabular: for **Solexa** data. Only take the maximal value from each 4 quality scores.
If *minimal length of trimmed reads to report* is set to a number larger than zero, the tool ouputs any substrings that are longer than this threshold.
If there are more than one substring, the tool outputs all of them in separated fasta format.
- If *minimal length of trimmed reads to report* is set to a number larger than zero, the tool ouputs any substrings that are longer than this threshold.
- If *minimal length of trimmed reads to report* is set to zero, and several substrings are of the same maximal length, the tool outputs all of them in separated fasta format.
-----
@@ -100,7 +96,7 @@ For 454 and SOLiD files
- if use **20** as *minimal quality score*, and
- use **0** as *minimal length of trimmed reads to report*
- the output will return the longest sub-sequence::
- the output will return the longest substring::
&gt;seq1
GGTATC
@@ -124,12 +120,34 @@ For Solexa files
- if use **20** as *minimal quality score*, and
- use **0** as *minimal length of trimmed reads to report*
- the output will return the longest sub-sequence::
- the output will return the longest substring::
&gt;15774
TTCTTACCTATTAGTGGTTGAACA
Note: The fasta title will be replaced by a unique number in Solexa output.
Note: The fasta title will be a unique number automatically generated by the tool.
-----
**Note**. Fasta title changed if there was more than one candidate.
If several substrings are of the same maximal length or longer than the non-zero threshold, all the substrings will be outputted. The new fasta title will be the original title with an underscore symbol and a number representing the order of the substrings.
- The length of the longest substring is 6 and three substrings are of this length::
&gt;seq1_0
GGTATC
&gt;seq1_1
AATATC
&gt;seq1_2
GGTATC
- The *minimal length of trimmed reads to report* is 10 and two substrings are longer than this threshold::
&gt;15774_0
TTCTTACCTAGTGGTATTAGTGGTTGAACA
&gt;15774_1
ACCTATTACCTAGGTGGTTGA
</help>
</tool>