de novo genome assembly from PacBio HiFi (High Fidelity) long reads

tags: assembly, purge_dups, BUSCO, hifiasm, small

This pipeline was developed for mammal genomes and includes some steps that are necessary for assembling these larger (~3Gb) genomes but that may not be needed for smaller genomes.

• Table of Contents
  
  • de novo genome assembly from PacBio HiFi (High Fidelity) long reads
    • Step 1: Run Hifiasm
    • Step 2: Re-run HiFiasm with different parameters and compare the quality of the assembly
      • 2a. Evaluate assembly with Busco
      • 2b. Evaluate assembly with assembly_stats
    • Step 3: Run purge_dups to improve the assembly
      • 3a. Generate an alignment and evaluate read coverage
      • 3b. Visualize the coverage histogram to evaluate the default cutoff values
      • 3c. Set your cutoffs manually
      • 3d. Split your assembly
      • 3e. Now you can purge!
    • Step 4: Evaluate your assemblies
    • Step 5: Get rid of contamination
      • 5a. Mask repeats with RepeatMasker

Step 1: Run Hifiasm

The first step is to use the HiFiasm software to assemble a genome from the raw PacBio reads. (You'll start with the slurm part at the top, as always). Then, the general structure for assembly of heterozygous genomes is:

hifiasm -o prefix -t numberThreads input.fasta.gz --write-paf --write-ec /dev/null

so the whole script looks something like this

#!/bin/bash
#SBATCH -N 1
#SBATCH -n 48
#SBATCH -t 24:00:00
#SBATCH -o assembly_ORCAI-homcov19s25_022724.out 
#SBATCH -e assembly_ORCAI-homcov19s25_022724.err

#I installed HiFiasm via anaconda
source sackettl/miniconda3/etc/profile.d/conda.sh
source activate hifiasm_env 

cd /sackettl/squirrels/ORCA_GMGS/

hifiasm -o ORCAI-homcov19s25v2 -t 44 purge_dups-1.2.6/m64086e_230206_201241.hifi_reads.fasta.gz purge_dups-1.2.6/m64086e_230210_154702.hifi_reads.fasta.gz purge_dups-1.2.6/m64086e_230212_010601.hifi_reads.fasta.gz --hom-cov 19 -s 0.25 --write-paf --write-ec /dev/null

#get primary contigs in FASTA
awk '/^S/{print ">"$2;print $3}' ORCAI-homcov19s25v2.bp.p_ctg.gfa > ORCAI-homcov19s25v2.p_ctg.fa

The last line is to translate the HiFiasm results to a fasta assembly, which is what you will want for most downstream processes.

---

Step 2: Re-run HiFiasm with different parameters and compare the quality of the assembly

You want to end up with an assembly that is close to the expected genome size and has the largest scaffold N50 and smallest number of scaffolds that you can achieve. Different combinations of parameters (e.g., s, D) can help you optimize the assembly quality.

You should also check the quality of each assembly with busco and assembly_stats and compare their outputs.

2a. Evaluate assembly with Busco

First, evaluate with busco using the correct database for your taxon:

#!/bin/bash
#SBATCH -N 1
#SBATCH -n 48
#SBATCH -t 12:00:00
#SBATCH -o busco_ORCAI-homcov19s25D10_041823.out 
#SBATCH -e busco_ORCAI-homcov19s25D10_041823.err

# I ran busco from its own conda environment
source /sackettl/miniconda3/etc/profile.d/conda.sh
conda activate busco5_env 
#source activate busco_env

#depending on how busco installs, you may need to use 'python run_BUSCO.py' and then the command
busco -i /sackettl/squirrels/ORCA_homcov19s25D10/ORCAI-homcov19s25D10.p_ctg.fa \
        -l mammalia_odb10 -c 44 -m genome -o ORCAI-homcov19s25D10_busco -f

2b. Evaluate assembly with assembly_stats

Next, evaluate your assembly using assembly_stats, which you can also install using anaconda:

#!/bin/bash
#SBATCH -N 1
#SBATCH -n 20
#SBATCH -t 1:30:00
#SBATCH -o assemblystats_ORCAI-homcov19s25D10_041823.out 
#SBATCH -e assemblystats_ORCAI-homcov1925D10_041823.err

source /sackettl/miniconda3/etc/profile.d/conda.sh
conda activate assemblystats_env

assembly-stats /sackettl/squirrels/ORCA_homcov19s25D10/ORCAI-homcov19s25D10.p_ctg.fa > ORCAI-homcov19s25D10_assembly-stats.out

Based on the busco scores and assembly stats for each variation of the assembly, you can choose the best one for use in purge_dups.

---

Step 3: Run purge_dups to improve the assembly

This mini-pipeline is based off the suggested pipeline for purge_dups and starts with the genome generated by HiFiasm above.

3a. Generate an alignment and evaluate read coverage

First, we need to map the raw PacBio HiFi reads to the draft genome from HiFiasm. We will use minimap2 with the map-hifi setting. However, map-hifi is not available as part of the purge_dups conda package, so we will need to install it separately.

git clone https://github.com/lh3/minimap2
cd minimap2 && make

Now we can run the mapping step:

#SBATCH lines

export JOBS_PER_NODE=44

cd /sackettl/ORCA_GMGS/

for i in ./*fasta.gz; do ./software/minimap2/minimap2 -I 10G -x map-hifi -t 46 ./ORCA-homcov19s25l2.p_ctg.fa $i > $i.paf; done

This will generate .paf alignment files for each reads file and use those as input for the next step, which will both generate reads statistics and estimate the cutoffs to use when purging:

#SBATCH lines

export JOBS_PER_NODE=44
source /project/sackettl/miniconda3/etc/profile.d/conda.sh
conda activate purgedups126_ENV

cd /sackettl/ORCA_GMGS/

pbcstat ./*.paf
calcuts PB.stat > automatic.cutoffs 2>calcuts_auto.log

The pbcstat step will generate PB.base.cov and PB.stat files. You can use these same input files to test different levels of purging, since they are specific to your initial assembly+reads. (but rename and generate new files if you try different mapping parameters.)

3b. Visualize the coverage histogram to evaluate the default cutoff values

Now we need to look at the histogram of read coverage using the alignment of our reads to our draft genome. It is possible that the automatically generated cutoffs are fine, but it is a good idea to try a few different values and evaluate them.

First, let purge_dups generate the automatic cutoff values:

#SBATCH lines

export JOBS_PER_NODE=44

source /miniconda3/etc/profile.d/conda.sh
conda activate purgedups126_ENV

cd /sackettl/ORCA_GMGS/

calcuts PB.stat > ORCA-homcov19s25l2_hifi_auto.cutoffs 2>calcuts_auto_hifi.log

Now, take a look at the histogram:

#SBATCH lines

export JOBS_PER_NODE=44

source /miniconda3/etc/profile.d/conda.sh
conda activate purgedups126_ENV

cd /sackettl/ORCA_GMGS/

pip install matplotlib

python3 ../purge_dups-1.2.6/scripts/hist_plot.py -c ORCA-homcov19s25l2_manual_l2m13u59.cutoffs PB.asm20.stat PB.cov.asm20_cut2-13-59.png 

Now download this to your computer with scp [e-mail]:/location/of/png ./local/folder/ and take a look at it. It should look something like this

There is some guidance on how to set the thresholds for your data here, and you should trust the purge_dups authors' recommendations over mine.

Essentially, you want to find the lower limit of coverage that you think represents junk (things that were sequenced at very low coverage), which is the point at the far-left side of the plot before the histogram line starts to increase. In my data, I compared results with -l1 and -l2 and ultimately decided on a cutoff of 2, although the results were nearly identical (the auto-generated suggestion for lower cutoff was 5).

For the middle threshold, you want to look for the valley if you have two peaks. After running some trials and evaluating purged assemblies with BUSCO & assembly-stats (see Step 2), I chose -m15 .

The upper regions of the coverage plot represent repeated regions, but choosing the upper threshold seems less straightforward to me. In many examples presented by users here, the suggested upper threshold was midway on the right tail. The automatic cutoffs for my assembly chose 45. I opted to start with a conservative upper limit of 60 and compare this to lower values (the lower -u values even below 45 produced assemblies with good BUSCO scores and high scaffold N50 values). You should try some different values and then check the busco scores to evaluate the number of complete, single-copy, duplicated, and missing buscos (to make sure you don't purge too much).

3c. Set your cutoffs manually

You will need to change the cutoffs to some other values to explore the results of different degrees of purging. This just requires flags for lower, middle, and upper cutoffs:

calcuts -l 2 -m 15 -u 50 PB.stat > ORCA-homcov19s25l2_manual_l2m15u50.cutoffs 2> calcuts_manual-l2m15u50.log 

Personally, although super-long filenames are annoying, I want them to include the specific parameters I used because I know I will be creating several different purged assemblies using different cutoff values. So I keep the cutoff values in all filenames, and rename the purged assemblies as soon as they are complete.

3d. Split your assembly

#SBATCH lines

export JOBS_PER_NODE=44

source /miniconda3/etc/profile.d/conda.sh
conda activate purgedups126_ENV

cd /sackettl/ORCA_GMGS/

split_fa ./ORCAI-homcov19s25l2.p_ctg.fa > ORCAI-homcov19s25l2.p_ctg.fa.split

Now you will use that split assembly in minimap2, outside of the purge_dups environment if you have HiFi data:

#SBATCH lines

export JOBS_PER_NODE=44

cd /sackettl/ORCA_GMGS/

./software/minimap2/minimap2 -x map-hifi -DP ORCAI-homcov19s25l2.p_ctg.fa.split ORCAI-homcov19s25l2.p_ctg.fa.split | gzip -c - > ORCAI-homcov19s25l2.p_ctg.fa.split.self.paf.gz

3e. Now you can purge!

Now that you have determined what the appropriate levels of purging may be (but check the results of multiple cutoff values, because you still won't know for sure!), you can purge the duplicates and junk, and keep the rest. You do so in these two steps:

purge_dups -2 -T ORCA-homcov19s25l2_manual_l2m15u50.cutoffs -c PB.cov ORCAI-homcov19s25l2.p_ctg.fa.split.self.paf.gz > ORCA-homcov19s25l2_manualcutoffs-l2m15u50_dups.bed 2> ORCA-homcov19s25l2_manualcutoffs-l2m15u50_purge_dups.log

get_seqs -e ORCA-homcov19s25l2_manualcutoffs-l2m15u50_dups.bed ./ORCAI-homcov19s25l2.p_ctg.fa

This will produce two output files: one called purge.fa, which is your purged assembly, and another called hap.fa, which is the removed scaffolds with a tag stating the reason they were removed (OVLP, REPEAT, etc.).

Rename those with your purging parameters (or move them to a folder titled with those parameters) so you can keep track:

mv purged.fa ORCA-homcov19s25l2_manual_l2m15u50_purged.fa
mv hap.fa ORCA-homcov19s25l2_manual_l2m15u50_hap.fa

Step 4: Evaluate your assemblies

You should check each assembly with assemblystats and busco as in Step 2.

My strategy was to use the assembly that had the largest scaffold N50/ smallest number of scaffolds/ smallest genome size that preserved the maximum number of BUSCOs.

---

Step 5: Get rid of contamination

Many assemblies are contaminated with bacterial sequences, so we need to get rid of these. If we do it now, it will make scaffolding and annotation easier.

5a. Mask repeats with RepeatMasker

Genome masking means hiding certain regions of the genome, such as repeats. Masking a genome improves the efficacy of the decontamination by emphasizing the regions that contribute the most to the classification (Saini et al. 2016, Caetano-Anolles).

We will both hard-mask (convert all repeat regions to Ns) and soft-mask (make repeat sequences lower-case) because we will need both versions for downstream analyses. Before proceeding, make sure the sequences in your assembly are uppercase.

conda activate repeatmasker417_env

RepeatMasker -species rodentia -s -parallel 18 -xsmall -gff ORCA-homcov19s25l2_x2cutsl5m15u50_purged.fasta

Next, you can change all lowercase letters in your assembly to Ns using sed:

sed /^>/!y/atcgn/NNNNN/ assembly_purged.softmasked.fasta > assembly_purged.hardmasked.fasta

Explanation: In sequence lines (not starting with “>”), replace each instance of lower case bases (a t c or g), including undefined ones (n), to Ns.