diff --git a/tools/metag_tools/short_reads_figure_high_quality_length.py b/tools/metag_tools/short_reads_figure_high_quality_length.py
new file mode 100644
index 00000000000..2d6644399fd
--- /dev/null
+++ b/tools/metag_tools/short_reads_figure_high_quality_length.py
@@ -0,0 +1,188 @@
+#! /usr/bin/python
+
+import os, sys, math, tempfile, zipfile, re
+from rpy import *
+
+def stop_err(msg):
+
+ sys.stderr.write(msg)
+ sys.stderr.write("\n")
+ sys.exit()
+
+
+def unzip(zip_file):
+
+ zip_inst = zipfile.ZipFile(zip_file, 'r')
+
+ tmpfilename = tempfile.NamedTemporaryFile().name
+ for name in zip_inst.namelist():
+ file(tmpfilename,'a').write(zip_inst.read(name))
+
+ zip_inst.close()
+
+ return tmpfilename
+
+
+
+def __main__():
+
+ # I/O
+ infile_score_name = sys.argv[1].strip()
+ outfile_R_name = sys.argv[2].strip()
+
+ try:
+ score_threshold = int(sys.argv[3].strip())
+ except:
+ stop_err('Please enter a threshold for quality score.')
+
+ # unzip infile
+ if (zipfile.is_zipfile(infile_score_name)):
+ unzip_infile = unzip(infile_score_name)
+ else: unzip_infile = infile_score_name
+
+ # detect whether it's tabular or fasta format
+ seq_method = None
+ score_file = unzip_infile
+ test_fh = open(score_file,'r')
+ while seq_method is None:
+ read_scorefile = test_fh.readline()
+ if read_scorefile.startswith('#'):
+ continue
+ if not read_scorefile:
+ continue
+ elif read_scorefile.startswith(">"):
+ read_next_line = test_fh.readline()
+ fields = read_next_line.split()
+ for score in fields:
+ try:
+ x = int(score)
+ seq_method = '454'
+ except:
+ seq_method = 'Failed'
+ break
+ elif len(read_scorefile.split('\t')) > 0:
+ fields = read_scorefile.split()
+ for score in fields:
+ try:
+ x = int(score)
+ seq_method = 'solexa'
+ except:
+ seq_method = 'Failed'
+ break
+ else:
+ stop_err('Your input file format does not fit the requirement. Please use either fasta format or tabular format')
+
+ if seq_method == 'Failed':
+ stop_err('Unable to determine the file format. Please use either fasta-like format (except title lines, the file contains only numeric values) or tabular format')
+
+ test_fh.close()
+
+
+ # R
+ cont_high_quality = []
+
+ invalid_lines = 0
+ invalid_scores = 0
+
+ if (seq_method == 'solexa'):
+ for i, line in enumerate(open(score_file)):
+ line = line.rstrip('\r\n')
+ if line.startswith('#'):
+ continue
+ if not line:
+ continue
+
+ try:
+ each_loc = line.split('\t')
+ except:
+ invalid_lines += 1
+ continue
+
+ for j, each_base in enumerate(each_loc):
+ each_nuc_error = each_base.split()
+
+ try:
+ each_nuc_error[0] = int(each_nuc_error[0])
+ each_nuc_error[1] = int(each_nuc_error[1])
+ each_nuc_error[2] = int(each_nuc_error[2])
+ each_nuc_error[3] = int(each_nuc_error[3])
+ big = max(each_nuc_error)
+ except:
+ invalid_scores += 1
+ big = 0
+
+ if j == 0:
+ cont_high_quality.append(1)
+ else:
+ if big >= score_threshold:
+ cont_high_quality[len(cont_high_quality)-1] += 1
+ else:
+ cont_high_quality.append(1)
+ else:
+ tmp_score = ''
+ for i, line in enumerate(open(score_file)):
+ if line.startswith('#'):
+ continue
+ if not line:
+ continue
+ if line.startswith('>'):
+ if len(tmp_score) > 0:
+ each_loc = tmp_score.split()
+ for j, each_base in enumerate(each_loc):
+ try:
+ each_base = int(each_base)
+ except:
+ invalid_scores += 1
+ each_base = 0
+ if j == 0:
+ cont_high_quality.append(1)
+ else:
+ if each_base >= score_threshold:
+ cont_high_quality[len(cont_high_quality)-1] += 1
+ else:
+ cont_high_quality.append(1)
+ tmp_score = ''
+ else:
+ tmp_score = tmp_score + ' ' + line
+
+ if len(tmp_score) > 0:
+ each_loc = tmp_score.split()
+ for j, each_base in enumerate(each_loc):
+ try:
+ each_base = int(each_base)
+ except:
+ invalid_scores += 1
+ each_base = 0
+ if j == 0:
+ cont_high_quality.append(1)
+ else:
+ if each_base >= score_threshold:
+ cont_high_quality[len(cont_high_quality)-1] += 1
+ else:
+ cont_high_quality.append(1)
+
+ # throw messages of invalid values
+ if invalid_lines > 0:
+ print >> sys.stdout, 'Skipped %d lines due to invalid format' %(invalid_lines)
+ if invalid_scores > 0:
+ print >> sys.stdout, 'Skipped %d scores due to invalid values' %(invalid_scores)
+
+ # generate pdf figures
+ cont_high_quality = array (cont_high_quality)
+
+ outfile_R_pdf = outfile_R_name
+ r.pdf(outfile_R_pdf)
+
+ title = "histogram of continueous high quality scores"
+ xlim_range = [1,max(cont_high_quality)]
+ nclass = max(cont_high_quality)
+ if nclass > 100: nclass = 100
+ r.hist(cont_high_quality,probability=True, xlab="Continueous High Quality Score length (bp)", ylab="Frequency (%)", xlim=xlim_range, main=title, nclass=nclass)
+
+ if zipfile.is_zipfile(infile_score_name) and os.path.exists(unzip_infile):
+ os.remove(unzip_infile)
+
+ r.dev_off()
+ r.quit(save = "no")
+
+if __name__=="__main__":__main__()
diff --git a/tools/metag_tools/short_reads_figure_high_quality_length.xml b/tools/metag_tools/short_reads_figure_high_quality_length.xml
new file mode 100644
index 00000000000..7581ce659f7
--- /dev/null
+++ b/tools/metag_tools/short_reads_figure_high_quality_length.xml
@@ -0,0 +1,59 @@
+
+ of high quality score reads
+
+short_reads_figure_high_quality_length.py $input1 $output1 $input2
+
+
+
+
+
+
+
+
+
+
+
+
+
+.. class:: warningmark
+
+**TIP**. To use this tool your 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 Quality Files generated by Roche (454), Illumina (Solexa), or ABI SOLiD machines and builds a histogram of lengths of high quality reads.
+
+-----
+
+**Examples of Quality Data**
+
+Roche (454) or ABI SOLiD data::
+
+ >seq1
+ 23 33 34 25 28 28 28 32 23 34 27 4 28 28 31 21 28
+
+Illumina (Solexa) data::
+
+ -40 -40 40 -40 -40 -40 -40 40
+
+-----
+
+**Note**
+
+- Quality score data::
+
+ >seq1
+ 23 33 34 25 28 28 28 32 23 34 27 4 28 28 31 21 28
+
+- If the threshold was set to 20::
+
+ lengths of continuous quality scores that are higher than 20 is 11, 5 (a low quality score 4 in the middle)
+
+ The histogram will be built based on the numbers (11, 5).
+
+- For Illumina (Solexa) data, only the maximal of 4 values will be considered.
+
+
+
diff --git a/tools/metag_tools/short_reads_figure_length.py b/tools/metag_tools/short_reads_figure_length.py
deleted file mode 100644
index dba6a07d493..00000000000
--- a/tools/metag_tools/short_reads_figure_length.py
+++ /dev/null
@@ -1,206 +0,0 @@
-#! /usr/bin/python
-"""
-Galaxy
-Input:
- sequence file(s): zip or text file --> for 454 and Solexa
-Output:
- an array of lengths of the reads
-----
-Wen-Yu Chung
-"""
-import os, sys, math, tempfile, zipfile, re
-from rpy import *
-
-# default and initialize values
-number_of_points = 20
-read_seqfile = []
-read_scorefile = []
-title_keys = []
-seq_hash = {}
-score_hash = {}
-trim_seq_hash = {}
-tmp_seq = ''
-tmp_score = [] # change variables here
-
-length_before_trim = []
-length_after_trim = []
-score_points = []
-
-database_tmp = "/tmp/" # default dir: current directory
-if (not os.path.isdir(database_tmp)):
- os.mkdir(database_tmp)
-
-# functions
-def stop_err(msg):
-
- sys.stderr.write(msg)
- sys.stderr.write("\n")
- sys.exit()
-
-
-def read_input_files(file_list):
-
- read_file = []
-
- for file_name in file_list:
- infile = open(file_name,'r')
- read_file = infile.readlines()
- infile.close()
-
- return read_file
-
-
-def unzip_files(file_name):
-
- read_file = []
-
- temp_dir_name = tempfile.mkdtemp() + '/' #dir=database_tmp
- temp_new_dest = temp_dir_name + file_name.split('/')[-1]
- command_line = 'cp ' + file_name + ' ' + temp_dir_name + '\nunzip ' + temp_new_dest + ' -d ' + temp_dir_name
-
- os.system(command_line)
- os.remove(temp_new_dest)
-
- dest = os.listdir(temp_dir_name)
- if ((len(dest) == 1) and (os.path.isdir(dest[0]))):
- new_file_dir = temp_dir_name + dest[0] + '/'
- dest = os.listdir(new_file_dir)
- else:
- new_file_dir = temp_dir_name
-
- for filex in dest:
- filex = new_file_dir + filex
- read_file.append(filex)
-
- return new_file_dir, read_file
-
-
-def parse_solexa_files(result, type):
- # parse solexa seq and probability files only
- # not fasta format
- # assign numeric id
-
- tmp_hash = {}
- tmp_title = ''
- tmp_seq = ''
- tmp_id = 0
-
- if (type == 'seq'):
- for each_line in result:
- tmp_seq = ''
- each_line.strip('\r\n')
- tmp_id += 1
- fasta_id = '>' + str(tmp_id)
- (run_id, file_id, read_x, read_y, read) = each_line.split()
- tmp_hash[(fasta_id)] = read
- else: # type == score
- for each_line in result:
- tmp_seq = []
- each_line.strip('\r\n')
- tmp_id += 1
- fasta_id = '>' + str(tmp_id)
- each_loc = each_line.split('\t')
- for each_base in each_loc:
- each_nuc_error = each_base.split()
- each_nuc_error[0] = int(each_nuc_error[0])
- each_nuc_error[1] = int(each_nuc_error[1])
- each_nuc_error[2] = int(each_nuc_error[2])
- each_nuc_error[3] = int(each_nuc_error[3])
- big = max(each_nuc_error)
- #baseIndex = each_nuc_error.index(big)
- tmp_seq.append(big)
- tmp_hash[(fasta_id)] = tmp_seq
-
- return tmp_hash
-
-
-def parse_fasta_format(result):
- # detect whether it's score or seq files
- # return a hash: key = title and value = seq
-
- tmp_hash = {}
- tmp_title = ''
- tmp_seq = ''
-
- for each_line in result:
- each_line = each_line.strip('\r\n')
- if (each_line[0] == '>'):
- if (len(tmp_seq) > 0):
- tmp_hash[(tmp_title)] = tmp_seq
- tmp_title = each_line
- tmp_seq = ''
- else:
- tmp_seq = tmp_seq + each_line
- if (each_line.split()[0].isdigit()):
- tmp_seq = tmp_seq + ' '
- if (len(tmp_seq) > 0):
- tmp_hash[(tmp_title)] = tmp_seq
-
- return tmp_hash
-
-
-def generate_hist_figure():
- # R module and code
-
- tmp_array = []
-
- # data for R hist
- # write the length only!
- outfile = open(outfile_R_name,'w')
- title_keys = seq_hash.keys()
- i = 0
- for read_title in title_keys:
- i += 1
- tmp_seq = seq_hash[(read_title)]
- length_before_trim.append(len(tmp_seq))
- print >> outfile,"%d\t%d" %(i, len(tmp_seq))
- outfile.close()
-
- max_length_before_trim = max(length_before_trim)
-
- # generate pdf figures
- """
- outfile_R_pdf = outfile_R_name
- r.pdf(outfile_R_pdf)
-
- title = infile_seq_name.split('/')[-1]
- xlim_range = [1,max_length_before_trim]
- r.hist(length_before_trim,prob=1, xlab="Read length (bp)", ylab="Frequency (%)", xlim=xlim_range, main=title, nclass=100)
- r.dev_off()
- """
- return 0
-
-
-# I/O
-infile_seq_name = sys.argv[1].strip()
-outfile_R_name = sys.argv[2].strip()
-
-# to unzip or not unzip file
-tmp_seq_dir = ''
-seq_file_list = []
-if (zipfile.is_zipfile(infile_seq_name)): (tmp_seq_dir, seq_file_list) = unzip_files(infile_seq_name)
-else: seq_file_list = [infile_seq_name]
-read_seqfile = read_input_files(seq_file_list)
-if (os.path.isdir(tmp_seq_dir)):
- for file_name in seq_file_list:
- os.remove(file_name)
- os.removedirs(tmp_seq_dir)
-
-# detect whether it's tabular or fasta format
-if read_seqfile[0].startswith(">"):
- seq_method = '454'
-else:
- seq_method = 'solexa'
-
-
-# parse files
-if (seq_method == 'solexa'):
- seq_hash = parse_solexa_files(read_seqfile, 'seq')
-else: # the other two are both fasta format
- seq_hash = parse_fasta_format(read_seqfile)
-
-# R
-status = generate_hist_figure()
-
-# close and remove temporary files
-r.quit(save = "no")
diff --git a/tools/metag_tools/short_reads_figure_length.xml b/tools/metag_tools/short_reads_figure_length.xml
deleted file mode 100644
index eb4cbe16721..00000000000
--- a/tools/metag_tools/short_reads_figure_length.xml
+++ /dev/null
@@ -1,68 +0,0 @@
-
-on sequence file
-
-short_reads_figure_length.py $input1 $output1
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-**What it does**
-
-This tool outputs an array of lengths of reads.
-
-Currently, the tool accepts two file formats:
-
-- fasta: for **454** data
-- tabular: for **Solexa** data. Only take the maximal value from each 4 quality scores.
-
-You can also upload your files in a zip file format and use this tool.
-To see the distribution, please use the Histogram tool in Graph/Display data.
-
------
-
-**Example**
-
-454 data::
-
- >EYKX4VC01B65GS length=54 xy=0784_1754 region=1 run=R_2007_11_07_16_15_57_
- CCGGTATCCGGGTGCCGTGATGAGCGCCACCGGAACGAATTCGACTATGCCGAA
- >EYKX4VC01BNCSP length=187 xy=0558_3831 region=1 run=R_2007_11_07_16_15_57_
- CTTACCGGTCACCACCGTGCCTTCAGGATTGATCGCCAGATCGGTCGGTGCGTCAGGCGG
- GGTGACATCGCCCACCACGGTACTCACTGGCTGGCTCTGGTTCCCGGCGGCATCGGAGGC
- CACCACGTTGAGGGTATTCCCCTCGGTTTGTGGCTCGGTGAGAACCACGTTGTAGTCGCC
- ATTGGTC
-
-Solexa data::
-
- 5 300 902 419 GACTCATGATTTCTTACCTATTAGTGGTTGAACATC
- 5 300 880 431 GTGATATGTATGTTGACGGCCATAAGGCTGCTTCTT
- 5 300 896 461 GTTGTCGATAGAACTTCATGTGCCTGTAAAACAAGT
- 5 300 890 751 ACCAACCAGAACGTGAAAAAGCGTCCTGCGTGTAGC
- 5 300 897 443 GTTTATGTTGGTTTCATGGTTTTGTCTAACTTTATC
- 5 300 906 879 GCTTTACCGTCTTTCCAGAAATTGTTCCAAGTATCG
-
-An array of lengths of the reads::
-
- Please see Graph/Display Data --> Histogram.
-
-