Mercurial > repos > drosofff > msp_oases
changeset 4:dc684e37f668 draft
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/msp_oases commit 47a8a06faad1aa3c36503934ee2fd9e6031d3b75
| author | drosofff |
|---|---|
| date | Thu, 12 May 2016 10:45:03 -0400 |
| parents | f08d9c814ffb |
| children | 68ed36eed6c5 |
| files | oases_optimiser.py oases_optimiser.xml |
| diffstat | 2 files changed, 102 insertions(+), 113 deletions(-) [+] |
line wrap: on
line diff
--- a/oases_optimiser.py Tue Sep 29 06:29:42 2015 -0400 +++ b/oases_optimiser.py Thu May 12 10:45:03 2016 -0400 @@ -6,50 +6,48 @@ 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 * +import os, sys +import subprocess -def stop_err( msg ): - sys.stderr.write( "%s\n" % msg ) +def stop_err(msg): + sys.stderr.write("%s\n" % msg) sys.exit() -def oases_optimiser(starthash, endhash, input, job_dir): + +def oases_optimiser(starthash, endhash, input): ''' Replaces call to oases_optimiser.sh. For all k-mers between starthash and endhash run velvet and oases. ''' for i in xrange(starthash, endhash, 2): - cmd1="velveth {0}/outputFolder_{1} {1} {2} && ".format(job_dir, i, input) - cmd2="velvetg {0}/outputFolder_{1} -read_trkg yes && ".format(job_dir, i) - cmd3="oases {0}/outputFolder_{1}".format(job_dir, i) - proc = subprocess.call( args=cmd1+cmd2+cmd3, shell=True, stdout=sys.stdout, stderr=sys.stderr ) - cmd4="velveth {0}/MergedAssemblyFolder 27 -long outputFolder_*/transcripts.fa && ".format(job_dir) - cmd5="velvetg {0}/MergedAssemblyFolder -read_trkg yes -conserveLong yes && ".format(job_dir) - cmd6="oases {0}/MergedAssemblyFolder -merge yes".format(job_dir) - proc = subprocess.call( args=cmd4+cmd5+cmd6, shell=True, stdout=sys.stdout, stderr=sys.stderr ) + cmd1 = "velveth outputFolder_{0} {0} {1} && ".format(i, input) + cmd2 = "velvetg outputFolder_{0} -read_trkg yes && ".format(i) + cmd3 = "oases outputFolder_{0}".format(i) + proc = subprocess.call(args=cmd1 + cmd2 + cmd3, shell=True, stdout=sys.stdout, stderr=sys.stdout) + if not proc == 0: + print("Oases failed at k-mer %s, skipping" % i) + continue + cmd4 = "velveth MergedAssemblyFolder 27 -long outputFolder_*/transcripts.fa && " + cmd5 = "velvetg MergedAssemblyFolder -read_trkg yes -conserveLong yes && " + cmd6 = "oases MergedAssemblyFolder -merge yes" + proc = subprocess.call(args=cmd4 + cmd5 + cmd6, shell=True, stdout=sys.stdout, stderr=sys.stdout) + if not proc == 0: + raise Exception("Oases could not merge assembly") def __main__(): - job_dir= os.getcwd() - #Parse Command Line starthash = int(sys.argv[1]) endhash = int(sys.argv[2]) input = sys.argv[3] transcripts = sys.argv[4] - transcripts_path = '' - print >> sys.stdout, "PATH = %s" % (os.environ['PATH']) try: - oases_optimiser(starthash, endhash, input, job_dir) + oases_optimiser(starthash, endhash, input) except Exception, e: - stop_err( 'Error running oases_optimiser.py' + str( e ) ) - out = open(transcripts,'w') - transcript_path = os.path.join(job_dir, "MergedAssemblyFolder", 'transcripts.fa') - print >> sys.stdout, transcript_path - for line in open(transcript_path): - out.write( "%s" % (line) ) - out.close() - + stop_err('Error running oases_optimiser.py\n' + str(e)) + with open(transcripts, 'w') as out: + transcript_path = os.path.join("MergedAssemblyFolder", 'transcripts.fa') + for line in open(transcript_path): + out.write("%s" % (line)) + + if __name__ == "__main__": __main__()
--- a/oases_optimiser.xml Tue Sep 29 06:29:42 2015 -0400 +++ b/oases_optimiser.xml Thu May 12 10:45:03 2016 -0400 @@ -1,64 +1,57 @@ -<tool id="oasesoptimiserv" name="Oases_optimiser" version="1.1.4"> - <description>Auto optimise de novo RNA-seq Oases/Velvet assembly</description> - <requirements> - <requirement type="package" version="0.2.08_7a32460a60929b510037952ae815bb6e29b68123">oases</requirement> - </requirements> - <command interpreter="python"> - oases_optimiser.py '$start_hash_length' '$end_hash_length' - '#for $i in $inputs - ${i.input} - #end for - ' - '$transcripts' - </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="11" selected="yes">11</option> - <option value="13">13</option> - <option value="15">15</option> - <option value="17">17</option> - <option value="19">19</option> - <option value="21">21</option> - <option value="23">23</option> - <option value="25">25</option> - <option value="35">35</option> - <option value="45">45</option> - <option value="55">55</option> - <option value="65">65</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. Value has to be higher than the Start Hash Length value"> - <option value="25">25</option> - <option value="27">27</option> - <option value="29">29</option> - <option value="31">31</option> - <option value="33">33</option> - <option value="35">35</option> - <option value="45">45</option> - <option value="55">55</option> - <option value="63">63</option> - <option value="69">69</option> - </param> - - - - <repeat name="inputs" title="Input Files"> - <param name="input" type="data" format="fasta" label="Dataset: short reads, fasta format"/> - </repeat> - </inputs> - - <outputs> +<tool id="oasesoptimiserv" name="Oases_optimiser" version="1.1.5"> + <description>Auto optimise de novo RNA-seq Oases/Velvet assembly</description> + <requirements> + <requirement type="package" version="0.2.08_7a32460a60929b510037952ae815bb6e29b68123">oases</requirement> + </requirements> + <command><![CDATA[ + python "$__tool_directory__"/oases_optimiser.py "$start_hash_length" "$end_hash_length" + #for $i in $inputs + "${i.input}" + #end for + "$transcripts" + ]]></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="11" selected="yes">11</option> + <option value="13">13</option> + <option value="15">15</option> + <option value="17">17</option> + <option value="19">19</option> + <option value="21">21</option> + <option value="23">23</option> + <option value="25">25</option> + <option value="35">35</option> + <option value="45">45</option> + <option value="55">55</option> + <option value="65">65</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. Value has to be higher than the Start Hash Length value"> + <option value="25">25</option> + <option value="27">27</option> + <option value="29">29</option> + <option value="31">31</option> + <option value="33">33</option> + <option value="35">35</option> + <option value="45">45</option> + <option value="55">55</option> + <option value="63">63</option> + <option value="69">69</option> + </param> + <repeat name="inputs" title="Input Files"> + <param name="input" type="data" format="fasta" label="Dataset: short reads, fasta format"/> + </repeat> + </inputs> + <outputs> <data format="fasta" name="transcripts" label="${tool.name} on ${on_string}: Denovo assembled transcripts"/> - </outputs> - <tests> - <test> - <param name="input" value="input.fa" ftype="fasta"/> - <param name="start_hash_length" value="15" /> - <param name="end_hash_length" value="35" /> - <output name="transcripts" file="output.fa"/> - </test> - </tests> - + </outputs> + <tests> + <test> + <param name="input" value="input.fa" ftype="fasta"/> + <param name="start_hash_length" value="15" /> + <param name="end_hash_length" value="35" /> + <output name="transcripts" file="output.fa"/> + </test> + </tests> <help> **Oases Optimiser Overview** @@ -85,7 +78,7 @@ Read the Oases `documentation`__ for details on using the Oases Assembler. -.. _Velvet: http://www.ebi.ac.uk/~zerbino/oases/ +.. _Oases: http://www.ebi.ac.uk/~zerbino/oases/ .. __: http://www.ebi.ac.uk/~zerbino/oases/Manual.txt @@ -120,28 +113,21 @@ 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] - + -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** @@ -157,6 +143,7 @@ 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. +Occasionally oases crashes at specific k-mers, in that case this Galaxy wrapper will skip this k-mer and continue the analysis without this k-mer. **Velvetg options** @@ -175,5 +162,9 @@ short (default) - </help> + </help> + <citations> + <citation type="doi">10.1093/bioinformatics/bts094</citation> + <citation type="doi">10.1101/gr.074492.107</citation> + </citations> </tool>