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.
This commit is contained in:
Wen-Yu Chung
2008-03-31 19:29:00 +00:00
parent 163e4bc200
commit 67986d25f6
5 changed files with 304 additions and 1 deletions
-1
View File
@@ -270,6 +270,5 @@
<tool file="metag_tools/megablast_xml_parser.xml" />
<tool file="metag_tools/blat_coverage_report.xml" />
<tool file="metag_tools/blat_mapping.xml" />
<tool file="metag_tools/convert_SOLiD_color2nuc.xml" />
</section>
</toolbox>
+98
View File
@@ -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__()
+42
View File
@@ -0,0 +1,42 @@
<tool id="blat2wig" name="Show Coverage of the Reads">
<description>in wiggle format</description>
<command interpreter="python">blat_mapping.py $input1 $output1</command>
<inputs>
<param name="input1" type="data" format="tabular" label="Alignment Result"/>
</inputs>
<outputs>
<data name="output1" format="wig"/>
</outputs>
<tests>
<test>
<param name="input1" value="blat_mapping_test1.txt" ftype="tabular" />
<output name="output1" file="blat_mapping_test1.out" />
</test>
</tests>
<help>
.. 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
</help>
</tool>
@@ -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__()
@@ -0,0 +1,72 @@
<tool id="color2nuc" name="Convert Color Space" version="1.0.0">
<description> to Nucleotides </description>
<command interpreter="python">convert_SOLiD_color2nuc.py $input1 $input2 $output1 </command>
<inputs>
<param name="input1" type="data" format="txt" label="SOLiD Color Coding File" />
<param name="input2" type="select" label="Keep Prefix Nucleotide">
<option value="yes">Yes</option>
<option value="no">No</option>
</param>
</inputs>
<outputs>
<data name="output1" format="fasta" />
</outputs>
<!--
<tests>
<test>
<param name="input1" value="convert_SOLiD_color2nuc_test1.txt" ftype="txt" />
<param name="input2" value="no" />
<output name="output1" file="convert_SOLiD_color2nuc_test1.out" />
</test>
</tests>
-->
<help>
.. 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::
&gt;seq1
A013
&gt;seq2
T011213122200221123032111221021210131332222101
- If you would like to **keep** the leading nucleotide::
&gt;seq1
AACG
&gt;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)::
&gt;seq1
ACG
&gt;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
</help>
</tool>