Smart converter for SOLiD, which takes care of F3-R3 pairing problem. Requires python2.5

This commit is contained in:
Anton Nekrutenko
2010-01-05 14:33:39 -05:00
parent d4302ae6b9
commit 73db93f7bb
6 changed files with 430 additions and 36 deletions
+58 -33
View File
@@ -6,15 +6,26 @@
<meta name="generator" content="Docutils 0.3.9: http://docutils.sourceforge.net/" />
<link rel="stylesheet" href="style/base.css" type="text/css" />
<style type="text/css">
.quickie {
text-align: center;
background: black;
margin: 10px;
}
.current-quickie {
width: 100%;
background: black;
width: 300px;
background: white;
margin: auto;
}
.current-quickie img {
padding: 15px;
border: 1px solid #ccc;
margin: auto;
background-color: white;
}
.previous {
width: 100%;
overflow: auto;
@@ -30,42 +41,35 @@
margin-left: auto;
margin-right: auto;
}
</style>
<script type="text/javascript" src="http://galaxy.psu.edu/welcome_img/jquery.min.js"></script>
<script type="text/javascript" src="http://galaxy.psu.edu/welcome_img/jquery.cycle.all.2.72.js"></script>
<script type="text/javascript">
$(document).ready(function() {
$('.current-quickie').cycle({
fx: 'fade',
pause: 1,
speed: 100
});
});
</script>
</head>
<body>
<div class="document">
<div class="infomessagelarge">
<strong>We are hiring!</strong>
<hr>
Thanks to your support and an unprecedented level of usage, we are looking for an experienced software developer to join our team. For more information about this position, please, click <a target="_blank" href="https://www.bx.psu.edu/cgi-bin/trac.cgi/galaxy/wiki/PythonDeveloper">here</a>.
<div class="document">
<h3 align="center">Galaxy in 2010...</h3>
<div align="center" class="current-quickie">
<img src="http://galaxy.psu.edu/welcome_img/welcome_images.001.png" width="300" height="200" />
<img src="http://galaxy.psu.edu/welcome_img/welcome_images.002.png" width="300" height="200" />
<img src="http://galaxy.psu.edu/welcome_img/welcome_images.003.png" width="300" height="200" />
<img src="http://galaxy.psu.edu/welcome_img/welcome_images.004.png" width="300" height="200" />
<img src="http://galaxy.psu.edu/welcome_img/welcome_images.005.png" width="300" height="200" />
</div>
<br/>
<br>
<hr>
<div id="screencasts">
<h2>Introducing Galactic Quickies</h2>
<p>
Galactic quickies are <i>super-short</i> screencasts that are <i>always</i> under 5 minutes. We thought it may be a good way to spread the word about Galaxy's functionality while keeping the &quot;annoyance factor&quot; to the minimum. The quickies will be updated weekly.
</p>
<div class="current-quickie">
<table border="0" cellpadding="0" cellspacing="0" width="100%">
<tr>
<td align="center">
<a href="javascript:parent.show_in_overlay({url:'http://screencast.g2.bx.psu.edu/galaxy/quickie5_join/flow.html',width:640,height:500,scroll:'no'})">
<img src="images/qk/quickie5_lrg.png" border="0">
<br>
</a>
</td>
</tr>
</table>
</div>
<h3>Previous Quickies</h3>
<h3 align="center">Galaxy Quickies</h3>
<div class="previous" id="previous">
<table border="0" cellpadding="0" cellspacing="0" width="100%">
@@ -98,6 +102,27 @@
</div>
</a>
</td>
<td>
<a href="javascript:parent.show_in_overlay({url:'http://screencast.g2.bx.psu.edu/galaxy/quickie5_join/flow.html',width:640,height:500,scroll:'no'})">
<div class="quickie">
<img src="images/qk/quickie5_small.png" border="0">
</div>
</a>
</td>
<td>
<a href="javascript:parent.show_in_overlay({url:'http://screencast.g2.bx.psu.edu/galaxy/quickie6_share/flow.html',width:640,height:500,scroll:'no'})">
<div class="quickie">
<img src="images/qk/quickie6_small.png" border="0">
</div>
</a>
</td>
<td>
<a href="javascript:parent.show_in_overlay({url:'http://screencast.g2.bx.psu.edu/galaxy/quickie7_sr_beta/flow.html',width:640,height:500,scroll:'no'})">
<div class="quickie">
<img src="images/qk/quickie7_small.png" border="0">
</div>
</a>
</td>
</tr>
</table>
</div>
+1 -1
View File
@@ -185,7 +185,7 @@
<tool file="metag_tools/short_reads_figure_score.xml" />
<tool file="metag_tools/short_reads_trim_seq.xml" />
<label text="AB-SOLiD data" id="solid" />
<tool file="next_gen_conversion/solid_to_fastq.xml" />
<tool file="next_gen_conversion/solid2fastq.xml" />
<tool file="solid_tools/solid_qual_stats.xml" />
<tool file="solid_tools/solid_qual_boxplot.xml" />
</section>
+1 -1
View File
@@ -1,6 +1,6 @@
<tool id="cshl_fastq_to_fasta" name="FASTQ to FASTA">
<description>converter</description>
<command>gunzip -cf $input | fastq_to_fasta $SKIPN $RENAMESEQ -o $output -v </command>
<command>gunzip -cf $input | fastq_to_fasta -Q 33 $SKIPN $RENAMESEQ -o $output -v </command>
<inputs>
<param format="fastq" name="input" type="data" label="FASTQ Library to convert" />
+1 -1
View File
@@ -1,6 +1,6 @@
<tool id="cshl_fastx_collapser" name="Collapse">
<description>sequences</description>
<command>zcat -f '$input' | fastx_collapser -v -o '$output' </command>
<command>zcat -f '$input' | fastx_collapser -Q 33 -v -o '$output' </command>
<inputs>
<param format="fastqsolexa,fasta" name="input" type="data" label="Library to collapse" />
+210
View File
@@ -0,0 +1,210 @@
#!/usr/bin/env python
import sys
import string
import optparse
import tempfile
import sqlite3
def stop_err( msg ):
sys.stderr.write( msg )
sys.exit()
def solid2sanger( quality_string, min_qual = 0 ):
sanger = ""
quality_string = quality_string.rstrip( " " )
for qv in quality_string.split(" "):
try:
if int( qv ) < 0:
qv = '0'
if int( qv ) < min_qual:
return False
break
sanger += chr( int( qv ) + 33 )
except:
pass
return sanger
def Translator(frm='', to='', delete='', keep=None):
allchars = string.maketrans('','')
if len(to) == 1:
to = to * len(frm)
trans = string.maketrans(frm, to)
if keep is not None:
delete = allchars.translate(allchars, keep.translate(allchars, delete))
def callable(s):
return s.translate(trans, delete)
return callable
def merge_reads_qual( f_reads, f_qual, f_out, trim_name=False, out='fastq', double_encode = False, trim_first_base = False, pair_end_flag = '', min_qual = 0, table_name=None ):
# Reads from two files f_csfasta (reads) and f_qual (quality values) and produces output in three formats depending on out parameter,
# which can have three values: fastq, txt, and db
# fastq = fastq format
# txt = space delimited format with defline, reads, and qvs
# dp = dump data into sqlite3 db.
# IMPORTNAT! If out = db two optins must be provided:
# 1. f_out must be a db connection object initialized with sqlite3.connect()
# 2. table_name must be provided
if out == 'db':
cursor = f_out.cursor()
sql = "create table %s (name varchar(50) not null, read blob, qv blob)" % table_name
cursor.execute(sql)
lines = []
line = " "
while line:
for f in [ f_reads, f_qual ]:
line = f.readline().rstrip( '\n\r' )
while line.startswith( '#' ):
line = f.readline().rstrip( '\n\r' )
lines.append( line )
if lines[0].startswith( '>' ):
defline = lines[0][1:]
if trim_name and ( defline[ len( defline )-3: ] == "_F3" or defline[ len( defline )-3: ] == "_R3" ):
defline = defline[ : len( defline )-3 ]
else:
if trim_first_base:
lines[0] = lines[0][1:]
if double_encode:
de = Translator(frm="0123.", to="ACGTN")
lines[0] = de(lines[0])
qual = solid2sanger( lines[1], int( min_qual ) )
if qual:
if out == 'fastq':
f_out.write( "@%s%s\n%s\n+\n%s\n" % ( defline, pair_end_flag, lines[0], qual ) )
if out == 'txt':
f_out.write( '%s %s %s\n' % (defline, lines[0], qual ) )
if out == 'db':
cursor.execute('insert into %s values("%s","%s","%s")' % (table_name, defline, lines[0], qual ) )
lines = []
def main():
usage = "%prog --fr F3.csfasta --fq R3.csfasta --fout fastq_output_file [option]"
parser = optparse.OptionParser(usage=usage)
parser.add_option(
'--fr','--f_reads',
metavar="F3_CSFASTA_FILE",
dest='fr',
help='Name of F3 file with color space reads')
parser.add_option(
'--fq','--f_qual',
metavar="F3_QUAL_FILE",
dest='fq',
help='Name of F3 file with color quality values')
parser.add_option(
'--fout','--f3_fastq_output',
metavar="F3_OUTPUT",
dest='fout',
help='Name for F3 output file')
parser.add_option(
'--rr','--r_reads',
metavar="R3_CSFASTA_FILE",
dest='rr',
default = False,
help='Name of R3 file with color space reads')
parser.add_option(
'--rq','--r_qual',
metavar="R3_QUAL_FILE",
dest='rq',
default = False,
help='Name of R3 file with color quality values')
parser.add_option(
'--rout',
metavar="R3_OUTPUT",
dest='rout',
help='Name for F3 output file')
parser.add_option(
'-q','--min_qual',
dest='min_qual',
default = '-1000',
help='Minimum quality threshold for printing reads. If a read contains a single call with QV lower than this value, it will not be reported. Default is -1000')
parser.add_option(
'-t','--trim_name',
dest='trim_name',
action='store_true',
default = False,
help='Trim _R3 and _F3 off read names. Default is False')
parser.add_option(
'-f','--trim_first_base',
dest='trim_first_base',
action='store_true',
default = False,
help='Remove the first base of reads in color-space. Default is False')
parser.add_option(
'-d','--double_encode',
dest='de',
action='store_true',
default = False,
help='Double encode color calls as nucleotides: 0123. becomes ACGTN. Default is False')
options, args = parser.parse_args()
if not ( options.fout and options.fr and options.fq ):
parser.error("""
One or more of the three required paremetrs is missing:
(1) --fr F3.csfasta file
(2) --fq F3.qual file
(3) --fout name of output file
Use --help for more info
""")
fr = open ( options.fr , 'r' )
fq = open ( options.fq , 'r' )
f_out = open ( options.fout , 'w' )
if options.rr and options.rq:
rr = open ( options.rr , 'r' )
rq = open ( options.rq , 'r' )
if not options.rout:
parser.error("Provide the name for f3 output using --rout option. Use --help for more info")
r_out = open ( options.rout, 'w' )
db = tempfile.NamedTemporaryFile()
print db.name
try:
con = sqlite3.connect(db.name)
cur = con.cursor()
except:
stop_err('Cannot connect to %s\n') % db.name
merge_reads_qual( fr, fq, con, trim_name=options.trim_name, out='db', double_encode=options.de, trim_first_base=options.trim_first_base, min_qual=options.min_qual, table_name="f3" )
merge_reads_qual( rr, rq, con, trim_name=options.trim_name, out='db', double_encode=options.de, trim_first_base=options.trim_first_base, min_qual=options.min_qual, table_name="r3" )
cur.execute('create index f3_name on f3( name )')
cur.execute('create index r3_name on r3( name )')
cur.execute('select * from r3,f3 where f3.name = r3.name')
for item in cur:
f_out.write( "@%s%s\n%s\n+\n%s\n" % (item[0], "/1", item[1], item[2]) )
r_out.write( "@%s%s\n%s\n+\n%s\n" % (item[3], "/2", item[4], item[5]) )
else:
merge_reads_qual( fr, fq, f_out, trim_name=options.trim_name, out='fastq', double_encode = options.de, trim_first_base = options.trim_first_base, min_qual=options.min_qual )
f_out.close()
if __name__ == "__main__":
main()
+159
View File
@@ -0,0 +1,159 @@
<tool id="solid2fastq" name="Convert">
<description>SOLiD output to fastq</description>
<command interpreter="python">
#if $is_run.paired == "no" #solid2fastq.py --fr=$input1 --fq=$input2 --fout=$out_file1 -q $qual $trim_name $trim_first_base $double_encode
#elif $is_run.paired == "yes" #solid2fastq.py --fr=$input1 --fq=$input2 --fout=$out_file1 --rr=$input3 --rq=$input4 --rout=$out_file2 -q $qual $trim_name $trim_first_base $double_encode
#end if#
</command>
<inputs>
<param name="input1" type="data" format="csfasta" label="Select Forward reads"/>
<param name="input2" type="data" format="qualsolid" label="Select Forward qualities"/>
<conditional name="is_run">
<param name="paired" type="select" label="Is this a mate-pair run?">
<option value="no" selected="true">No</option>
<option value="yes">Yes</option>
</param>
<when value="yes">
<param name="input3" type="data" format="csfasta" label="Select Reverse reads"/>
<param name="input4" type="data" format="qualsolid" label="Select Reverse qualities"/>
</when>
<when value="no">
</when>
</conditional>
<param name="qual" label="Remove reads containing color qualities below this value" type="integer" value="0"/>
<param name="trim_name" type="select" label="Trim trailing &quot;_F3&quot; and &quot;_R3&quot; ?">
<option value="-t" selected="true">Yes</option>
<option value="">No</option>
</param>
<param name="trim_first_base" type="select" label="Trim first base?">
<option value="-f">Yes</option>
<option value="" selected="true">No</option>
</param>
<param name="double_encode" type="select" label="Double encode?">
<option value="-d">Yes</option>
<option value="" selected="true">No</option>
</param>
</inputs>
<outputs>
<data format="fastqsanger" name="out_file1"/>
<data format="fastqsanger" name="out_file2">
<filter>is_run['paired'] == 'yes'</filter>
</data>
</outputs>
<tests>
<test>
<param name="input1" value="fr.csf" ftype="csfasta"/>
<param name="input2" value="fr.qual" ftype="qualsolid" />
<output name="out_file1" file="f.fastq"/>
<param name="paired" value="no"/>
<param name="qual" value="0" />
<param name="trim_first_base" value="No" />
<param name="trim_name" value="No" />
<param name="double_encode" value="No"/>
</test>
<test>
<param name="input1" value="fr.csf" ftype="csfasta"/>
<param name="input2" value="fr.qual" ftype="qualsolid" />
<param name="input3" value="rr.csf" ftype="csfasta"/>
<param name="input4" value="rr.qual" ftype="qualsolid" />
<output name="out_file1" file="fr_file1.fastq"/>
<param name="paired" value="yes"/>
<param name="qual" value="0" />
<param name="trim_first_base" value="No" />
<param name="trim_name" value="Yes" />
<param name="double_encode" value="No"/>
</test>
</tests>
<help>
**What it does**
Converts output of SOLiD instrument (versions 3.5 and earlier) to fastq format suitable for bowtie, bwa, and PerM mappers.
--------
**Input datasets**
Below are examples of forward (F3) reads and quality scores:
Reads::
>1831_573_1004_F3
T00030133312212111300011021310132222
>1831_573_1567_F3
T03330322230322112131010221102122113
Quality scores::
>1831_573_1004_F3
4 29 34 34 32 32 24 24 20 17 10 34 29 20 34 13 30 34 22 24 11 28 19 17 34 17 24 17 25 34 7 24 14 12 22
>1831_573_1567_F3
8 26 31 31 16 22 30 31 28 29 22 30 30 31 32 23 30 28 28 31 19 32 30 32 19 8 32 10 13 6 32 10 6 16 11
**Mate pairs**
If your data is from a mate-paired run, you will have one additional read and quality datasets that will look similar to the ones above with one exception: the names of reads will be ending with &quot;_R3&quot;.
In this case choose **Yes** from the *Is this a mate-pair run?* drop down and you will be able to select R reads. When processing mate pairs this tool generated two output files: one for F3 reads and the other for R3 reads.
The reads are guaranteed to be paired -- mated reads will be in the same position in F3 and R3 fastq file. However, because pairing is verified it may take a while to process an entire SOLiD runs (several hours).
------
**Explanation of parameters**
**Remove reads containing color qualities below this value** - any read that contains as least one color call with quality lower than the specified value **will not** be reported.
**Trim trailing &quot;_F3&quot; and &quot;_R3&quot;?** - does just that. Not necessary for bowtie. Required for BWA.
**Trim first base?** - SOLiD reads contain an adapter base such as the first T in this read::
>1831_573_1004_F3
T00030133312212111300011021310132222
this option removes this base leaving only color calls. Not necessary for bowtie. Required for BWA.
**Double encode?** - converts color calls (0123.) to pseudo-nucleotides (ACGTN). Not necessary for bowtie. Required for BWA.
------
**Examples of output**
When all parameters are left &quot;as-is&quot; you will get this (using reads and qualities shown above)::
@1831_573_1004
T00030133312212111300011021310132222
+
%>CCAA9952+C>5C.?C79,=42C292:C(9/-7
@1831_573_1004
T03330322230322112131010221102122113
+
);@@17?@=>7??@A8?==@4A?A4)A+.'A+'1,
Setting *Trim first base from reads* to **Yes** will produce this::
@1831_573_1004
00030133312212111300011021310132222
+
%>CCAA9952+C>5C.?C79,=42C292:C(9/-7
@1831_573_1004
03330322230322112131010221102122113
+
);@@17?@=>7??@A8?==@4A?A4)A+.'A+'1,
Finally, setting *Double encode* to **Yes** will yield::
@1831_573_1004
TAAATACTTTCGGCGCCCTAAACCAGCTCACTGGGG
+
%>CCAA9952+C>5C.?C79,=42C292:C(9/-7
@1831_573_1004
TATTTATGGGTATGGCCGCTCACAGGCCAGCGGCCT
+
);@@17?@=>7??@A8?==@4A?A4)A+.'A+'1,
</help>
</tool>