Assembly With Canu

Author:Brant C. Faircloth
Copyright:This documentation is available under a Creative Commons (CC-BY) license.

Modification History

See Assembly With Canu.

Purpose

There are several options to assemble PacBio long-read data, and one of those (potentially the easier to install) is canu (another is to the the SMRTAnalysis pipeline and/or pb-assembly ). canu works reasonably well on @QB2 - I’ve just generally learned that it’s easier to run in single-threaded mode rather than try to make the grid mode work (it seems as if grid mode most likely will NOT work on the queueing system that we use.

Steps

  1. Because canu is compute intensive, the following steps have been documented assuming you are using @QB2

  2. Compile canu according to Compiling Canu

  3. Create a pacbio environment for conda (after installing miniconda and configuring for bioconda)

    conda create -n pacbio python=2.7 bam2fastx
    
  4. Create a working directory for your data (here, I’m just using the _Arabidopsis_ test data):

    mkdir arabidopsis-pacbio && cd arabidopsis-pacbio
    
  5. Download Arabidopsis test data from PacBio. Be sure to get the *.pbi files because we need them to convert the bam data to fastq format

    wget -P pacbio-raw https://downloads.pacbcloud.com/public/SequelData/ArabidopsisDemoData/SequenceData/1_A01_customer/m54113_160913_184949.subreads.bam
    wget -P pacbio-raw https://downloads.pacbcloud.com/public/SequelData/ArabidopsisDemoData/SequenceData/1_A01_customer/m54113_160913_184949.subreads.bam.pbi
    
    wget -P pacbio-raw https://downloads.pacbcloud.com/public/SequelData/ArabidopsisDemoData/SequenceData/3_C01_customer/m54113_160914_092411.subreads.bam
    wget -P pacbio-raw https://downloads.pacbcloud.com/public/SequelData/ArabidopsisDemoData/SequenceData/3_C01_customer/m54113_160914_092411.subreads.bam.pbi
    
  6. canu requires data in fastq format, so convert each bam file to fastq.

    #!/bin/bash
    #PBS -q single
    #PBS -A <allocation>
    #PBS -l walltime=06:00:00
    #PBS -l nodes=1:ppn=2
    #PBS -V
    #PBS -N bam_to_fastq
    #PBS -o bam_to_fastq.out
    #PBS -e bam_to_fastq.err
    
    # load the parallel module to run files in parallel (up to 4 cores in single queue)
    module load gnuparallel/20170122
    
    # activate our conda env
    source activate pacbio
    
    
    cd $PBS_O_WORKDIR
    mkdir pacbio-fastq && cd pacbio-fastq
    find ../pacbio-raw/ -name *.bam | parallel "bam2fastq -o {/.} {}"
    

    Note

    You may need to adjust queues and cores to suit your needs. Here, I’m using the single queue because I only have 2 files to convert and we can use up to 4 CPUs in single. Also note that you may need to adjust the time needed for each run - particularly for larger bam files you are converting.

  7. Once those data are converted, we can kick off the canu assembly job. Again, I’ve found that we need to keep these assembly jobs “local”, meaning that we’re not going to run in grid mode. However, you do want to run them using the bigmem queue on @QB2. Also note here that we’re redirecting stdout and stderr to files - we’re doing this so that we can check on job status as the runs go along (since the queuing system typically keeps these in temp files until the end of the run):

    #!/bin/bash
    #PBS -q bigmem
    #PBS -A <allocation>
    #PBS -l walltime=72:00:00
    #PBS -l nodes=1:ppn=48
    #PBS -V
    #PBS -N canu_config
    #PBS -o canu_config.out
    #PBS -e canu_config.err
    #PBS -m abe
    #PBS -M [email protected]
    
    module load gcc/6.4.0
    module load java/1.8.0
    
    cd $PBS_O_WORKDIR
    mkdir -p canu-assembly && cd canu-assembly
    
    canu \
        -p arabidopsis \
        -d arabidopsis-pacbio \
        genomeSize=123m \
        useGrid=false \
        -pacbio-raw ../pacbio-fastq/*.fastq.gz 1>canu-assembly.stdout 2>canu-assembly.stderr
    

    Warning

    You will need to adjust the genome size of your organism in the above to something that’s suitable. Gb-size genome size is set using genomeSize=1.1g, which would be appropriate for a bird.

  8. If you need to restart the job at any time (e.g., you run out of walltime, which is likely), you may want to rename canu-assembly.stdout and canu-assembly.stderr, so they are not overwritten:

    for i in canu-assembly/*.std*; do mv $i $i.old; done
    
  9. Then you simply need to resubmit the qsub script and the job will restart from where it last started.