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

modify data structure used in short read tools.

Browse Source 
update functional test data.
This commit is contained in:
Wen-Yu Chung
2008-02-25 18:42:39 +00:00
parent 737f5b49cd
commit 5b4e5a772b
4 changed files with 354 additions and 459 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tool_conf.xml.sample
+3 -2
View File
@@ -318,6 +318,7 @@
</section>
<section name="Short Read Analysis" id="short_read_analysis">
<tool file="metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml" />
<tool file="metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.xml" />
</section>
<tool file="metag_tools/blat_wrapper/1.0.0/blat_wrapper.xml" />
<tool file="metag_tools/megablast_wrapper/1.0.0/megablast_wrapper.xml" />
</section>
</toolbox>
tools/metag_tools/quality_score_distribution/1.0.0/short_reads_figure_score.py
+182 -197
View File
@@ -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__()
tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.py
+139 -240
View File
@@ -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__()
tools/metag_tools/trim_reads/1.0.0/short_reads_trim_seq.xml
+30 -20
View File
@@ -1,8 +1,8 @@
<tool id="trim_reads" name="Trim Reads">
<tool id="trim_reads" name="Trim Reads" version="2.0.0">
<description>based on quality scores</description>
<command interpreter="python">#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
<command interpreter="python">#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
</command>
@@ -10,16 +10,21 @@
<page>
<conditional name="sequencing_method_choice">
<param name="sequencer" type="select" label="short read sequencing method">
<option value="454">454</option>
<option value="454">454 and SOLiD</option>
<option value="Solexa">Solexa</option>
</param>
<when value="454">
<param name="input1" type="data" format="fasta,txtseq.zip" label="Sequence file" />
<param name="input2" type="data" format="txt,txtseq.zip" label="Score file" />
<param name="input2" type="data" format="txt,txtseq.zip" label="Score file" />
<param name="input3" type="select" label="To keep homopolymers" >
<option value="yes">Yes</option>
<option value="no">No</option>
</param>
</when>
<when value="Solexa">
<param name="input1" type="data" format="tabular,txtseq.zip" label="Sequence file" />
<param name="input2" type="data" format="tabular,txtseq.zip" label="Score file" />
<param name="input3" type="integer" size="5" value="0" label="Trim reads after this position" help="0 as do nothing" />
</when>
</conditional>
<param name="trim" type="integer" size="5" value="20" label="Minimal quality score" />
@@ -32,22 +37,24 @@
</outputs>
<tests>
<test>
<param name="sequencer" value="Solexa" />
<param name="input1" value="solexa.fna" />
<param name="input2" value="solexa.qual" />
<param name="trim" value="20" />
<param name="length" value="0" />
<output name="output1" file="solexaTest.fa" />
</test>
<test>
<param name="sequencer" value="454" />
<param name="input1" value="454.fna" />
<param name="input2" value="454.qual" />
<param name="input3" value="no" />
<param name="trim" value="20" />
<param name="length" value="0" />
<output name="output1" file="454Test.fa" />
</test>
<test>
<param name="sequencer" value="Solexa" />
<param name="input1" value="solexa.fna" />
<param name="input2" value="solexa.qual" />
<param name="input3" value="0" />
<param name="trim" value="20" />
<param name="length" value="0" />
<output name="output1" file="solexaTest.fa" />
</test>
</tests>
<help>
@@ -55,6 +62,12 @@
.. class:: warningmark
Please select correct sequencing method.
<!--
.. class:: warningmark
Your **Quality Score** file needs to be a specific Galaxy **qualityscore** file format. Please click on the pencil symbol in the history panel to change the file format.
-->
-----
@@ -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::
&gt;seq1
CCGGTATCCG
454 quality score::
&gt;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)::
&gt;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)::
&gt;15774
TTCTTACCTATTAGTGGTTGAACA