# HG changeset patch # User drosofff # Date 1463064303 14400 # Node ID dc684e37f66873103a3820c1758911c32444515f # Parent f08d9c814ffb378dea6dd897a4e5517dd6842a07 planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/msp_oases commit 47a8a06faad1aa3c36503934ee2fd9e6031d3b75 diff -r f08d9c814ffb -r dc684e37f668 oases_optimiser.py --- 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__() diff -r f08d9c814ffb -r dc684e37f668 oases_optimiser.xml --- 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 @@ - - Auto optimise de novo RNA-seq Oases/Velvet assembly - - oases - - - oases_optimiser.py '$start_hash_length' '$end_hash_length' - '#for $i in $inputs - ${i.input} - #end for - ' - '$transcripts' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + Auto optimise de novo RNA-seq Oases/Velvet assembly + + oases + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - + + + + + + + + + **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) - + + + 10.1093/bioinformatics/bts094 + 10.1101/gr.074492.107 +