Long Read Correction and other Correction tools

This package is a loose collection of scripts. To run the correction routine see the section below. Descriptions of the other scripts are at the bottom of this file.

Contact: gurtowsk@cshl.edu

Correction

In short, the correction algorithm takes as input the unitigs from a short read assembly and uses them to correct long read data. More background information for the algorithm can be found: http://schatzlab.cshl.edu/presentations/2013-06-18.PBUserMeeting.pdf

Running Correction.

ECTOOLS_HOME=/path/to/this/directory

1. Create Celera Unitigs from Illumina (or other high identity) data

   Output: organism.utg.fasta

2. Create a new directory to house the correction

   $> mkdir organism_correct

3. Simlink your organism.utg.fasta into the correction dir

   $> cd organism_correct $> ln -s /path/to/organism.utg.fasta

4. Generally we only correct reads greater than 1kb. Filter out shorter reads if possible. If you have it, 20-30x of raw long read coverage is preferred so that you have ~20x remaining after correction.

5. Partition pacbio reads into small batches

   $> python ${ECTOOLS_HOME}/partition.py 20 500 pbreads.legnth_filtered.fa

   Output: Series of directories and a ReadIndex.txt which maps reads to partitions

   *Note: You probably only want a few reads per file (20 is usually good)

6. Copy correction script to working dir

   $> cp ${ECTOOLS_HOME}/correct.sh .

7. Modify the global variables a the top of correct.sh

8. Ensure Nucmer suite is in your PATH

   $> nucmer

9. Launch Correction The correct.sh script should be run in each of the partition directories. The number of files in directory should be specified as the array parameter to grid engine. This simple for loop works well:

   $> for i in {0001..000N}; do cd $i; qsub -cwd -j y -t 1:${NUM_FILES_PER_PARTITION} ../correct.sh; cd ..; done

   Where N is the number of partitions (directories) and NUM_FILES_PER_PARTITION is the number of files per partition specified in the partitions script (500 in this case)

10. Wait for correction to complete. Concatenate all corrected output files from all parititons into a single fa file:

    $> cat ????/*.cor.fa > organism.cor.fa

11. Run Celera on the output file.

    *Notes: Use convert-fasta-to-v2.pl to make celera frg file from organism.cor.fa

    example: $> convert-fasta-to-v2.pl -l organism_pbcor -s organism.cor.fa \ -q <( python ${ECTOOLS_HOME}/qualgen.py organism.cor.fa) \

    organism.cor.frg

    Celera's default parameters work pretty well with corrected PB data, you can also drop the utgGraphErrorRate to 0.01.

    Make sure you use OverlapBasedTrimming (it is on by default).

12. (Optional) If SGE is not available in your compute environment, you can use GNU parallel to run ECTools. Although not the recommended method, as alignment can be very compute intensive, for small genomes (bacteria), this method can be used on a single machine.

For each of the directories created by the partition script (0001..000N), cd into the directory and run:

$>for j in {1..500}; do echo "SGE_TASK_ID=$j TMPDIR=/tmp ../correct.sh"; done | \ parallel -j <# of compute cores>

Script Descriptions

@@ schtats @@ Calculates read length statistics for sequence files.

@@ simulate_reads.py @@ Simple read simulator. Takes a genome and a file of line separated read lengths. A read is produced for each given read length at the specified error rate.

@@ gccontent.py @@ Builds an ascii graph of gc content for each item in a fasta file

@@ io.py @@ Has various parsers for parsing delta files

@@ join.py @@ general purpose script to do database like join on flat file.

@@ group_by.py @@ general purpose script for grouping sorted line oriented data by a column.

@@ filter.py @@ general purpose script for filtering line separated data. Takes a database and uses it to filter the query file.

@@ gc_count.py @@ summarizes gc content from a set of alignments

@@ fix_delta_header.py @@ Because the correct.sh writes all output to a temporary directory in SGE and the copies it back to the original working directory, delta files will have the header of the temporary directory. This script just fixes the header so that downstream mummer tools work after the pipeline has completed.