From 67986d25f60394f720f69a67bd0dff992cfe245a Mon Sep 17 00:00:00 2001 From: Wen-Yu Chung Date: Mon, 31 Mar 2008 19:29:00 +0000 Subject: [PATCH] sorry, hit the wrong button (want to click cancel but click ok instead). add two tools and functional test data for short reads. update tool_conf.xml.sample. --- tool_conf.xml.sample | 1 - tools/metag_tools/blat_mapping.py | 98 +++++++++++++++++++ tools/metag_tools/blat_mapping.xml | 42 ++++++++ tools/metag_tools/convert_SOLiD_color2nuc.py | 92 +++++++++++++++++ tools/metag_tools/convert_SOLiD_color2nuc.xml | 72 ++++++++++++++ 5 files changed, 304 insertions(+), 1 deletion(-) create mode 100644 tools/metag_tools/blat_mapping.py create mode 100644 tools/metag_tools/blat_mapping.xml create mode 100644 tools/metag_tools/convert_SOLiD_color2nuc.py create mode 100644 tools/metag_tools/convert_SOLiD_color2nuc.xml diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 3f7fcd8d8a8..56f212dc2c2 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -270,6 +270,5 @@ - diff --git a/tools/metag_tools/blat_mapping.py b/tools/metag_tools/blat_mapping.py new file mode 100644 index 00000000000..ad02263b3eb --- /dev/null +++ b/tools/metag_tools/blat_mapping.py @@ -0,0 +1,98 @@ +#! /usr/bin/python + +import os, sys + +assert sys.version_info[:2] >= (2.4) + +def reverse_complement(s): + + complement_dna = {"A":"T", "T":"A", "C":"G", "G":"C", "a":"t", "t":"a", "c":"g", "g":"c", "N":"N", "n":"n" , ".":"."} + reversed_s = [] + for i in s: + reversed_s.append(complement_dna[i]) + reversed_s.reverse() + + return "".join(reversed_s) + + +def __main__(): + + nuc_index = {'a':0,'t':1,'c':2,'g':3,'n':4} + coverage = {} # key = (chrom, index) + + invalid_lines = 0 + invalid_chrom = 0 + + infile = sys.argv[1] + outfile = sys.argv[2] + + for i, line in enumerate(open(infile)): + + line = line.rstrip('\r\n') + fields = line.split() + + if line.startswith('#'): continue + if not line: continue + + if (len(fields) < 21): # standard number of pslx columns + invalid_lines += 1 + continue + if (not fields[0].isdigit()): + invalid_lines += 1 + continue + + chrom = fields[13] + + try: + assert chrom.startswith('chr') is True + except: + invalid_chrom += 1 + continue + + try: + block_count = int(fields[17]) + except: + invalid_lines += 1 + continue + + block_size = fields[18].split(',') + chrom_start = fields[20].split(',') + + + for j in range(block_count): + + try: + this_block_size = int(block_size[j]) + this_chrom_start = int(chrom_start[j]) + except: + continue + + # brut force coverage + for k in range(this_block_size): + cur_index = this_chrom_start+k + if coverage.has_key((chrom,cur_index)): + coverage[(chrom, cur_index)] += 1 + else: + coverage[(chrom, cur_index)] = 1 + + # generate a index file + outputfh = open(outfile, 'w') + keys = coverage.keys() + keys.sort() + previous_chrom = '' + for i in keys: + (chrom, location) = i + sum = coverage[(i)] + if (chrom != previous_chrom): + print >> outputfh, 'variableStep chrom=%s' %(chrom) + previous_chrom = chrom + print >> outputfh, location, sum + outputfh.close() + + if invalid_lines: + print >> sys.stdout, "Skip %d invalid lines. These lines could be headers or have fewer columns than standard output." %(invalid_lines) + + if invalid_chrom: + print >> sys.stdout, "Skip %d invalid lines with errors in chromosome id. The chromosome id must begin with \'chr\' to be correctly mapped to ucsc genome browser." + +if __name__ == '__main__': __main__() \ No newline at end of file diff --git a/tools/metag_tools/blat_mapping.xml b/tools/metag_tools/blat_mapping.xml new file mode 100644 index 00000000000..d093c67a638 --- /dev/null +++ b/tools/metag_tools/blat_mapping.xml @@ -0,0 +1,42 @@ + + in wiggle format + blat_mapping.py $input1 $output1 + + + + + + + + + + + + + + +.. class:: warningmark + +**TIP**. To generate acceptable files, please use alignment program **BLAT** with option **-out=pslx**. + +.. class:: warningmark + +**TIP**. Please edit the database information by click on the pencil symbol in your history panel. Select the corresponding genome build. + +----- + +**What it does** + + This tool takes **BLAT pslx** output and returns a wig-like file showing the number of reads (coverage) mapped at each chromosome location. Use **Graph/Display Data --> Build custom track** tool to show the coverage mapping in UCSC Genome Browser. + +----- + +**Example** + + Showing reads coverage on human chromosome 22 (partial result) in UCSC Genome Browser Custom Track (black lines) with SNPs: + + .. image:: ../static/images/blat_mapping_example.png + :width: 600 + + + diff --git a/tools/metag_tools/convert_SOLiD_color2nuc.py b/tools/metag_tools/convert_SOLiD_color2nuc.py new file mode 100644 index 00000000000..f679615660f --- /dev/null +++ b/tools/metag_tools/convert_SOLiD_color2nuc.py @@ -0,0 +1,92 @@ +#! /usr/bin/python +""" +convert SOLiD calor-base data to nucleotide sequence +example: T011213122200221123032111221021210131332222101 + TTGTCATGAGAAAGACAGCCGACACTCAAGTCAACGTATCTCTGGT +""" + +import sys, os + +assert sys.version_info[:2] >= (2.4) + +def stop_err(msg): + + sys.stderr.write(msg) + sys.stderr.write('\n') + sys.exit() + +def color2base(color_seq): + + first_nuc = ['A','C','G','T'] + code_matrix = {} + code_matrix['0'] = ['A','C','G','T'] + code_matrix['1'] = ['C','A','T','G'] + code_matrix['2'] = ['G','T','A','C'] + code_matrix['3'] = ['T','G','C','A'] + + overlap_nuc = '' + nuc_seq = '' + + seq_prefix = prefix = color_seq[0].upper() + color_seq = color_seq[1:] + + try: + assert (seq_prefix in first_nuc) is True + except: + stop_err('The leading nucleotide is invalid. Must be one of the four nucleotides: A, T, C, G.\nThe file contains a %s' %seq_prefix ) + + for code in color_seq: + + try: + assert (code in ['0','1','2','3']) is True + except: + stop_err('Expect digits (0, 1, 2, 3) in the color-cading data. File contains numbers other than the set.\nThe file contains a %s' %code) + + second_nuc = code_matrix[code] + overlap_nuc = second_nuc[first_nuc.index(prefix)] + nuc_seq += overlap_nuc + prefix = overlap_nuc + + return seq_prefix, nuc_seq + +def __main__(): + + infilename = sys.argv[1] + keep_prefix = sys.argv[2].lower() + outfilename = sys.argv[3] + + outfile = open(outfilename,'w') + prefix = '' + color_seq = '' + for i, line in enumerate(file(infilename)): + line = line.rstrip('\r\n') + if not line: continue + if line.startswith("#"): continue + if line.startswith(">"): + + if color_seq: + prefix, nuc_seq = color2base(color_seq) + + if keep_prefix == 'yes': + nuc_seq = prefix + nuc_seq + + print >> outfile, title + print >> outfile, nuc_seq + + title = line + color_seq = '' + else: + color_seq += line + + if color_seq: + prefix, nuc_seq = color2base(color_seq) + + if keep_prefix == 'yes': + nuc_seq = prefix + nuc_seq + + print >> outfile, title + print >> outfile, nuc_seq + + outfile.close() + +if __name__=='__main__': __main__() diff --git a/tools/metag_tools/convert_SOLiD_color2nuc.xml b/tools/metag_tools/convert_SOLiD_color2nuc.xml new file mode 100644 index 00000000000..94ffe30d86c --- /dev/null +++ b/tools/metag_tools/convert_SOLiD_color2nuc.xml @@ -0,0 +1,72 @@ + + to Nucleotides +convert_SOLiD_color2nuc.py $input1 $input2 $output1 + + + + + + + + + + + + + + +.. class:: warningmark + +**TIP**. The tool was designed for color space files generated from AB SOLiD sequencer. The file format must be fasta-like: the title starts with a ">" sign, and each color-space sequence starts with a leading nucleotide. + +----- + +**What it does** + + This tool convert a color-space sequence to nucleotides. The leading character must be one of the nucleotides (A, C, G, T). + +----- + +**Example** + +- If the color-space file looks like this:: + + >seq1 + A013 + >seq2 + T011213122200221123032111221021210131332222101 + +- If you would like to **keep** the leading nucleotide:: + + >seq1 + AACG + >seq2 + TTGTCATGAGAAAGACAGCCGACACTCAAGTCAACGTATCTCTGGT + +- If you **do not want to keep** the leading nucleotide (the length of nucleotide sequence will be one less than the color-space sequence):: + + >seq1 + ACG + >seq2 + TGTCATGAGAAAGACAGCCGACACTCAAGTCAACGTATCTCTGGT + +----- + +**SOLiD Color Coding Alignment matrix** + + Each di-nucleotide is represented by a single digit: 0 to 3. The matrix is symmetric, thus the leading nucleotide is necessary to determine the sequence (otherwise there are four possibilities). + + + .. image:: ../static/images/dualcolorcode.png + + + +