cleanup SPAdes and expand readme

README.md

@@ -22,7 +22,7 @@ optional arguments:
                         [union]
   --kmer KMER, -k KMER  Kmer to use, depends on read length [81]
   --min_contig_len MIN_CONTIG_LEN, -l MIN_CONTIG_LEN
-                        Minimum contig length in final assembly [600]
+                        Minimum contig length in final assembly [1000]
   --min_kmers_threshold MIN_KMERS_THRESHOLD, -m MIN_KMERS_THRESHOLD
                         Exclude k-mers occurring less than this [25]
   --max_kmers_threshold MAX_KMERS_THRESHOLD, -x MAX_KMERS_THRESHOLD
@@ -35,6 +35,35 @@ optional arguments:
   --version             show program's version number and exit
 ```
 
+##Input parameters
+The following parameters change the results:
+
+__action__: There are two fundamental methods of operation. The default is 'union', where kmers which occur in ANY trait sample, but are absent from the nontrait samples, get used to filter the reads. So in effect you are assembling the whole accessory genome of the trait samples. This leads to larger end assemblies and more false positives, but will capture greater regions of the accessory genome. It is tolerant to situations where you have a plasmid which can vary substantially with different backbones or payloads.  The next is 'intersection', where kmers must occur in ALL trait samples and not in the nontrait samples. This leads to smaller end assemblies and more fragmentation, with less false positives.  It is less tolerant to variation. 
+
+__kmer__: Choosing a kmer is not an exact science, and have greatly influence the final results. This kmer size is used by KMC for counting and filtering, and by SPAdes for assembly.  Ideally it should be between about 60-90% of the read length, should be an odd number and between 21 and 127 (SPAdes restriction).  If choose a kmer which is too small, you will get a lot more false positives. If you choose a kmer too big, you will use a lot more RAM and potentially get too little data returned. Quite often with Illumina data the beginning and end of the reads have higher sequencing error rates. Ideally you want a kmer size which sits nicely inside the good cycles of the read. Trimming with Trimmomatic can help if the quality collapses quite badly at the end of the read.
+
+__min_contig_len__: This needs to be about a minimum of twice the mean fragment size (insert size) of your library to reduce the impact of false positives. For example if you have a single kmer which randomly occurs in the genome, using it will then allow for reads upstream and downstream, plus their mates, to be assembled. This variable can control this noise. Setting this too high will lead to valuable information being lost (e.g. small plasmids).
+
+__min_kmers_threshold__: This value lets you set a minimum threshold for the occurance of a kmer. Ideally you need about 30X depth of coverage to perform de novo assembly. This value default to 25, so excluding kmers below this level eliminates kmers where you wont get a good assembly, thus reducing false positives. The maximum value is 254, but the results poor.
+
+__max_kmers_threshold__: This value lets you set a maximum threshold for the occurance of a kmer. The occurance of kmers forms a Poisson distribution, with a very long tail. With KMC, there is a catchall bin for occurances of 255 and greater (so 255 is the maximum value). By default it is set to 254 which excludes this catchall bin for kmers, and thus the long tail of very common kmers. This reduces the false positives. You need to be careful when setting this lower because you could exclude all of the interesting kmers.
+
+
+The following parameters have no impact on the results:
+
+__threads__: This sets the number of threads available to KMC and SPAdes. It should never be more than the number of CPUs available on the server. If you use a compute cluster, make sure to request the same number of threads on a single server. It defaults to 1 and you will get a reasonable speed increase by adding a few CPUs, but the benefit tails off quite rapidly since the I/O becomes the limiting factor (speed of reading files from a disk or network).
+
+__spades_exec__: By default SPAdes is assumed to be in your PATH and called spades.py. You can set this to point to a different executable, which might be required if you have multiple versions of SPAdes installed.
+
+__verbose__: By default the output is limited and all intermediate files are deleted. Setting this flag allows you output more details of the software as it runs and it keeps the intermediate files.
+
+##Required resources
+###RAM (memory)
+The largest consumer of RAM (memory) is SPAdes. Assembling a whole bacteria takes approximately 4GB of RAM. If the filtering allows everything through then this worst case will occur, but generally less than 1GB of RAM is required. Poor quality sequencing data will increase the amount of RAM required. In this instance running Trimmomatic first will help greatly.
+
+###Disk space
+By default all of the intermediate files are cleaned up at the end, so the overall disk space usage is quite low. As an example, an input of 800 Mbytes of compressed reads created 40 Mbytes of output data at the end. While the algorithm is running the disk usage will never exceed the size of the input reads. The intermediate files can be kept if you use the 'verbose' option. 
+
 #Outputs 
 For every trait sample you will get an assembly of nucleotide sequences in FASTA format. You will also get a text file describing the process, with versions of software, parameters used and references.
 

VERSION

@@ -1 +1 @@
-0.0.2
+0.0.3

plasmidtron/PlasmidTron.py

@@ -50,9 +50,11 @@ def run(self):
 			kmc_filter.filter_fastq_file_against_kmers()
 			kmc_filters.append(kmc_filter)
 
+		spades_assemblies = []
 		for sample in trait_samples:
 			spades_assembly = SpadesAssembly( sample, self.output_directory, self.threads, self.kmer, self.spades_exec, self.min_contig_len)
 			spades_assembly.run()
+			spades_assemblies.append(spades_assembly)
 			print(spades_assembly.filtered_spades_assembly_file() + '\n')
 			sample.cleanup()
 
@@ -69,4 +71,7 @@ def run(self):
 
 			for kmc_filter in kmc_filters:
 				kmc_filter.cleanup()
+
+			for spades_assembly in spades_assemblies:
+				spades_assembly.cleanup()
 

plasmidtron/SpadesAssembly.py

@@ -2,6 +2,7 @@
 import logging
 import tempfile
 import subprocess
+import shutil
 from Bio import SeqIO
 
 '''Assemble a filtered sample with SPAdes'''
@@ -20,10 +21,10 @@ def spades_command(self):
 		return ' '.join([self.spades_exec, '--careful', '--only-assembler','-k', str(self.kmer), '-1', self.sample.filtered_forward_file, '-2', self.sample.filtered_reverse_file, '-o', self.spades_output_directory])
 
 	def spades_assembly_file(self):
-		return os.path.join(self.spades_output_directory,'contigs.fasta')
+		return os.path.join(self.spades_output_directory,'scaffolds.fasta')
 
 	def filtered_spades_assembly_file(self):
-		return os.path.join(self.spades_output_directory,'filtered_contigs.fasta')
+		return os.path.join(self.spades_output_directory,'filtered_scaffolds.fasta')
 
 	def remove_small_contigs(self,input_file, output_file):
 		with open(input_file, "r") as spades_input_file, open(output_file, "w") as spades_output_file:
@@ -38,3 +39,10 @@ def run(self):
 		self.logger.info("Assembling sample" )
 		subprocess.call(self.spades_command(), shell=True)
 		self.remove_small_contigs(self.spades_assembly_file(), self.filtered_spades_assembly_file())
+
+	def cleanup():
+		shutil.rmtree(os.path.join(self.spades_output_directory, 'tmp' ))
+		shutil.rmtree(os.path.join(self.spades_output_directory, 'mismatch_corrector' ))
+		shutil.rmtree(os.path.join(self.spades_output_directory, 'misc' ))
+		shutil.rmtree(os.path.join(self.spades_output_directory, 'K'+self.kmer ))
+

plasmidtron/tests/PlasmidTron_test.py

@@ -31,9 +31,9 @@ def test_small_valid_chrom_plasmid(self):
 		plasmid_tron = PlasmidTron(options)
 		plasmid_tron.run()
 
-		final_assembly = os.path.join(data_dir,'out/spades_S_typhi_CT18_chromosome_pHCM2/filtered_contigs.fasta')
+		final_assembly = os.path.join(data_dir,'out/spades_S_typhi_CT18_chromosome_pHCM2/filtered_scaffolds.fasta')
 
-		self.assertTrue(os.path.isfile(os.path.join(data_dir,'out/spades_S_typhi_CT18_chromosome_pHCM2/contigs.fasta')))
+		self.assertTrue(os.path.isfile(os.path.join(data_dir,'out/spades_S_typhi_CT18_chromosome_pHCM2/scaffolds.fasta')))
 		self.assertTrue(os.path.isfile(final_assembly))
 		'''The final assembly should be about 6k so leave some margin for variation in SPAdes'''
 		self.assertTrue(os.path.getsize(final_assembly) > 5000)

scripts/plasmidtron

@@ -23,7 +23,7 @@ parser.add_argument('file_of_trait_fastqs', help='File of filenames of trait (ca
 parser.add_argument('file_of_nontrait_fastqs', help='File of filenames of nontrait (control) FASTQs', type=InputTypes.is_file_of_nontrait_fastqs_valid)
 parser.add_argument('--action',	 '-a', help='Control how the traits kmers are filtered for assembly [%(default)s]', choices=['intersection','union'],  default = 'union')
 parser.add_argument('--kmer',	 '-k', help='Kmer to use, depends on read length [%(default)s]', type=InputTypes.is_kmer_valid,  default = 81)
-parser.add_argument('--min_contig_len',	 '-l', help='Minimum contig length in final assembly [%(default)s]', type=InputTypes.is_min_contig_len_valid, default = 600)
+parser.add_argument('--min_contig_len',	 '-l', help='Minimum contig length in final assembly [%(default)s]', type=InputTypes.is_min_contig_len_valid, default = 1000)
 parser.add_argument('--min_kmers_threshold',	 '-m', help='Exclude k-mers occurring less than this [%(default)s]', type=InputTypes.is_min_kmers_threshold_valid,  default = 25)
 parser.add_argument('--max_kmers_threshold',	 '-x', help='Exclude k-mers occurring more than this [%(default)s]', type=InputTypes.is_max_kmers_threshold_valid,  default = 254)
 parser.add_argument('--threads',  '-t', help='Number of threads [%(default)s]', type=InputTypes.is_threads_valid,  default = 1)