Skip to content

Genome Assembly

Jump to bottom
Zane edited this page Oct 25, 2026 · 30 revisions

Genome Assembly - October 27, 2026

Biological Background:

  • Organism: Dogwood Anthracnose (Discula destructiva)
  • Host: Dogwoods (Cornus spp.)
  • Lifestyle: Hemibiotrophic pathogen
  • Genome size: 49 Mb
  • Ploidy: Haploid
image image image image
  • Project Set-Up
# Directory
cd /lustre/isaac24/proj/UTK0505/$USER
mkdir -p module3_genome_assembly
cd /lustre/isaac24/proj/UTK0505/$USER/module3_genome_assembly
  • Directory structure:
.
├── analysis
│   ├── 01_fastplong
│   ├── 02_hifiasm
│   ├── 03_assembly_qc
│   └── 04_oatk
└── raw_data
    ├── 23ddas111_1327065_S3HiC_R1.fastq.gz
    ├── 23ddas111_1327065_S3HiC_R2.fastq.gz
    └── dd_as111_100X.fastq.gz
  • Installing the necessary conda environments. Typically, I install nearly all of my software in separate environments, so feel free to do that if you like!
# General tools used in this lab.
conda create -n assembly fastp fastplong hifiasm oatk samtools gfatools bbmap

# Environment for compleasm because it clashes with oatk.
conda create -n compleasm compleasm
0. Read QC

0. Read QC

  1. Nanopore (ONT) and PacBio Instruments (HiFi) typically give really nice QC reports out-of-the-box. Let's take a peek at some!

ONT

HiFi

01. Read Trimming (fastplong)

01. Read Trimming (fastplong)

  1. Create a conda environment.
conda create -n assembly fastp fastplong hifiasm oatk samtools gfatools bbmap
  1. Link your files.
ln -s ../../raw_data/dd_as111_100X.fastq.gz .
  1. Trim your Hi-C (short) reads.
#!/bin/bash
#SBATCH --job-name=fastp
#SBATCH --nodes=1
#SBATCH --ntasks=16
#SBATCH --mem=32G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH -o outfiles/%x_%A_%a.out
#SBATCH -e outfiles/%x_%A_%a.err
#SBATCH --time=3:00:00

# activate conda environment
eval "$(conda shell.bash hook)"
conda activate fastp

# Read 1 and 2
READ_1=23ddas111_1327065_S3HiC_R1.fastq.gz
READ_2=$(basename $READ_1 | sed 's/_R1.fastq.gz/_R2.fastq.gz/')
PREFIX=$(basename $READ_1 | sed 's/_R1.fastq.gz//')

# designate trimmed reads files
TRIMMED_READ_1=$(basename $READ_1 | sed 's/_R1.fastq.gz/_trimmed_R1.fastq.gz/')
TRIMMED_READ_2=$(basename $READ_2 | sed 's/_R2.fastq.gz/_trimmed_R2.fastq.gz/')

# designate unpaired read file names
UNPAIRED_READ_1=$(basename $READ_1 | sed 's/_R1.fastq.gz/_unpaired_R1.fastq.gz/')
UNPAIRED_READ_2=$(basename $READ_2 | sed 's/_R2.fastq.gz/_unpaired_R2.fastq.gz/')

# Echo in- and outfiles in the outfiles
echo "Read 1: $READ_1"
echo "Read 2: $READ_2"
echo "Trimmed Read 1: $TRIMMED_READ_1"
echo "Trimmed Read 2: $TRIMMED_READ_2"

# Run fastp
fastp -i ${READ_1} -o ${TRIMMED_READ_1} \
        -I ${READ_2} -O ${TRIMMED_READ_2} \
        --detect_adapter_for_pe \
        --trim_poly_g \
        --thread 16 \
        --json ${PREFIX}.json \
        --html ${PREFIX}.html
  1. Trim your long reads.
#!/bin/bash
#SBATCH --job-name=fastplong
#SBATCH --nodes=1
#SBATCH --ntasks=16
#SBATCH --mem=60G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=3:00:00
#SBATCH -o outfiles/%x_%A.out
#SBATCH -e outfiles/%x_%A.err
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH --mail-user=zsmith10@vols.utk.edux

eval "$(conda shell.bash hook)"
conda activate fastplong
module load samtools

# Define variables
FASTQ=dd_as111_100X.fastq.gz
TRIMMED_FASTQ=$(basename $FASTQ | sed 's/.fastq.gz/.trimmed.fastq/')
FAILED_FASTQ=$(basename $FASTQ | sed 's/.fastq.gz/.FAILED.fastq/')
PREFIX=$(basename $FASTQ | sed 's/.fastq.gz//')

# Run fastplong
fastplong \
        --in $FASTQ \
        --out $TRIMMED_FASTQ \
        --disable_quality_filtering \
        --thread 16 \
        --json ${PREFIX}.fastplong-report.json \
        --html ${PREFIX}.fastplong-report.html

bgzip $TRIMMED_FASTQ
  1. Run multiqc!
multiqc .
02. Assembly (hifiasm)

02. Assembly (hifiasm)

  1. Link your trimmed reads.
ln -s ../01_fastplong/*trimmed*fastq.gz .
  1. Run hifiasm to assemble your genome.
#!/bin/bash
#SBATCH --job-name=hifiasm
#SBATCH --nodes=1
#SBATCH --ntasks=40
#SBATCH --mem=120G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=3:00:00
#SBATCH -o outfiles/%x_%A.out
#SBATCH -e outfiles/%x_%A.err
#SBATCH --time=3:00:00

eval "$(conda shell.bash hook)"
conda activate hifiasm

HIFI_READS=dd_as111_100X.trimmed.fastq.gz
HIC_R1=23ddas111_1327065_S3HiC_trimmed_R1.fastq.gz
HIC_R2=23ddas111_1327065_S3HiC_trimmed_R2.fastq.gz
PREFIX=dd_as111_assembly

hifiasm \
        -o $PREFIX \
        -t 40 \
        --hg-size 49m  \
        --n-hap 1 \
        --n-weight 10 \
        --f-perturb 100000 \
        --n-perturb 0.25 \
        --h1 $HIC_R1 \
        --h2 $HIC_R2 \
        $HIFI_READS \
        >& hifiasm_hgsize49m_nhap1_nweight10_fperturb10k_nperturb_0.25.log

# -o $outfile : assembly prefix
# -t 32 : threads
# --hg-size 49m  : genome size (usually not necessary)
# --n-hap 1 : number of haplotypes
# --n-weight 10 : Number of rounds of Hi-C reweighting (helps with phasing)
# --f-perturb 100000 : Number of rounds to "flip" when checking for misjoins
# --n-perturb 0.25 : Percent of the genome to "flip" when checking for misjoin during each perturbation
# --h1 $hic_r1 : Hi-C read 1
# --h2 $hic_r2 : Hi-C read 2
  1. Look at your assembly k-mer spectrum. This will give you an idea if hifiasm is finding the expected number of peaks. Usually, you should have 1 peak per ploidy. Discula destructiva is haploid, and we see one peak as expected.
[M::ha_analyze_count] lowest: count[12] = 541
[M::ha_analyze_count] highest: count[99] = 1015657
[M::ha_hist_line]     2: **************************** 286208
[M::ha_hist_line]     3: **** 36871
[M::ha_hist_line]     4: * 12640
[M::ha_hist_line]     5: * 7587
[M::ha_hist_line]     6:  4157
[M::ha_hist_line]     7:  2608
[M::ha_hist_line]     8:  2079
[M::ha_hist_line]     9:  1378
[M::ha_hist_line]    10:  791
[M::ha_hist_line]    11:  759
[M::ha_hist_line]    12:  541
[M::ha_hist_line]    13:  637
[M::ha_hist_line]    14:  430
[M::ha_hist_line]    15:  229
[M::ha_hist_line]    16:  164
[M::ha_hist_line]    17:  134
[M::ha_hist_line]    18:  156
[M::ha_hist_line]    19:  48
[M::ha_hist_line]    20:  246
[M::ha_hist_line]    21:  450
[M::ha_hist_line]    22:  364
[M::ha_hist_line]    23:  871
[M::ha_hist_line]    24:  1225
[M::ha_hist_line]    25:  3704
[M::ha_hist_line]    26:  3397
[M::ha_hist_line]    27:  2468
[M::ha_hist_line]    28:  4187
[M::ha_hist_line]    29: * 6392
[M::ha_hist_line]    30: * 6563
[M::ha_hist_line]    31: * 7401
[M::ha_hist_line]    32: * 7865
[M::ha_hist_line]    33: * 7299
[M::ha_hist_line]    34: * 5851
[M::ha_hist_line]    35: * 6536
[M::ha_hist_line]    36: * 7448
[M::ha_hist_line]    37: * 7477
[M::ha_hist_line]    38: * 7914
[M::ha_hist_line]    39: * 7306
[M::ha_hist_line]    40: * 8569
[M::ha_hist_line]    41: * 9447
[M::ha_hist_line]    42: * 9413
[M::ha_hist_line]    43: * 9606
[M::ha_hist_line]    44: * 12675
[M::ha_hist_line]    45: ** 17733
[M::ha_hist_line]    46: ** 16578
[M::ha_hist_line]    47: ** 17371
[M::ha_hist_line]    48: ** 20574
[M::ha_hist_line]    49: ** 24436
[M::ha_hist_line]    50: *** 28769
[M::ha_hist_line]    51: *** 28035
[M::ha_hist_line]    52: *** 33200
[M::ha_hist_line]    53: **** 36466
[M::ha_hist_line]    54: **** 37490
[M::ha_hist_line]    55: **** 37508
[M::ha_hist_line]    56: **** 38128
[M::ha_hist_line]    57: **** 38512
[M::ha_hist_line]    58: **** 41539
[M::ha_hist_line]    59: **** 44336
[M::ha_hist_line]    60: **** 45464
[M::ha_hist_line]    61: ***** 54460
[M::ha_hist_line]    62: ****** 58069
[M::ha_hist_line]    63: ****** 63647
[M::ha_hist_line]    64: ******* 70231
[M::ha_hist_line]    65: ******** 77575
[M::ha_hist_line]    66: ******** 83252
[M::ha_hist_line]    67: ********* 90008
[M::ha_hist_line]    68: ********** 98279
[M::ha_hist_line]    69: ************ 117064
[M::ha_hist_line]    70: ************** 139608
[M::ha_hist_line]    71: *************** 153385
[M::ha_hist_line]    72: **************** 167350
[M::ha_hist_line]    73: ******************* 189184
[M::ha_hist_line]    74: ******************** 203788
[M::ha_hist_line]    75: ********************** 222691
[M::ha_hist_line]    76: ************************* 248917
[M::ha_hist_line]    77: **************************** 283751
[M::ha_hist_line]    78: ******************************* 313645
[M::ha_hist_line]    79: *********************************** 351158
[M::ha_hist_line]    80: ************************************* 373163
[M::ha_hist_line]    81: **************************************** 403863
[M::ha_hist_line]    82: ****************************************** 428870
[M::ha_hist_line]    83: ********************************************** 470303
[M::ha_hist_line]    84: ************************************************* 493654
[M::ha_hist_line]    85: ****************************************************** 552549
[M::ha_hist_line]    86: ********************************************************** 589057
[M::ha_hist_line]    87: ************************************************************** 634094
[M::ha_hist_line]    88: ******************************************************************** 694233
[M::ha_hist_line]    89: ************************************************************************* 736953
[M::ha_hist_line]    90: ****************************************************************************** 795498
[M::ha_hist_line]    91: *********************************************************************************** 842439
[M::ha_hist_line]    92: ************************************************************************************* 868240
[M::ha_hist_line]    93: *************************************************************************************** 887811
[M::ha_hist_line]    94: **************************************************************************************** 895923
[M::ha_hist_line]    95: ******************************************************************************************* 928766
[M::ha_hist_line]    96: *********************************************************************************************** 962281
[M::ha_hist_line]    97: ************************************************************************************************ 976522
[M::ha_hist_line]    98: *************************************************************************************************** 1001864
[M::ha_hist_line]    99: **************************************************************************************************** 1015657
[M::ha_hist_line]   100: *************************************************************************************************** 1008436
[M::ha_hist_line]   101: ************************************************************************************************** 992095
[M::ha_hist_line]   102: ************************************************************************************************ 977306
[M::ha_hist_line]   103: *********************************************************************************************** 963001
[M::ha_hist_line]   104: ********************************************************************************************** 953401
[M::ha_hist_line]   105: ******************************************************************************************** 935992
[M::ha_hist_line]   106: *************************************************************************************** 884211
[M::ha_hist_line]   107: ********************************************************************************** 830597
[M::ha_hist_line]   108: ****************************************************************************** 787771
[M::ha_hist_line]   109: ************************************************************************* 737663
[M::ha_hist_line]   110: ********************************************************************* 698739
[M::ha_hist_line]   111: ***************************************************************** 655752
[M::ha_hist_line]   112: ************************************************************* 621422
[M::ha_hist_line]   113: ********************************************************** 589590
[M::ha_hist_line]   114: ***************************************************** 539818
[M::ha_hist_line]   115: ************************************************* 493792
[M::ha_hist_line]   116: ******************************************** 441823
[M::ha_hist_line]   117: ************************************** 387669
[M::ha_hist_line]   118: ********************************** 349436
[M::ha_hist_line]   119: ******************************* 319712
[M::ha_hist_line]   120: *************************** 277357
[M::ha_hist_line]   121: *********************** 236866
[M::ha_hist_line]   122: ******************** 202427
[M::ha_hist_line]   123: ****************** 178630
[M::ha_hist_line]   124: **************** 158024
[M::ha_hist_line]   125: ************** 137863
[M::ha_hist_line]   126: *********** 114951
[M::ha_hist_line]   127: ********** 99818
[M::ha_hist_line]   128: ******** 85601
[M::ha_hist_line]   129: ******* 75671
[M::ha_hist_line]   130: ****** 62108
[M::ha_hist_line]   131: ***** 51300
[M::ha_hist_line]   132: **** 42754
[M::ha_hist_line]   133: *** 34086
[M::ha_hist_line]   134: *** 27961
[M::ha_hist_line]   135: ** 23675
[M::ha_hist_line]   136: ** 18259
[M::ha_hist_line]   137: * 14572
[M::ha_hist_line]   138: * 13556
[M::ha_hist_line]   139: * 13054
[M::ha_hist_line]   140: * 9927
[M::ha_hist_line]   141: * 8842
[M::ha_hist_line]   142: * 6514
[M::ha_hist_line]  rest: ******************* 191062
03. Assembly QC: (bbmap & compleasm)

03. Assembly QC: BBMap & Compleasm

  1. First, we need a new environment because compleasm and oatk clash.
conda create -n compleasm compleasm
  1. Link your assembly!
ln -s ../02_hifiasm/*
  1. Run bbmap to check basic assembly statistics. Note: I ran this as a for-loop, but there is only one assembly we need to assess in this case.
#!/bin/bash
#SBATCH --job-name=bbmap_stats
#SBATCH --nodes=1
#SBATCH --ntasks=2
#SBATCH --mem=8G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=00:30:00
#SBATCH -o outfiles/%x_%A.out
#SBATCH -e outfiles/%x_%A.err
#SBATCH --mail-type=BEGIN,END,FAIL
#SBATCH --mail-user=zsmith10@vols.utk.edux

eval "$(conda shell.bash hook)"
conda activate bbmap

for ASM in *p_ctg.fasta
do
        STATS=$(basename $ASM | sed 's/.fasta/.stats.txt/')

        stats.sh \
                -Xmx8g \
                in=$ASM \
                > $STATS
done
  1. Run compleasm to check BUSCO completeness.
#!/bin/bash
#SBATCH --job-name=compleasm
#SBATCH --nodes=1
#SBATCH --ntasks=12
#SBATCH --mem=50G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=3:00:00
#SBATCH -o outfiles/%x_%A.out
#SBATCH -e outfiles/%x_%A.err

eval "$(conda shell.bash hook)"
conda activate compleasm

ASM="dd_as111_assembly.hic.p_ctg.fasta"
LINEAGE="sordariomycetes_odb12"
OUTDIR=$(basename $ASM | sed 's/.fasta/.COMPLEASM/')

echo "$ASM"
echo "$OUTDIR"
echo "$LINEAGE"


# Run compleasm -- make sure you have version 0.2.8
compleasm --version

compleasm \
        run \
        -a $ASM \
        -o $OUTDIR \
        -t 12 \
        -l $LINEAGE
04. Organelle Assembly (oatk)

  1. Here, we'll link our trimmed HiFi reads back over.
ln -s ../01_fastplong/*trimmed.fastq.gz .
  1. Next, let's make sure we have the Oatk databases.
git clone https://github.com/c-zhou/OatkDB.git
  1. Now, we can assemble our mitochondrion! Plant plastids can also be assembled with oatk. Typically, you would expect your mitochondrial read depth to be approximately ~5-10x higher than your nuclear reads, while plastid reads would be expect to be ~50x-100x higher.

Here, we have 100x nuclear genome coverage, so we would expect ~500x-1000x mitochondrial read coverage. Using a different assembler (like Spades), you can expect similar read ratios in Illumina WGS reads.

These are useful rules of thumb for the -c parameter, which is minimum read coverage of a unitig to be considered a candidate for organelle assembly. However, it's best to be conservative and underestimate.

#!/bin/bash
#SBATCH --job-name=oatk
#SBATCH --nodes=1
#SBATCH --ntasks=20
#SBATCH --mem=100G
#SBATCH -A ACF-UTK0032
#SBATCH --partition=short
#SBATCH --qos=short
#SBATCH --time=3:00:00
#SBATCH -o outfiles/%x_%A.out
#SBATCH -e outfiles/%x_%A.err

eval "$(conda shell.bash hook)"
conda activate oatk

# HiFi reads
HIFI_READS="dd_as111_100X.trimmed.fastq.gz"

# Set output prefix: 'dd_as111_100X'
PREFIX="${HIFI_READS%%.*}"

# Specify your database
MITO_DB="${PWD}/OatkDB/v20230921/dikarya_mito.fam"

# Run oatk
oatk \
    -k 1001 \
    -c 60 \
    -t 20 \
    --nhmmscan nhmmscan \
    -m "${MITO_DB}" \
    -o "${PREFIX}_oatk_default" \
    "$HIFI_READS"
  1. Now, we can visually inspect our mitochondrial genome using the Bandage software, which you can download here: https://github.com/rrwick/Bandage

To view in Bandage:

  • Secure copy your GFA to your desktop
  • Load bandage
  • Load your GFA: File => Load Graph => Draw Graph
# Insert your NetID below.
scp '<YOUR_NETID>@dtn1.isaac.utk.edu:/lustre/isaac24/proj/UTK0505/<YOUR_NETID>/module3_genome_assembly/analysis/04_oatk/*mito.gfa' .
image

Assembly Scaffolding - October 29, 2026

Hi-C Contact Map for Discula destructiva

image

Genome Annotation - November 3, 2026

Clone this wiki locally