From 35a07c8b43e1c2800eba24d9b85696dc3f9f3aba Mon Sep 17 00:00:00 2001 From: Wen-Yu Chung Date: Tue, 25 Mar 2008 22:35:26 +0000 Subject: [PATCH] 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. --- .../short_reads_figure_high_quality_length.py | 188 ++++++++++++++++ ...short_reads_figure_high_quality_length.xml | 59 +++++ .../metag_tools/short_reads_figure_length.py | 206 ------------------ .../metag_tools/short_reads_figure_length.xml | 68 ------ 4 files changed, 247 insertions(+), 274 deletions(-) create mode 100644 tools/metag_tools/short_reads_figure_high_quality_length.py create mode 100644 tools/metag_tools/short_reads_figure_high_quality_length.xml delete mode 100644 tools/metag_tools/short_reads_figure_length.py delete mode 100644 tools/metag_tools/short_reads_figure_length.xml diff --git a/tools/metag_tools/short_reads_figure_high_quality_length.py b/tools/metag_tools/short_reads_figure_high_quality_length.py new file mode 100644 index 00000000000..2d6644399fd --- /dev/null +++ b/tools/metag_tools/short_reads_figure_high_quality_length.py @@ -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__() diff --git a/tools/metag_tools/short_reads_figure_high_quality_length.xml b/tools/metag_tools/short_reads_figure_high_quality_length.xml new file mode 100644 index 00000000000..7581ce659f7 --- /dev/null +++ b/tools/metag_tools/short_reads_figure_high_quality_length.xml @@ -0,0 +1,59 @@ + + of high quality score reads + +short_reads_figure_high_quality_length.py $input1 $output1 $input2 + + + + + + + + + + + + + +.. 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. + + + diff --git a/tools/metag_tools/short_reads_figure_length.py b/tools/metag_tools/short_reads_figure_length.py deleted file mode 100644 index dba6a07d493..00000000000 --- a/tools/metag_tools/short_reads_figure_length.py +++ /dev/null @@ -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") diff --git a/tools/metag_tools/short_reads_figure_length.xml b/tools/metag_tools/short_reads_figure_length.xml deleted file mode 100644 index eb4cbe16721..00000000000 --- a/tools/metag_tools/short_reads_figure_length.xml +++ /dev/null @@ -1,68 +0,0 @@ - -on sequence file - -short_reads_figure_length.py $input1 $output1 - - - - - - - - - - - - - - - - - - - - - - - - -**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. - -