diff --git a/tools/metag_tools/blat_coverage_report.xml b/tools/metag_tools/blat_coverage_report.xml index 215b5b745a9..1e96e1c7355 100644 --- a/tools/metag_tools/blat_coverage_report.xml +++ b/tools/metag_tools/blat_coverage_report.xml @@ -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. + diff --git a/tools/metag_tools/blat_wrapper.py b/tools/metag_tools/blat_wrapper.py index 72e192a4a34..0fb7c30e931 100644 --- a/tools/metag_tools/blat_wrapper.py +++ b/tools/metag_tools/blat_wrapper.py @@ -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) diff --git a/tools/metag_tools/blat_wrapper.xml b/tools/metag_tools/blat_wrapper.xml index a7099014cdf..07db56e7cc3 100644 --- a/tools/metag_tools/blat_wrapper.xml +++ b/tools/metag_tools/blat_wrapper.xml @@ -21,7 +21,7 @@ - + @@ -40,29 +40,33 @@ +.. 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:: >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. diff --git a/tools/metag_tools/megablast_wrapper.py b/tools/metag_tools/megablast_wrapper.py index 4b2c93b57e7..a2dd1eff357 100644 --- a/tools/metag_tools/megablast_wrapper.py +++ b/tools/metag_tools/megablast_wrapper.py @@ -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') diff --git a/tools/metag_tools/megablast_xml_parser.py b/tools/metag_tools/megablast_xml_parser.py index def8ba3aca1..fa63f855d42 100644 --- a/tools/metag_tools/megablast_xml_parser.py +++ b/tools/metag_tools/megablast_xml_parser.py @@ -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 tag + if event == "end" and elem.tag == "Iteration": + query = elem.findtext("Iteration_query-def") + qLen = elem.findtext("Iteration_query-len") + # for every within + for hit in elem.findall("Iteration_hits/Hit/"): + subject = hit.findtext("Hit_id") + sLen = hit.findtext("Hit_len") + # for every within + 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 tag - if event == "end" and elem.tag == "Iteration": - query = elem.findtext("Iteration_query-def") - qLen = elem.findtext("Iteration_query-len") - # for every within - for hit in elem.findall("Iteration_hits/Hit/"): - subject = hit.findtext("Hit_id") - sLen = hit.findtext("Hit_len") - # for every within - 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__() diff --git a/tools/metag_tools/megablast_xml_parser.xml b/tools/metag_tools/megablast_xml_parser.xml new file mode 100644 index 00000000000..3bcaa14e123 --- /dev/null +++ b/tools/metag_tools/megablast_xml_parser.xml @@ -0,0 +1,54 @@ + + +megablast_xml_parser.py $input1 $output1 + + + + + + + + + + + + + + +.. 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:: + + >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. + + + + \ No newline at end of file diff --git a/tools/metag_tools/short_reads_trim_seq.py b/tools/metag_tools/short_reads_trim_seq.py index 8baae26fb0b..26d63a14996 100644 --- a/tools/metag_tools/short_reads_trim_seq.py +++ b/tools/metag_tools/short_reads_trim_seq.py @@ -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 = [] diff --git a/tools/metag_tools/short_reads_trim_seq.xml b/tools/metag_tools/short_reads_trim_seq.xml index 1dcada9effc..fb00d1f6da4 100644 --- a/tools/metag_tools/short_reads_trim_seq.xml +++ b/tools/metag_tools/short_reads_trim_seq.xml @@ -28,7 +28,7 @@ - + @@ -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:: >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:: >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:: + + >seq1_0 + GGTATC + >seq1_1 + AATATC + >seq1_2 + GGTATC + +- The *minimal length of trimmed reads to report* is 10 and two substrings are longer than this threshold:: + + >15774_0 + TTCTTACCTAGTGGTATTAGTGGTTGAACA + >15774_1 + ACCTATTACCTAGGTGGTTGA