Mercurial > repos > konradpaszkiewicz > oasesoptimiser

changeset 0:d35f713f5dc4 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 1.0.0 from old tool shed archive to new tool shed repository
author konradpaszkiewicz
date Tue, 07 Jun 2011 17:09:27 -0400
parents
children
files README oases_optimiser.py oases_optimiser.sh oases_optimiser.xml
diffstat 4 files changed, 403 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README	Tue Jun 07 17:09:27 2011 -0400
@@ -0,0 +1,15 @@
+#Created 07/01/2011
+#Konrad Paszkiewicz, Exeter Sequencing Service, University of Exeter 
+
+Oases optimiser
+
+This Galaxy tool is intended to provide a crude but effective method of producing semi-optimised transcriptome assemblies. It operates by performing a number of user defined Velvet assemblies upon which Oases is then run. The results of all of these runs are then put through a final Oases step where they act as scaffolds for the main assembly.
+
+Prerequisites:
+
+1. Enclosed scripts
+2. Working installation of Velvet and Oases
+
+Limitations:
+
+This is a crude optimisation step which DOES NOT try to evaluate the quality of each assembly or map to any reference cDNA. 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/oases_optimiser.py	Tue Jun 07 17:09:27 2011 -0400
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+
+"""
+VelvetOptimiser Wrapper
+Adapted from velveth and velvetg tools in Galaxy
+Konrad Paszkiewicz	University of Exeter, UK.
+
+"""
+import pkg_resources;
+import logging, os, string, sys, tempfile, glob, shutil, types, urllib
+import shlex, subprocess
+from optparse import OptionParser, OptionGroup
+from stat import *
+
+
+log = logging.getLogger( __name__ )
+
+assert sys.version_info[:2] >= ( 2, 4 )
+
+def stop_err( msg ):
+    sys.stderr.write( "%s\n" % msg )
+    sys.exit()
+
+def __main__():
+    #Parse Command Line
+    s = 'oases_optimiser.py:  argv = %s\n' % (sys.argv)
+    #print >> sys.stderr, s # so will appear as blurb for file
+    argcnt = len(sys.argv)
+    starthash = sys.argv[1]
+    endhash = sys.argv[2]
+    inputs = sys.argv[3]
+    contigs = sys.argv[4]
+    LastGraph = sys.argv[5]
+    afgFile = sys.argv[6]
+    unused_reads_fasta = sys.argv[7]
+    stats = sys.argv[8]
+    othervelvetgoptions = sys.argv[9]
+    otheroasesoptions = sys.argv[10]
+    transcripts = sys.argv[11]
+    splicingevents = sys.argv[12]
+    contigordering = sys.argv[13]
+    #contigs.extra_files_path
+
+    
+    working_dir = '' 
+    
+    cmdline = '/users/galaxy/galaxyscripts/oases_optimiser.sh %s %s \'%s\' %s %s  2&1>/dev/null' % (starthash, endhash, inputs, othervelvetgoptions, otheroasesoptions)
+    #print >> sys.stderr, cmdline # so will appear as blurb for file
+    try:
+        proc = subprocess.Popen( args=cmdline, shell=True, stderr=subprocess.PIPE )
+        returncode = proc.wait()
+        # get stderr, allowing for case where it's very large
+        stderr = ''
+        buffsize = 1048576
+        try:
+            while True:
+                stderr += proc.stderr.read( buffsize )
+                if not stderr or len( stderr ) % buffsize != 0:
+                    break
+        except OverflowError:
+            pass
+        if returncode != 0:
+            raise Exception, stderr
+    except Exception, e:
+        stop_err( 'Error running oases_optimiser.py' + str( e ) )
+    out = open(transcripts,'w')
+    transcript_path = os.path.join(working_dir,'transcripts.fa')
+    #print >> sys.stderr, path
+    for line in open(transcript_path):
+        out.write( "%s" % (line) )
+    out.close()
+ 
+    out = open(splicingevents,'w')
+    path = os.path.join(working_dir,'splicing_events.txt')
+    #print >> sys.stderr, contigs_path
+    for line in open(path ):
+        out.write( "%s" % (line) )
+    out.close()
+
+    
+
+    out = open(contigs,'w')
+    contigs_path = os.path.join(working_dir,'contigs.fa')
+    #print >> sys.stderr, contigs_path
+    for line in open(contigs_path ):
+        out.write( "%s" % (line) )
+    out.close()
+    out = open(stats,'w')
+    stats_path = os.path.join(working_dir,'stats.txt')
+    for line in open( stats_path ):
+        out.write( "%s" % (line) )
+    out.close()
+    if LastGraph != 'None':
+        out = open(LastGraph,'w')
+        LastGraph_path = os.path.join(working_dir,'LastGraph')
+        for line in open( LastGraph_path ):
+            out.write( "%s" % (line) )
+        out.close()
+    if afgFile != 'None':
+        out = open(afgFile,'w')
+        afgFile_path = os.path.join(working_dir,'oases_asm.afg')
+        try:
+            for line in open( afgFile_path ):
+                out.write( "%s" % (line) )
+        except:
+            logging.warn( 'error reading %s' %(afgFile_path))
+            pass
+        out.close()
+    if unused_reads_fasta != 'None':
+        out = open(unused_reads_fasta,'w')
+      	unused_reads_fasta_path = os.path.join(working_dir,'UnusedReads.fa')
+        try:
+            for line in open( unused_reads_fasta_path ):
+                out.write( "%s" % (line) )
+        except:
+            logging.info( 'error reading %s' %(unused_reads_fasta_path))
+            pass
+        out.close()
+
+if __name__ == "__main__": __main__()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/oases_optimiser.sh	Tue Jun 07 17:09:27 2011 -0400
@@ -0,0 +1,26 @@
+#!/bin/bash
+
+min_k=$1
+max_k=$2
+input=$3
+args=$#
+if [$args -gt 3]; then
+	velvetoptions=$4
+	oasesoptions=$5
+else
+	velvetoptions=''
+	oasesoptions=''
+fi
+
+for i in `seq $min_k 2 $max_k`
+do
+     velveth . $i $input
+     velvetg . -read_trkg yes 
+     oases . $oasesoptions
+     mv transcripts.fa transcripts_$i.fa 
+done
+
+velveth . $max_k -long transcripts_*.fa $input 
+velvetg . -read_trkg yes $velvetoptions 
+oases . -amos_file yes -unused_reads yes $oasesoptions 
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/oases_optimiser.xml	Tue Jun 07 17:09:27 2011 -0400
@@ -0,0 +1,242 @@
+<tool id="oasesoptimiser" name="oasesoptimiser" version="1.0.0">
+	<description>Auto optimise a denovo RNA-seq Oases/Velvet assembly</description>
+	<command interpreter="python">
+	oases_optimiser.py  '$start_hash_length' '$end_hash_length' 
+	'#for $i in $inputs
+		${i.file_format}
+		${i.read_type}
+		${i.input}
+           #end for
+	'
+	'$contigs'
+	'$LastGraph'
+	'$velvet_asm'
+	'$unused_reads_fasta'
+	'$stats'
+	'$othervelvetgoptions'
+	'$otheroasesoptions'
+	'$transcripts'
+	'$splicingevents'
+	'$contigordering'
+	'$contigs.extra_files_path'
+	</command>
+
+        <inputs>
+          <param label="Start Hash Length" name="start_hash_length" type="select" help="k-mer length in base pairs of the words being hashed. Shorter hash lengths (i.e. less than 31) may cause out-of-memory problems.">
+
+                        <option value="23">23</option>
+                        <option value="25">25</option>
+                        <option value="27">27</option>
+                        <option value="29">29</option>
+			<option value="31" selected="yes">31</option>
+                        <option value="33">33</option>
+                        <option value="35">35</option>
+                        <option value="37">37</option>
+                        <option value="39">39</option>
+                        <option value="41">41</option>
+
+          </param>
+	 <param label="End Hash Length" name="end_hash_length" type="select" help="k-mer length in base pairs of the words being hashed.">
+                        <option value="43">43</option>
+                        <option value="45">45</option>
+                        <option value="47">47</option>
+                        <option value="49">49</option>
+                        <option value="51">51</option>
+                        <option value="53">53</option>
+                        <option value="55">55</option>
+                           <option value="57">57</option>
+                         <option value="59">59</option>
+                        <option value="61">61</option>
+                        <option value="63">63</option>
+                        <option value="65">65</option>
+                        <option value="67">67</option>
+			 <option value="69" selected="yes">69</option>
+
+	 </param>
+
+
+	<param label="Advanced velvetg options" name="othervelvetgoptions" title="Advanced velvetg options" help="Other Velvetg options - see below for details" type="text" size="30">
+		
+	</param>
+
+	
+        <param label="Advanced Oases options" default="-ins_length 50 -ins_length_sd 100" name="otheroasesoptions" title="Advanced oases options" help="Other oases options - see below for details" type="text" size="30">
+
+        </param>
+
+    <param name="strand_specific" type="boolean" checked="false" truevalue="-strand_specific" falsevalue="" label="Use strand specific transcriptome sequencing" help="If you are using a strand specific transcriptome sequencing protocol, you may wish to use this option for better results."/> 
+          <repeat name="inputs" title="Input Files">
+                      <param label="file format" name="file_format" type="select">
+                              <option value="-fasta" selected="yes">fasta</option>
+                              <option value="-fastq">fastq</option>
+                              <option value="-eland">eland</option>
+                              <option value="-gerald">gerald</option>
+                      </param>
+                      <param label="read type" name="read_type" type="select">
+                              <option value="-short" selected="yes">short reads</option>
+                              <option value="-shortPaired">shortPaired reads</option>
+                              <option value="-short2">short2 reads</option>
+                              <option value="-shortPaired2">shortPaired2 reads</option>
+                              <option value="-long">long reads (reads >200bp e.g. 454 reads)</option>
+                              <option value="-longPaired">longPaired reads</option>
+                      </param>
+
+
+            <param name="input" type="data" format="fasta,fastq,eland,gerald" label="Dataset"/>
+          </repeat>
+        </inputs>
+	
+	<outputs>
+
+	        <data format="fasta" name="transcripts" label="${tool.name} on ${on_string}: Denovo assembled transcripts"/>
+		<data format="tabular" name="splicingevents" label="${tool.name} on ${on_string}: Splicing events in assembled transcripts"/>
+              <data format="tabular" name="contigordering" label="${tool.name} on ${on_string}: Contig ordering in assembled transcripts"/> 
+	       <data format="fasta" name="contigs" label="${tool.name} on ${on_string}: Denovo assembled contigs"/>
+
+                <data format="fasta" name="unused_reads_fasta" label="${tool.name} on ${on_string}: Unused Reads">
+                 <!--  <filter>unused_reads['generate_unused'] == "yes"</filter>-->
+                </data>
+                <data format="tabular" name="stats" label="${tool.name} on ${on_string}: Stats"/>
+   		  <data format="afg" name="velvet_asm" label="${tool.name} on ${on_string}: AMOS.afg">
+                  <!-- <filter>generate_amos['afg'] == "yes"</filter> -->
+                </data>
+
+  		  <data format="txt" name="LastGraph" label="${tool.name} on ${on_string}: LastGraph">
+		    <!-- <filter>last_graph['generate_graph'] == "yes"</filter> -->
+                </data>
+
+            
+	</outputs>
+	<requirements>
+		<requirement type="package">velvet</requirement>
+	        <requirement type="package">oases</requirement>
+
+	</requirements>
+	
+	<help>
+**Oases Optimiser Overview**
+
+Oases_ is a de novo transcriptome assembler specially designed to work with short reads. It is an extension of Velvet developed by Daniel Zerbino and Ewan Birney at the European Bioinformatics Institute (EMBL-EBI), near Cambridge, in the United Kingdom.
+
+Provide all the information about insert lengths and their standard deviation as
+possible (identical to Velvet). In the advanced oases options you should set the -ins_length flag when using paired-end reads (e.g. -ins_length 200 -ins_length_sd 100). If you do not do this, then paired-end information will not be used by Oases!
+
+----------------------------------------------------------------------------------
+
+Oases produces a number of output files, which correspond to the different algorithms
+being run succesively on the data. In the above example, you would find:
+
+transcripts.fa
+	A FASTA file containing the transcripts imputed directly from trivial
+	clusters of contigs (loci with less than two transcripts and Confidence Values = 1)
+	and the highly expressed transcripts imputed by dynamic
+	programming (loci with more than 2 transcripts and Confidence Values less than 1).
+
+splicing_events.txt
+	A hybrid file which describes the contigs contained in each locus in FASTA
+	format, interlaced with one line descriptions of splicing events using the
+	AStalavista nomenclature*.
+
+Read the Oases `documentation`__ for details on using the Oases Assembler.
+
+.. _Velvet: http://www.ebi.ac.uk/~zerbino/oases/
+
+.. __: http://www.ebi.ac.uk/~zerbino/oases/Manual.txt
+
+------
+
+**Other outputs (as per Velvet)**
+
+
+**Contigs**
+
+The *contigs.fa* file.
+This fasta file contains the sequences of the contigs longer than 2k, where k is the word-length used in velveth. If you have specified a min contig lgth threshold, then the contigs shorter than that value are omitted.
+Note that the length and coverage information provided in the header of each contig should therefore be understood in k-mers and in k-mer coverage (cf. 5.1) respectively.
+The N's in the sequence correspond to gaps between scaffolded contigs. The number of N's corresponds to the estimated length of the gap. For reasons of compatibility with the archives, any gap shorter than 10bp is represented by a sequence of 10 N's.
+
+**Stats**
+
+The *stats.txt* file.
+This file is a simple tabbed-delimited description of the nodes. The column names are pretty much self-explanatory. Note however that node lengths are given in k-mers. To obtain the length in nucleotides of each node you simply need to add k - 1, where k is the word-length used in velveth.
+The in and out columns correspond to the number of arcs on the 5' and 3' ends of the contig respectively.
+The coverages in columns short1 cov, short1 Ocov, short2 cov, and short2 Ocov are provided in k-mer coverage (5.1).
+Also, the difference between # cov and # Ocov is the way these values are computed. In the first count, slightly divergent sequences are added to the coverage tally. However, in the second, stricter count, only the sequences which map perfectly onto the consensus sequence are taken into account.
+
+**LastGraph**
+
+The *LastGraph* file.
+This file describes in its entirety the graph produced by Velvet.
+
+**AMOS.afg**
+The *velvet_asm.afg* file.
+This file is mainly designed to be read by the open-source AMOS genome assembly package. Nonetheless, a number of programs are available to transform this kind of file into other assembly file formats (namely ACE, TIGR, Arachne and Celera). See http://amos.sourceforge.net/ for more information.
+The file describes all the contigs contained in the contigs.fa file (cf 4.2.1).
+
+**Advanced options**
+	    -ins_length integer	   : expected distance between two paired-end reads in the first short-read dataset (default: no read pairing)
+	    -ins_length2 integer          : expected distance between two paired-end reads in the second short-read dataset (default: no read pairing)
+        -ins_length_long integer      : expected distance between two long paired-end reads (default: no read pairing)
+        -ins_length*_sd  integer      : est. standard deviation of respective dataset (default: 10% of corresponding length)
+        -scaffolding  yes|no             : scaffolding of contigs used paired end information (default: on)
+        -max_branch_length  integer      : maximum length in base pair of bubble (default: 100)
+        -max_divergence  floating-point  : maximum divergence rate between two branches in a bubble (default: 0.2)
+        -max_gap_count  integer          : maximum number of gaps allowed in the alignment of the two branches of a bubble (default: 3)
+        -min_pair_count  integer         : minimum number of paired end connections to justify the scaffolding of two long contigs (default: 10)
+        -max_coverage  floating point    : removal of high coverage nodes AFTER tour bus (default: no removal)
+        -long_mult_cutoff  int           : minimum number of long reads required to merge contigs (default: 2)
+	-min_trans_length
+	simple threshold on output transcript length
+	-cov_cutoff 
+	minimum number of times a k-mer has to be observed to be used in the 
+	assembly (just like in Velvet) [default=3]
+	-min_pair_cov
+	minimum number of times two contigs must be connected by reads or read pairs to be clustered together [default=4]
+	-paired_cutoff 
+	minimum ratio between the numbers of observed and expected connecting
+	read pairs between two contigs [default=0.1]
+
+
+
+**Hash Length**
+
+The hash length, also known as k-mer length, corresponds to the length, in base pairs, of the words being hashed. 
+
+The hash length is the length of the k-mers being entered in the hash table. Firstly, you must observe three technical constraints::
+
+# it must be an odd number, to avoid palindromes. If you put in an even number, Velvet will just decrement it and proceed.
+# it must be below or equal to MAXKMERHASH length (cf. 2.3.3, by default 31bp), because it is stored on 64 bits
+# it must be strictly inferior to read length, otherwise you simply will not observe any overlaps between reads, for obvious reasons.
+
+Now you still have quite a lot of possibilities. As is often the case, it's a trade- off between specificity and sensitivity. Longer kmers bring you more specificity (i.e. less spurious overlaps) but lowers coverage (cf. below). . . so there's a sweet spot to be found with time and experience.
+We like to think in terms of "k-mer coverage", i.e. how many times has a k-mer been seen among the reads. The relation between k-mer coverage Ck and standard (nucleotide-wise) coverage C is Ck = C # (L - k + 1)/L where k is your hash length, and L you read length.
+Experience shows that this kmer coverage should be above 10 to start getting decent results. If Ck is above 20, you might be "wasting" coverage. Experience also shows that empirical tests with different values for k are not that costly to run! VelvetOptimiser automates these tests for you.
+
+
+**Velvetg options**
+
+
+
+**Input Files**
+
+Velvet works mainly with fasta and fastq formats. For paired-end reads, the assumption is that each read is next to its mate
+read. In other words, if the reads are indexed from 0, then reads 0 and 1 are paired, 2 and 3, 4 and 5, etc.
+
+Supported file formats are::
+
+  fasta
+  fastq 
+  eland
+  gerald
+
+Read categories are::
+
+  short (default)
+  shortPaired
+  short2 (same as short, but for a separate insert-size library - i.e. if you have two libraries of different lengths)
+  shortPaired2 (see above)
+  long (for Sanger, 454 or even reference sequences)
+  longPaired
+
+	</help>
+</tool>