mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
add a section: metagenomics
tools for running phred(need to change path later), trim sequence, generate boxplot and an array of length of the reads. test data are added. trim sequence and array of length passed functional test. phred (2 ouput files) and boxplot (pdf file) didn't pass functional test yet.
This commit is contained in:
Wen-Yu Chung
File diff suppressed because it is too large
Load Diff
@@ -294,4 +294,10 @@
|
||||
<tool file="emboss_5/emboss_wordcount.xml" />
|
||||
<tool file="emboss_5/emboss_wordmatch.xml" />
|
||||
</section> -->
|
||||
<section name="Metagenomics" id="metagenomics">
|
||||
<tool file="metag_tools/short_reads_run_phred.xml" />
|
||||
<tool file="metag_tools/short_reads_trim_seq.xml" />
|
||||
<tool file="metag_tools/short_reads_figure_score.xml" />
|
||||
<tool file="metag_tools/short_reads_figure_length.xml" />
|
||||
</section>
|
||||
</toolbox>
|
||||
|
||||
@@ -0,0 +1,209 @@
|
||||
#! /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:
|
||||
sequence file(s): zip or text file --> for 454 and Solexa
|
||||
Output:
|
||||
pdf files from R, length histogram --> before and after trim
|
||||
----
|
||||
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()
|
||||
for read_title in title_keys:
|
||||
tmp_seq = seq_hash[(read_title)]
|
||||
length_before_trim.append(len(tmp_seq))
|
||||
print >> outfile, 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")
|
||||
@@ -0,0 +1,44 @@
|
||||
<tool name="reads length array" id="length_array_short_reads">
|
||||
<description>on sequence file</description>
|
||||
|
||||
<command interpreter="python">short_reads_figure_length.py $input1 $output1</command>
|
||||
|
||||
<inputs>
|
||||
<page>
|
||||
<param name="input1" type="data" format="fasta,tabular,txtseq.zip" label="Sequence file" />
|
||||
</page>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data name="output1" format="text" />
|
||||
</outputs>
|
||||
|
||||
<tests>
|
||||
<test>
|
||||
<param name="input1" value="solexa.fna" />
|
||||
<output name="output1" file="solexaLength.txt" />
|
||||
</test>
|
||||
<test>
|
||||
<param name="input1" value="454.fna" />
|
||||
<output name="output1" file="454Length.txt" />
|
||||
</test>
|
||||
</tests>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
For **Sanger**, **454 or Solexa** sequencing technology only.
|
||||
|
||||
-----
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool takes a multi-fasta sequence file and outputs a text file:
|
||||
|
||||
- an array of the lengths in the file
|
||||
|
||||
To see the distribution, please use the Histogram tool in Graph/Display data.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
@@ -0,0 +1,266 @@
|
||||
#! /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)
|
||||
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_boxplot_figure():
|
||||
# R module and code
|
||||
|
||||
score_matrix = []
|
||||
tmp_array = []
|
||||
|
||||
# data for R boxplot
|
||||
title_keys = score_hash.keys()
|
||||
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
|
||||
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 = []
|
||||
|
||||
# reverse the matrix, for R
|
||||
for i in range(number_of_points-1):
|
||||
for j in range(len(score_points)):
|
||||
tmp_array.append(score_points[j][i])
|
||||
score_matrix.append(tmp_array)
|
||||
tmp_array = []
|
||||
|
||||
# generate pdf figures
|
||||
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'):
|
||||
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)
|
||||
x_old_range = []
|
||||
x_new_range = []
|
||||
step = 100/number_of_points
|
||||
for i in xrange(0,100,step):
|
||||
x_old_range.append((i/step))
|
||||
x_new_range.append(i)
|
||||
r.axis(1,x_old_range,x_new_range)
|
||||
|
||||
r.dev_off()
|
||||
|
||||
return 0
|
||||
|
||||
|
||||
# I/O
|
||||
infile_score_name = sys.argv[1].strip()
|
||||
outfile_R_name = sys.argv[2].strip()
|
||||
|
||||
# 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")
|
||||
@@ -0,0 +1,42 @@
|
||||
<tool name="quality score distribution" id="score_boxplot_short_reads">
|
||||
<description>on a boxplot</description>
|
||||
|
||||
<command interpreter="python">short_reads_figure_score.py $input1 $output1 </command>
|
||||
|
||||
<inputs>
|
||||
<page>
|
||||
<param name="input1" type="data" format="fasta,tabular,txtseq.zip" label="Quality score file" />
|
||||
</page>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data name="output1" format="pdf" />
|
||||
</outputs>
|
||||
<!--
|
||||
<tests>
|
||||
<test>
|
||||
<param name="input1" value="solexa.qual" />
|
||||
<output name="output1" file="solexaScore.pdf" />
|
||||
</test>
|
||||
<test>
|
||||
<param name="input1" value="454.qual" />
|
||||
<output name="output1" file="454Score.pdf" />
|
||||
</test>
|
||||
</tests>
|
||||
-->
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
For **Sanger**, **454 or Solexa** sequencing technology only.
|
||||
|
||||
-----
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool takes quality score files and outputs one file:
|
||||
|
||||
- figures in pdf file showing score distribution.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
@@ -0,0 +1,108 @@
|
||||
#! /usr/bin/python
|
||||
"""
|
||||
Galaxy
|
||||
to commit
|
||||
2. use galaxy extention
|
||||
Input:
|
||||
trace file(s): zip or single scf file --> for Sanger
|
||||
Output:
|
||||
1. sequcence file
|
||||
2. quality score file
|
||||
----
|
||||
Wen-Yu Chung
|
||||
"""
|
||||
import os, sys, math, tempfile, zipfile, re
|
||||
|
||||
# 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 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_trace_files(infile_name):
|
||||
|
||||
trace_file = []
|
||||
tmp_dir = ''
|
||||
seq_file_name = tempfile.mkstemp()[1] #tempfile.NamedTemporaryFile(dir=database_tmp).name
|
||||
score_file_name = tempfile.mkstemp()[1] #tempfile.NamedTemporaryFile(dir=database_tmp).name
|
||||
file_list_name = tempfile.mkstemp()[1] #tempfile.NamedTemporaryFile(dir=database_tmp).name
|
||||
file_list = open(file_list_name,'w')
|
||||
|
||||
if (zipfile.is_zipfile(infile_name)):
|
||||
(tmp_dir, trace_file) = unzip_files(infile_name)
|
||||
else:
|
||||
trace_file = [infile_name]
|
||||
for i in trace_file:
|
||||
print >> file_list,i
|
||||
file_list.close()
|
||||
|
||||
# call phred
|
||||
phred_command = '~/phred/phred -if ' + file_list_name + ' -sa ' + seq_file_name + ' -qa ' + score_file_name + ' -process_nomatch 2>/dev/null'
|
||||
os.system(phred_command)
|
||||
|
||||
# remove every scf file and the file_list itself
|
||||
if (os.path.isdir(tmp_dir)):
|
||||
for i in trace_file:
|
||||
os.remove(i)
|
||||
os.removedirs(tmp_dir)
|
||||
os.remove(file_list_name)
|
||||
|
||||
return seq_file_name, score_file_name
|
||||
|
||||
def __main__():
|
||||
# I/O
|
||||
infile_name = sys.argv[1].strip()
|
||||
outfile_seq_name = sys.argv[2].strip()
|
||||
outfile_score_name = sys.argv[3].strip()
|
||||
|
||||
(infile_seq_name, infile_score_name)=parse_trace_files(infile_name)
|
||||
os.system("mv " + infile_seq_name + " " + outfile_seq_name)
|
||||
os.system("mv " + infile_score_name + " " + outfile_score_name)
|
||||
|
||||
if __name__ == "__main__" : __main__()
|
||||
@@ -0,0 +1,27 @@
|
||||
<tool name="run phred on trace file" id="run_phred_on_trace_file">
|
||||
<description>return sequence and score files</description>
|
||||
<command interpreter="python">short_reads_run_phred.py $input1 $output1 $output2 </command>
|
||||
<inputs>
|
||||
<param name="input1" type="data" format="ab1,binseq.zip,scf" label="Trace file" />
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data name="output1" format="fasta" />
|
||||
<data name="output2" format="fasta" />
|
||||
</outputs>
|
||||
<help>
|
||||
.. class:: infomark
|
||||
|
||||
For **trace (abi and scf) file** only.
|
||||
|
||||
-----
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool takes raw data from sequencing facility and outputs two files:
|
||||
|
||||
- sequences in fasta format;
|
||||
|
||||
- quality score.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
@@ -0,0 +1,329 @@
|
||||
#! /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)
|
||||
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_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 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 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)]
|
||||
else:
|
||||
pass_nuc = ' '
|
||||
trim_seq = trim_seq + pass_nuc
|
||||
|
||||
# find the max substrings
|
||||
segments = trim_seq.split()
|
||||
max_segment = ''
|
||||
len_max_segment = 0
|
||||
|
||||
if (threshold_report == 0):
|
||||
for each_segment in segments:
|
||||
if (len_max_segment < len(each_segment)):
|
||||
max_segment = each_segment + ','
|
||||
len_max_segment = len(each_segment)
|
||||
elif (len_max_segment == len(each_segment)):
|
||||
max_segment = max_segment + each_segment + ','
|
||||
else:
|
||||
for each_segment in segments:
|
||||
if (len(each_segment) >= threshold_report):
|
||||
max_segment = max_segment + each_segment + ','
|
||||
|
||||
trim_seq_hash[(read_title)] = max_segment[0:-1]
|
||||
|
||||
return trim_seq_hash
|
||||
|
||||
def __main__():
|
||||
# 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)
|
||||
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__()
|
||||
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)
|
||||
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()
|
||||
|
||||
@@ -0,0 +1,72 @@
|
||||
<tool name="trim reads" id="trim_short_reads">
|
||||
<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
|
||||
#end if
|
||||
</command>
|
||||
|
||||
<inputs>
|
||||
<page>
|
||||
<conditional name="sequencing_method_choice">
|
||||
<param name="sequencer" type="select" label="short read sequencing method">
|
||||
<option value="454">fasta</option>
|
||||
<option value="Solexa">tabular</option>
|
||||
</param>
|
||||
<when value="454">
|
||||
<param name="input1" type="data" format="fasta,txtseq.zip" label="Sequence file" />
|
||||
<param name="input2" type="data" format="fasta,txtseq.zip" label="Score file" />
|
||||
</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" />
|
||||
</when>
|
||||
</conditional>
|
||||
<param name="trim" type="integer" size="5" value="20" label="Minimal quality score" />
|
||||
<param name="length" type="integer" size="5" value="0" label="Minimal length of trimmed reads to report" help="To output the read if its length is longer than this threshold. Use 0 to return the longest substring" />
|
||||
</page>
|
||||
</inputs>
|
||||
|
||||
<outputs>
|
||||
<data name="output1" format="fasta" />
|
||||
</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="trim" value="20" />
|
||||
<param name="length" value="0" />
|
||||
<output name="output1" file="454Test.fa" />
|
||||
</test>
|
||||
</tests>
|
||||
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
For **Sanger**, **454 or Solexa** sequencing technology only.
|
||||
Either trace file or both of the read file and quality score file must be provided.
|
||||
|
||||
-----
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool takes raw data from sequencing facility and outputs a multi-fasta file:
|
||||
|
||||
- trimmed sequences in fasta format;
|
||||
|
||||
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.
|
||||
If there are more than one substring, the tool outputs all of them in separated fasta format.
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
Reference in New Issue
Block a user