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