From 5b4e5a772bcc53c3d9e62741f2938c166990a607 Mon Sep 17 00:00:00 2001 From: Wen-Yu Chung Date: Mon, 25 Feb 2008 18:42:39 +0000 Subject: [PATCH] modify data structure used in short read tools. update functional test data. --- tool_conf.xml.sample | 5 +- .../1.0.0/short_reads_figure_score.py | 379 +++++++++--------- .../trim_reads/1.0.0/short_reads_trim_seq.py | 379 +++++++----------- .../trim_reads/1.0.0/short_reads_trim_seq.xml | 50 ++- 4 files changed, 354 insertions(+), 459 deletions(-) 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