Evaluating dual/diploid assemblies with Flagger

Overview

Here is a description of a read-based pipeline that can detect different types of mis-assemblies in a draft dual/diploid assembly. (What is a dual assembly? Read this page). The core component of this pipeline is Flagger. Flagger recieves the read alignments to a draft diploid assembly, detects anomalies in read coverage along the assembly and partitions the assembly into 4 main components; erroneous, (falsely) duplicated, haploid (correctly assembled) and collapsed.

This evaluation has 3 steps:

• Align long reads to the diploid assembly
• Phase and relocalize the reads with secondary alignments using secphase (Optional)
• Run Flagger

1. Align long reads

The ONT and HiFi reads can be aligned to a dual/diploid assembly (~ 6Gbases long in human) with a long read aligner like winnowmap and minimap2. Here are the commands for producing the alignments (taken from the winnowmap docs):

  # making the k-mer table with meryl
  meryl count k=15 output merylDB asm.fa
  meryl print greater-than distinct=0.9998 merylDB > repetitive_k15.txt

  # alignment with winnowmap
  winnowmap -W repetitive_k15.txt -ax [map-ont | map-pb] -Y -L --eqx --cs -I8g <(cat pat_asm.fa mat_asm.fa) reads.fq.gz | \
    samtools view -hb | samtools sort -@8 > read_alignment.bam

Any other appropriate long read aligner can be employed in this step.

2. Relocalize wrongly phased reads (Optional)

In this step we use Secphase to phase and relocalize the reads with multiple alignments. To be more precise all the secondary and primary alignments of the same read are scored based on marker consistency and the alignment with the highest score is selected as the primary alignment. The output of this section is a corrected version of the input bam file, in which the primary and secondary alignments are swapped whenever neccessary. Secphase can work properly only if the secondary alignments are available with their full sequence and base quality array.

More information about Secphase is available here

3. Run Flagger

The produced alignment file (in step 1) can be used as the input to Flagger. Flagger outputs a bed file with 5 labels; erroneous (Err), duplicated (Dup), haploid (Hap), collapsed (Col) and unkown (Unk). Any component other than "haploid" is pointing to unreliable blocks in assembly and unkown label is for the bases that couldn't be assigned confidently to the other components (for example when assigned blocks are extremely short). These components are explained in detail here.

More information about Flagger is available here

Running pipeline with WDL

It is easier to run the pipeline using the WDLs described below. A WDL file can be run locally using Cromwell, which is an open-source Workflow Management System for bioinformatics. The latest releases of Cromwell are available here and the documentation is available here.

It is recommended to run the whole pipeline using flagger_end_to_end_with_mapping.wdl.

Here is a list of input parameters for flagger_end_to_end_with_mapping.wdl (The parameters marked as "(Mandatory)" are mandatory to be defined in the input json):

Using CHM13 annotation files (Optional)

A set of annotations files in the coordinates of chm13v2.0 are prepared beforehand and they can be used as inputs to the workflow. These bed files are useful when there is no denovo annotation for the assemblies and we like to project annotations from chm13v2.0 to the assembly coordinates for two main purposes:

1. Detecting regions with coverage biases: Satellite repeat arrays might have coverage biases so before running Flagger the pipeline will detect potentially baised regions. Flagger will then fit a separate Gaussian model to each detected annotation.
2. Stratifying final results with the projected annotations.

Using these bed files are optional and the related parameters can be left undefined (being absent from the input json). In this case the workflow should still work properly. However users should be aware that in this case potential coverage biases in HSat arrays may mislead the pipeline. The final summary tsv files will contain zero values for any stratification whose bed file was left undefined.

All files with git and s3 urls are publicly accessible so if you are running the WDL on Terra/AnVIL platforms you can use the same urls. They are also available in the directories misc/annotations and misc/biased_regions of this repository for those who want to run locally.

Other major workflows

flagger_end_to_end_with_mapping.wdl is running two major workflows; long_read_aligner_scattered.wdl and flagger_end_to_end.wdl. The first one runs read mapping and the second one runs the end-to-end flagger pipeline (including annotation projection, bias detection, secphase and flagger).

Users can run each of them separately. For example if there is a read alignment file available beforehand users can run only flagger_end_to_end.wdl. The parameters for each of these workflows is a subset of the parameters listed above for flagger_end_to_end_with_mapping.wdl therefore this table can still be used as a reference for either long_read_aligner_scattered.wdl or flagger_end_to_end.wdl.

Running WDLs locally using Cromwell

Below are the main commands for running flagger_end_to_end_with_mapping.wdl locally using Cromwell.

wget https://github.com/broadinstitute/cromwell/releases/download/85/cromwell-85.jar
wget https://github.com/broadinstitute/cromwell/releases/download/85/womtool-85.jar

# Get version 0.4.0 of Flagger
wget https://github.com/mobinasri/flagger/archive/refs/tags/v0.4.0.zip

unzip v0.4.0.zip

# make a directory for saving outputs and json files
mkdir workdir 

cd workdir

java -jar ../womtool-85.jar inputs ../flagger-0.4.0/wdls/workflows/flagger_end_to_end_with_mapping.wdl > inputs.json

After modifying inputs.json, setting mandatory parameters: sampleName, suffixForFlagger, suffixForMapping, hap1AssemblyFasta, hap2AssemblyFasta, readfiles, preset and removing unspecified parameters (they will be set to default values), users can run the command below:

## run flagger workflow
java -jar ../cromwell-85.jar run ../flagger-0.4.0/wdls/workflows/flagger_end_to_end_with_mapping.wdl -i inputs.json -m outputs.json

The paths to output files will be saved in outputs.json. The instructions for running any other WDL is similar.

Running WDLs on Slurm using Toil

Instructions for running WDLs on Slurm are provided here , which includes some test data sets for each of the workflows:

• long_read_aligner_scattered.wdl
• flagger_end_to_end.wdl
• flagger_end_to_end_with_mapping.wdl

Dockstore links

All WDLs are uploaded to Dockstore for easier import into platforms like Terra or AnVIL.

• Dockstore link for long_read_aligner_scattered.wdl
• Dockstore link for flagger_end_to_end.wdl
• Dockstore link for flagger_end_to_end_with_mapping.wdl

Components

Each of these components has their own color when they are shown in the IGV or the UCSC Genome Browser.

Running Flagger on HPRC assemblies

The haplotype-resolved assemblies of the HPRC-Y1 samples and their corresponding data sets are available in

https://s3-us-west-2.amazonaws.com/human-pangenomics/index.html?prefix=working/HPRC/

For more details read this github page:

https://github.com/human-pangenomics/HPP_Year1_Assemblies

We have used the Genbank version of the HPRC-Y1 assemblies.

The v0.1 results are available in https://s3-us-west-2.amazonaws.com/human-pangenomics/index.html?prefix=submissions/e9ad8022-1b30-11ec-ab04-0a13c5208311--COVERAGE_ANALYSIS_Y1_GENBANK/FLAGGER/

Publications

Liao, Wen-Wei, Mobin Asri, Jana Ebler, Daniel Doerr, Marina Haukness, Glenn Hickey, Shuangjia Lu et al. "A draft human pangenome reference." Nature 617, no. 7960 (2023): 312-324.