diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index 27b370da5a5..266041134c5 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -318,6 +318,7 @@
+
+
+
diff --git a/tools/metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.py b/tools/metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.py
index 4bb62fb8df1..74829f8a7ac 100644
--- a/tools/metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.py
+++ b/tools/metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.py
@@ -1,42 +1,8 @@
#! /usr/bin/python
-"""
-Galaxy
-to commit
-4.2 for i, line in enumerate(file(filename)):
-4.4 stop passing read_seqfile and read_scorefile around
-4.5 use --> if __name__ == "__main__": __main__()
-5. use galaxy extention
-Input:
- quality score file(s): zip or text file
-Output:
-pdf files from R
- score distribution --> boxplot, Solexa as 36bp, 454 as 20 points
-----
-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)
@@ -44,23 +10,11 @@ def stop_err(msg):
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_dir_name = tempfile.mkdtemp() + '/'
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
@@ -81,125 +35,190 @@ def unzip_files(file_name):
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 merge_to_20_datapoints(tmp_score):
+ number_of_points = 20
+ read_length = len(tmp_score)
-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 = ''
+ step = int(math.floor((read_length-1)*1.0/number_of_points))
+ score = []
+ point = 1
+ point_sum = 0
+ step_average = 0
+ score_points_tmp = 0
+
+ for i in xrange(1,read_length):
+ if (i < (point * step)):
+ point_sum += int(tmp_score[i])
+ step_average += 1
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
+ point_avg = point_sum*1.0/step_average
+ score.append(point_avg)
+ point += 1
+ point_sum = 0
+ step_average = 0
+
+ if (step_average > 0):
+ point_avg = point_sum*1.0/step_average
+ score.append(point_avg)
-
-def generate_boxplot_figure():
- # R module and code
-
- score_matrix = []
- tmp_array = []
+ if (len(score) > number_of_points):
+ last_avg = 0
+ for j in xrange(number_of_points-1,len(score)):
+ last_avg += score[j]
+ last_avg = last_avg/(len(score)-number_of_points+1)
+ else:
+ last_avg = score[-1]
- # data for R boxplot
- title_keys = score_hash.keys()
+ score_points_tmp = []
+
+ for k in range(number_of_points-1):
+ score_points_tmp.append(score[k])
+
+ score_points_tmp.append(last_avg)
+
+ return score_points_tmp
+
+def __main__():
+ # I/O
+ infile_score_name = sys.argv[1].strip()
+ outfile_R_name = sys.argv[2].strip()
+
+ # unzip infile
+ tmp_score_dir = ''
+ score_file_list = []
+ if (zipfile.is_zipfile(infile_score_name)): (tmp_score_dir, score_file_list) = unzip_files(infile_score_name)
+ else: score_file_list = [infile_score_name]
+
+ # detect whether it's tabular or fasta format
+ seq_method = ''
+ test_file = score_file_list[0]
+ test_fh = open(test_file,'r')
+ while seq_method == '':
+ read_scorefile = test_fh.readline()
+ if read_scorefile.startswith('#'):
+ continue
+ elif read_scorefile.startswith(">"):
+ seq_method = '454'
+ else:
+ seq_method = 'solexa'
+ test_fh.close()
+
+ # quantile array
+ quality_score = {}
+
+ # R
+ score_points = []
+ score_matrix = []
+
+ tmp_read_length = 0
+ tmp_varied_length = False
+
+ test_file = score_file_list[0]
if (seq_method == 'solexa'):
- for read_title in title_keys:
- score_points.append(score_hash[(read_title)])
- else:
- for read_title in title_keys:
- tmp_score = '0 ' + score_hash[(read_title)]
- tmp_list = tmp_score.split()
- read_length = len(tmp_list)
- if (read_length > 100):
- step = int(math.floor((read_length-1)*1.0/number_of_points))
- score = []
- point = 1
- point_sum = 0
- step_average = 0
- for i in xrange(1,read_length):
- if (i < (point * step)):
- point_sum += int(tmp_list[i])
- step_average += 1
+ tmp_score = []
+ for i, line in enumerate(open(test_file)):
+ line = line.rstrip('\r\n')
+ tmp_score = line.split('\t')
+ if (tmp_read_length == 0): tmp_read_length = len(tmp_score)
+ if (tmp_read_length != len(tmp_score)):
+ tmp_varied_length = True
+ else:
+ # skip the last fasta sequence
+ tmp_score = ''
+ for i, line in enumerate(open(test_file)):
+ line = line.rstrip('\r\n')
+ if line.startswith('>'):
+ if len(tmp_score) > 0:
+ tmp_score = tmp_score.split()
+ if (tmp_read_length == 0): tmp_read_length = len(tmp_score)
+ if (tmp_read_length != len(tmp_score)):
+ tmp_varied_length = True
+ tmp_score = ''
+ else:
+ tmp_score = tmp_score + ' ' + line
+
+ if (tmp_varied_length): number_of_points = 20
+ else: number_of_points = tmp_read_length
+
+ # data for R boxplot
+ # data for quantile
+ if (seq_method == 'solexa'):
+ for score_file in score_file_list:
+ for i, line in enumerate(open(score_file)):
+ line = line.rstrip('\r\n')
+ each_loc = line.split('\t')
+ tmp_array = []
+ 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)
+ tmp_array.append(big)
+ score_points.append(tmp_array)
+ # quantile
+ for j,k in enumerate(tmp_array):
+ if quality_score.has_key((j,k)):
+ quality_score[(j, k)] += 1
else:
- point_avg = point_sum*1.0/step_average
- score.append(point_avg)
- point += 1
- point_sum = 0
- step_average = 0
-
- if (step_average > 0):
- point_avg = point_sum*1.0/step_average
- score.append(point_avg)
-
- if (len(score) > number_of_points):
- last_avg = 0
- for j in xrange(number_of_points-1,len(score)):
- last_avg += score[j]
- last_avg = last_avg/(len(score)-number_of_points+1)
- else:
- last_avg = score[-1]
- score_points_tmp = []
- for k in range(number_of_points-1):
- score_points_tmp.append(score[k])
- score_points_tmp.append(last_avg)
- score_points.append(score_points_tmp)
- score_points_tmp = []
+ quality_score[(j, k)] = 1
+ else:
+ tmp_score = ''
+ for score_file in score_file_list:
+ for i, line in enumerate(open(score_file)):
+ if line.startswith('>'):
+ if len(tmp_score) > 0:
+ tmp_score = ['0'] + tmp_score.split()
+ read_length = len(tmp_score)
+ tmp_array = []
+ if (tmp_varied_length is False):
+ tmp_score.pop(0)
+ score_points.append(tmp_score)
+ tmp_array = tmp_score
+ elif (read_length > 100):
+ score_points_tmp = merge_to_20_datapoints(tmp_score)
+ score_points.append(score_points_tmp)
+ tmp_array = score_points_tmp
+ for j, k in enumerate(tmp_array):
+ if quality_score.has_key((j,k)):
+ quality_score[(j,k)] += 1
+ else:
+ quality_score[(j,k)] = 1
+ tmp_score = ''
+ else:
+ tmp_score = tmp_score + ' ' + line
+ if len(tmp_score) > 0:
+ tmp_score = ['0'] + tmp_score.split()
+ read_length = len(tmp_score)
+ if (tmp_varied_length is False):
+ tmp_score.pop(0)
+ score_points.append(tmp_score)
+ elif (read_length > 100):
+ score_points_tmp = merge_to_20_datapoints(tmp_score)
+ score_points.append(score_points_tmp)
+ tmp_array = score_points_tmp
+ for j, k in enumerate(tmp_array):
+ if quality_score.has_key((j,k)):
+ quality_score[(j,k)] += 1
+ else:
+ quality_score[(j,k)] = 1
+
+ """
+ # quantile
+ keys = quality_score.keys()
+ keys.sort()
+ for key in keys:
+ print key, quality_score[key]
+ """
+
# reverse the matrix, for R
+ tmp_array = []
for i in range(number_of_points-1):
for j in range(len(score_points)):
- tmp_array.append(score_points[j][i])
+ tmp_array.append(int(score_points[j][i]))
score_matrix.append(tmp_array)
tmp_array = []
@@ -207,9 +226,8 @@ def generate_boxplot_figure():
outfile_R_pdf = outfile_R_name
r.pdf(outfile_R_pdf)
- #title = infile_score_name.split('/')[-1]
title = "boxplot of quality scores"
- if (seq_method=='solexa'):
+ if (tmp_varied_length is False):
r.boxplot(score_matrix,xlab="location in read length",main=title)
else:
r.boxplot(score_matrix,xlab="percentage in read length",xaxt="n",main=title)
@@ -222,45 +240,12 @@ def generate_boxplot_figure():
r.axis(1,x_old_range,x_new_range)
r.dev_off()
-
- return 0
+ if (os.path.isdir(tmp_score_dir)):
+ for file_name in score_file_list:
+ os.remove(file_name)
+ os.removedirs(tmp_score_dir)
-# I/O
-infile_score_name = sys.argv[1].strip()
-outfile_R_name = sys.argv[2].strip()
+ r.quit(save = "no")
-# to unzip or not unzip file
-tmp_score_dir = ''
-score_file_list = []
-if (zipfile.is_zipfile(infile_score_name)): (tmp_score_dir, score_file_list) = unzip_files(infile_score_name)
-else: score_file_list = [infile_score_name]
-read_scorefile = read_input_files(score_file_list)
-if (os.path.isdir(tmp_score_dir)):
- for file_name in score_file_list:
- os.remove(file_name)
- os.removedirs(tmp_score_dir)
-
-# detect whether it's tabular or fasta format
-if read_scorefile[0].startswith(">"):
- seq_method = '454'
-else:
- seq_method = 'solexa'
-
-# detect whether it's a sequence file
-if re.match('[A-Z]', read_scorefile[1]):
- stop_err("this is not a quality score file.")
-
-# parse files
-if (seq_method == 'solexa'):
- score_hash = parse_solexa_files(read_scorefile, 'score')
- number_of_points = 36
-else: # the other two are both fasta format
- score_hash = parse_fasta_format(read_scorefile)
- number_of_points = 20
-
-print seq_method, number_of_points
-
-# R
-status = generate_boxplot_figure()
-r.quit(save = "no")
+if __name__=="__main__":__main__()
diff --git a/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.py b/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.py
index 71c92e95569..26a71e9a4d3 100644
--- a/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.py
+++ b/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.py
@@ -1,49 +1,7 @@
#! /usr/bin/python
-"""
-Galaxy
-to commit
-1. filter: select read length longer than a threshold without trimming (minimum score = 0?)
-4. incorporate greg's changes.
-4.2 for i, line in enumerate(file(filename)):
-4.4 stop passing read_seqfile and read_scorefile around
-4.5 use --> if __name__ == "__main__": __main__()
-5. use galaxy extention
-6. a better trimming process for 454 (homopolymers)
-6.1 to keep the longest homoloplymers
-Input:
-1. sequence file(s): zip or text file --> for 454 and Solexa
-2. quality score file(s): zip or text file --> for 454 and Solexa
-3. trace file(s): zip or single scf file --> for Sanger
-4. threshold to trim the sequence
-5. threshold of the trimmed sequence length to be reported
-Output:
-trimmed sequence in one file
-----
-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)
@@ -51,23 +9,11 @@ def stop_err(msg):
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_dir_name = tempfile.mkdtemp() + '/'
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
@@ -87,98 +33,49 @@ def unzip_files(file_name):
return new_file_dir, read_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])
+ return
-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_values = each_line.split()
- read = tmp_values[-1]
- 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 trim_seq(seq, score, specific_argument, trim_score, threshold):
+ seq_method = '454'
+ trim_after_this_position = 0
+ keep_homopolymers = 'no'
-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 = ''
+ # trim after a certain position
+ if specific_argument.isdigit():
+ keep_homopolymers = 'no'
+ trim_after_this_position = int(specific_argument)
+ if (trim_after_this_position > 0 and trim_after_this_position < len(seq)):
+ seq = seq[0:trim_after_this_position]
+ else:
+ keep_homopolymers = specific_argument
+
+ new_trim_seq = ''
+
+ for i in range(len(seq)):
+ if (i >= len(score)): score.append(0)
+ if (int(score[i]) >= trim_score):
+ pass_nuc = seq[i:(i+1)]
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 trim_seq():
- # trim sequence
- title_keys = seq_hash.keys()
- for read_title in title_keys:
- tmp_seq = seq_hash[(read_title)]
- if (seq_method == 'solexa'):
- tmp_score = score_hash[(read_title)]
- else:
- tmp_score = score_hash[(read_title)].split()
-
- trim_seq = ''
-
- for i in range(len(tmp_seq)):
- if (int(tmp_score[i]) > threshold_trim):
- pass_nuc = tmp_seq[i:(i+1)]
+ # keep homopolymers?
+ if keep_homopolymers == 'yes' and (((i == 0) or (seq[i:(i+1)].lower() == seq[(i-1):i].lower()))):
+ pass_nuc = seq[i:(i+1)]
else:
- pass_nuc = ' '
- trim_seq = trim_seq + pass_nuc
+ pass_nuc = ' '
+ new_trim_seq = new_trim_seq + pass_nuc
# find the max substrings
- segments = trim_seq.split()
+ segments = new_trim_seq.split()
max_segment = ''
len_max_segment = 0
- if (threshold_report == 0):
+ if (threshold == 0):
for each_segment in segments:
if (len_max_segment < len(each_segment)):
max_segment = each_segment + ','
@@ -187,14 +84,13 @@ def trim_seq():
max_segment = max_segment + each_segment + ','
else:
for each_segment in segments:
- if (len(each_segment) >= threshold_report):
+ if (len(each_segment) >= threshold):
max_segment = max_segment + each_segment + ','
-
- trim_seq_hash[(read_title)] = max_segment[0:-1]
-
- return trim_seq_hash
+
+ return max_segment[0:-1]
def __main__():
+
# I/O
seq_method = sys.argv[1].strip().lower()
try:
@@ -202,128 +98,131 @@ def __main__():
except:
stop_err("Invalid value for minimal quality score")
try:
- threshold_report= int(sys.argv[3].strip())
+ threshold_report = int(sys.argv[3].strip())
except:
stop_err("Invalid value for minimal sequence length")
outfile_seq_name = sys.argv[4].strip()
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:
+ trim_after_this_position = int(special_argument)
- # to unzip or not unzip file
+ # unzip files
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)
tmp_score_dir = ''
score_file_list = []
if (zipfile.is_zipfile(infile_score_name)): (tmp_score_dir, score_file_list) = unzip_files(infile_score_name)
else: score_file_list = [infile_score_name]
- read_scorefile = read_input_files(score_file_list)
- if (os.path.isdir(tmp_score_dir)):
- for file_name in score_file_list:
- os.remove(file_name)
- os.removedirs(tmp_score_dir)
- # parse files
- if (seq_method == 'solexa'):
- seq_hash = parse_solexa_files(read_seqfile, 'seq')
- score_hash = parse_solexa_files(read_scorefile, 'score')
- if (threshold_report > 36):
- print "the read length from solexa is 36bp, your choice of read length is not valid.\n"
- print "the threshold is set to zero (return the longest substring).\n"
- threshold_report = 0
- number_of_points = 36
- else: # the other two are both fasta format
- seq_hash = parse_fasta_format(read_seqfile)
- score_hash = parse_fasta_format(read_scorefile)
- number_of_points = 20
-
- # trim sequence
- trim_seq_hash = trim_seq()
+ # open both sequence and score file
+ score_dirname_list = []
+ score_rootname_list = []
+ score_extname_list = []
+ for score_file in score_file_list:
+ (score_dirname, score_basename) = os.path.split(score_file)
+ (score_rootname, score_extname) = os.path.splitext(score_basename)
+ score_dirname_list.append(score_dirname)
+ score_rootname_list.append(score_rootname)
+ score_extname_list.append(score_extname)
- # output trimmed sequence to a fasta file
outfile_seq = open(outfile_seq_name,'w')
- title_keys = seq_hash.keys()
- for read_title in title_keys:
- tmp_seq = seq_hash[(read_title)]
- tmp_trim_seq = trim_seq_hash[(read_title)].split(',')
- if (len(tmp_trim_seq) > 1):
- for i in range(len(tmp_trim_seq)):
- print >> outfile_seq, "%s %d\n%s" % (read_title, i, tmp_trim_seq[i])
- elif (len(tmp_trim_seq[0]) > 0):
- print >> outfile_seq, "%s\n%s" % (read_title, tmp_trim_seq[0])
- outfile_seq.close()
-
-#if __name__ == "__main__": __main__()
-if 1:
- # I/O
- seq_method = sys.argv[1].strip().lower()
- try:
- threshold_trim = int(sys.argv[2].strip())
- except:
- stop_err("Invalid value for minimal quality score")
- try:
- threshold_report= int(sys.argv[3].strip())
- except:
- stop_err("Invalid value for minimal sequence length")
- outfile_seq_name = sys.argv[4].strip()
- infile_seq_name = sys.argv[5].strip()
- infile_score_name = sys.argv[6].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)
+ for seq_file in seq_file_list:
+ if (len(seq_file_list) > 1):
+ (seq_dirname, seq_basename) = os.path.split(seq_file)
+ (seq_rootname, seq_extname) = os.path.splitext(seq_basename)
+ if (seq_rootname in score_rootname_list):
+ file_index = score_rootname_list.index(seq_rootname)
+ score_file = os.path.join(score_dirname_list[file_index], seq_rootname+score_extname_list[file_index])
+ else:
+ score_file = score_file_list[0]
+ if (os.path.exists(seq_file) and os.path.exists(score_file)):
+ # read one sequence
+ to_find_score = True
+ seq = None
+ score = None
+ score_fh = open(score_file,'r')
+ if seq_method == '454':
+ for i, line in enumerate (open (seq_file) ):
+ line = line.rstrip('\r\n')
+ if (line.startswith('>')):
+ if seq:
+ to_find_score = True
+ score = None
+ while to_find_score:
+ score_line = score_fh.readline().rstrip('\r\n')
+ if (score_line.startswith('>')):
+ if score:
+ score = score.split()
+ new_trim_seq_segments = trim_seq(seq, score, special_argument, threshold_trim, threshold_report)
+ # output trimmed sequence to a fasta file
+ segments = new_trim_seq_segments.split(',')
+ show_output(outfile_seq, seq_title, segments)
+ to_find_score = False
+ score = None
+ else:
+ if not score: score = score_line
+ else:
+ score = score + ' ' + score_line
+ seq_title = line
+ seq = None
+ else:
+ if not seq: seq = line
+ else:
+ seq = seq + line
+ if seq:
+ score = None
+ while score_line:
+ score_line = score_fh.readline().rstrip('\r\n')
+ if ( not score_line.startswith('>')):
+ if not score: score = score_line
+ else:
+ score = score + ' ' + score_line
+ if score:
+ score = score.split()
+ new_trim_seq_segments = trim_seq(seq, score, special_argument, threshold_trim, threshold_report)
+
+ # output trimmed sequence to a fasta file
+ segments = new_trim_seq_segments.split(',')
+ show_output(outfile_seq, seq_title, segments)
+ else: # Solexa format
+ for i, line in enumerate (open (seq_file) ):
+ line = line.rstrip('\r\n')
+ seq_title = '>' + str(i)
+ seq = line.split()[-1]
+ score = score_fh.readline()
+ each_loc = score.split('\t')
+ score = []
+ 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)
+ score.append(big)
+ new_trim_seq_segments = trim_seq(seq, score, special_argument, threshold_trim, threshold_report)
+ # output trimmed sequence to a fasta file
+ segments = new_trim_seq_segments.split(',')
+ show_output(outfile_seq, seq_title, segments)
+ score_fh.close()
+ outfile_seq.close()
+
if (os.path.isdir(tmp_seq_dir)):
for file_name in seq_file_list:
os.remove(file_name)
os.removedirs(tmp_seq_dir)
-
- tmp_score_dir = ''
- score_file_list = []
- if (zipfile.is_zipfile(infile_score_name)): (tmp_score_dir, score_file_list) = unzip_files(infile_score_name)
- else: score_file_list = [infile_score_name]
- read_scorefile = read_input_files(score_file_list)
+
if (os.path.isdir(tmp_score_dir)):
for file_name in score_file_list:
os.remove(file_name)
os.removedirs(tmp_score_dir)
- # parse files
- if (seq_method == 'solexa'):
- seq_hash = parse_solexa_files(read_seqfile, 'seq')
- score_hash = parse_solexa_files(read_scorefile, 'score')
- if (threshold_report > 36):
- print "the read length from solexa is 36bp, your choice of read length is not valid.\n"
- print "the threshold is set to zero (return the longest substring).\n"
- threshold_report = 0
- number_of_points = 36
- else: # the other two are both fasta format
- seq_hash = parse_fasta_format(read_seqfile)
- score_hash = parse_fasta_format(read_scorefile)
- number_of_points = 20
-
- # trim sequence
- trim_seq_hash = trim_seq()
-
- # output trimmed sequence to a fasta file
- outfile_seq = open(outfile_seq_name,'w')
- title_keys = seq_hash.keys()
- for read_title in title_keys:
- tmp_seq = seq_hash[(read_title)]
- tmp_trim_seq = trim_seq_hash[(read_title)].split(',')
- if (len(tmp_trim_seq) > 1):
- for i in range(len(tmp_trim_seq)):
- print >> outfile_seq, "%s %d\n%s" % (read_title, i, tmp_trim_seq[i])
- elif (len(tmp_trim_seq[0]) > 0):
- print >> outfile_seq, "%s\n%s" % (read_title, tmp_trim_seq[0])
- outfile_seq.close()
-
+if __name__ == "__main__": __main__()
\ No newline at end of file
diff --git a/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml b/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml
index 4616cc7dfcf..124ce361282 100644
--- a/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml
+++ b/tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml
@@ -1,8 +1,8 @@
-
+
based on quality scores
-#if $sequencing_method_choice.sequencer=="454":#short_reads_trim_seq.py 454 $trim $length $output1 $input1 $input2
-#else:#short_reads_trim_seq.py solexa $trim $length $output1 $input1 $input2
+#if $sequencing_method_choice.sequencer=="454":#short_reads_trim_seq.py 454 $trim $length $output1 $sequencing_method_choice.input1 $sequencing_method_choice.input2 $sequencing_method_choice.input3
+#else:#short_reads_trim_seq.py solexa $trim $length $output1 $sequencing_method_choice.input1 $sequencing_method_choice.input2 $sequencing_method_choice.input3
#end if
@@ -10,16 +10,21 @@
-
+
-
+
+
+
+
+
+
@@ -32,22 +37,24 @@
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
@@ -55,6 +62,12 @@
.. class:: warningmark
Please select correct sequencing method.
+
+
-----
@@ -64,7 +77,7 @@ This tool takes raw data from sequencing facility and outputs a multi-fasta file
Accept the following two file formats:
-- fasta: for **454** data
+- fasta: for **454** and **SOLiD** data
- tabular: for **Solexa** data. Only take the maximal value from each 4 quality scores.
If *minimal trimmed length to report* is set to a number larger than zero, the tool ouputs any substrings that are longer than this threshold.
@@ -73,17 +86,15 @@ If there are more than one substring, the tool outputs all of them in separated
-----
-454 data::
+454 and SOLiD data::
>seq1
CCGGTATCCG
-454 quality score::
-
>seq1
33 19 34 25 28 28 28 32 15 34
-454 Output (trimmed by quality score 20 and returned the longest substring)::
+454 and SOLiD Output (trimmed by quality score 20 and returned the longest substring)::
>seq1
GGTATC
@@ -92,12 +103,11 @@ Solexa data::
5 300 902 419 TTCTTACCTATTAGTGGTTGAACATC
-
Solexa quality score (only showed the maximum value)::
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40
-Solexa Output::
+Solexa Output (trimmed by quality score 20 and returned the longest substring)::
>15774
TTCTTACCTATTAGTGGTTGAACA