Indian Institute of Science

Revisiting Graph-theoretic Models for
Genome Assembly in the Era of Long Reads
Chirag Jain
Assistant Professor
Computational and Data Sciences
Indian Institute of Science

Based on Bioinformatics paper
https://doi.org/10.1093/bioinformatics/btad124
1

CiE (Computability in Europe) 2023
July 24, 2023

Indian Institute of Science

A human genome is a set of 46 strings
Two “copies”, inheritedfrom father and mother(~99.9% identical)

Sum total length isabout 6 billion
characters

Image source: NHGRI
2

Indian Institute of Science

Human-genome sequencing(IDEAL)
Text file with the 46 strings
Human DNACCGenome assembly: Reconstruction of the original genome from reads

C

TCAC ATAC

Reads

3

T

TTAC

GTCC
C

CCGT

C

(ACTUAL)

TCAC

T

CCGT

GACC

G

A sequencing
instrument

Indian Institute of Science

Latest: Long and accurate sequencing• Enables de novo genome assembly of both maternal and paternal haplotypes• Was not feasible using previous technologies
State-of-the-art assemblies:
3 Gbp (collapsed)
6 Gbp haplotype-resolved

100%

Long-readaccuracy

PacBio HiFi: Accuracy = 99.8%, Avg. length > 10 kbp
Nanopore: Accuracy = ~99%, Avg. length > 40 kbp
2020
onwards

2010-2080%

10 kbp4Long-read length (kb)Indian Institute of Science

Latest: Long and accurate sequencing
100

2000

2005

2010

2015

PB HiFi

ONT

0.01

ONT Duplex

1

PromethION

Mean contig
N50 (Mbp)

PB Sequel

ONT MinION

Thanks to long-read assemblers !

2020

Assembly submission year (NCBI database)

5

Figure from Kovaka et al. 2023, “Approaching complete genomes, transcriptomes and epi-omes with accurate long-read sequencing”

Indian Institute of Science

Graph-theoretic models for assembly• Input: Set of reads R•

De Bruijn graph : Bk(R)[Idury and Waterman 1995]• Vertices are distinct k-mers observed in R
• An edge implies a su x-pre x overlap of length k − 1 between two k-mers• Overlap graph : Ok(R)• Vertices are input reads• An edge implies an exact su x-pre x match of length ≥ k• String graph : Sk(R)[Myers 1995, 2005]• Subgraph of Ok(R)
fi

ffi

fi

• Next slide…
ffi

6Indian Institute of Science

Graph-theoretic models for assembly• Input: Set of reads R•

De Bruijn graph : Bk(R)[Idury and Waterman 1995]• Vertices are distinct k-mers observed in R
• An edge implies a su x-pre x overlap of length k − 1 between two k-mers• Overlap graph : Ok(R)• Vertices are input reads• An edge implies an exact su x-pre x match of length ≥ k between two reads• String graph : Sk(R)[Myers 1995, 2005]• Subgraph of Ok(R)
fi

ffi

fi

• Next slide…
ffi

7

Indian Institute of Science

Graph-theoretic models for assembly• Input: Set of reads R•

De Bruijn graph : Bk(R)[Idury and Waterman 1995]• Vertices are distinct k-mers observed in R
• An edge implies a su x-pre x overlap of length k − 1 between two k-mers• Overlap graph : Ok(R)• Vertices are input reads• An edge implies an exact su x-pre x match of length ≥ k between two reads• String graph : Sk(R)[Myers 1995, 2005]• Subgraph of Ok(R)
fi

ffi

fi

• Next slide…
ffi

8Indian Institute of Science

String graph• Used in most long-read assemblers• Sub-graph of overlap graph[Myers 1995, 2005]• Keep only the longest su x-pre x overlap between a pair of reads• Remove contained reads• Remove transitive edges

fi

ffi

9

ACTGCTTAC
CTGCTT
ACTGCT
CTGCTGTGCTGA

ACTGCT

CTGCTGTGCTGA

Indian Institute of Science

Graph sparsification cannot be avoided in practiceOverlap graph obtained using a subset of simulated
nanopore reads from human chr20

Removecontained reads
Removetransitive edges

String graph10

Indian Institute of Science

Are graph models “coverage-preserving”?
• Suppose input reads cover the entire genome, do we have a guaranteethat the “true” chromosomes can be spelled as a walk in the graph?

11

Indian Institute of Science

Are graph models “coverage-preserving”?
• Suppose input reads cover the entire genome, do we have a guaranteethat each candidate chromosome can be spelled as a walk in the graph?

Say R = TATACA, CATATA

A

C

T

A
A

T

Candidate 1

12

A

C

T

A

A

T
T

A

Candidate 2

Indian Institute of Science

Are graph models “coverage-preserving”?
• Suppose input reads cover the entire genome, do we have a guaranteethat each candidate chromosome can be spelled as a walk in the graph?

Say R = TATACA, CATATA

A

C

T

A
A

T

Candidate 1

13

A

C

T

A

A

T
T

A

Candidate 2• Circular string z is a candidate if ∃l1, l2 ∈ ℕ, l2 > l1 such that all intervals oflength l1 in z include the starting position of at least one read of length l2

Indian Institute of Science

Theoretical evaluation• Input: set of reads R
Graph model

Guarantee?

de Bruijn graph Bk(R)YES

Proof technique∀ k ≤ l2 − l1 + 1

By contradiction∀ k ≤ l2 − l1

Algorithm to identify aclosed walk in the graphfor each candidate

Overlap graph Ok(R)YESString graph Sk(R)NO for any k

Counter-example
Consistent with prior works[e.g., Hui et al. ISIT 2016]

Jain “Coverage-preserving sparsi cation of overlap graphs for long-read assembly” Bioinformatics 2023
fi

14

Indian Institute of Science

Theoretical evaluation• Input: set of reads R
Graph model

Guarantee?

de Bruijn graph Bk(R)YES

Proof technique∀ k ≤ l2 − l1 + 1

By contradiction∀ k ≤ l2 − l1

Algorithm to identify aclosed walk in the graphfor each candidate

Overlap graph Ok(R)YESString graph Sk(R)NO for any k

Counter-example
Consistent with prior works[e.g., Hui et al. ISIT 2016]

Jain “Coverage-preserving sparsi cation of overlap graphs for long-read assembly” Bioinformatics 2023
fi

15

Indian Institute of Science

Proof sketch
Assembly-graph model Coverage-preserving?

Proof techniquede Bruijn graph Bk(R)

By contradiction

YES ∀ k ≤ l2 − l1 + 1

interval oflength l1• Assume there is a candidate chromosome z not spelled
by graph⇒ At least one k-mer in z is absent from the set of vertices

16zread of length ≥ l2

Indian Institute of Science

Proof sketch
Assembly-graph model Coverage-preserving?

Proof technique

Overlap graph Ok(R)Proposed algorithm to identify aclosed walk for each candidate

YES ∀ k ≤ l2 − l1intervalszz

Candidate circularchromosome z

All intervals oflength l1 in z

fi

ffi

17

reads

reads

z

Pull out reads of length ≥ l2that start in the intervalsz

Identify a subset of these reads
with su x-pre x overlaps

Indian Institute of Science

Counter example for string graph
Assembly-graph model Coverage-preserving?
String graph

Say R =TGTGCA
CACGTGCACGTGGTGGGCA18

NO for any k

G
T
G

C

Candiate 1

T

G

G
A

G

C

G

C

Candiate 2

T

A
C

G

Candidate 1 cannot be spelled in graphafter contained read CACGTG is removed

Indian Institute of Science

Further questions addressed• Does it really matter for genome assembly quality in reality?
• Is there an alternate method to sparsify overlap graph that ispractical and provably-good?

• Good heuristics to recover non-redundant contained reads?

19

Indian Institute of Science‘Safe’ rules to sparsify overlap graph

Practical?

• It is safe to remove a vertex (or edge) if the set of circular string walks
remains unchanged• Transitive edge reduction in [Myers 2005] is safe
• Removing a contained read is safe if it maps to only a single candidate
genome at a unique location
parentread pcontained
read x[formal proof in paper]

candiate 2

candiate 1

candiate 320

Indian Institute of Science

Heuristic - 1• By computing all-versus-all read alignments, we can estimate if acontained read maps uniquely within a single haplotype (either paternal ormaternal)

overlapping reads(using hi asm)

contained read

C
C
C
C
C
A
A
A
A
check for heterozygous mutationsfi

21

Indian Institute of Science

Heuristic - 2• Estimate if contained read contributes a “novel” string walk in the graph• κ1 = set of k-mers observed within bounded length string walks in theassembly graph from a contained read• κ2 = set of k-mers observed similarly from its “parent” reads• Remove contained read if κ1⊆ κ2
κ1read xCTGCTT
p1 AACTGCTTACTC
p2 ACTGCTTGG

κ2

graph

x
p1
p2

22

Indian Institute of Science

ContainX
Prototype implementation in C++
Input:(a) Long reads

github.com/at-cg/ContainX

Heuristic 1

Heuristic 2

Discard contained
reads with het
mutationsDiscard contained
reads which lack novel
walks

Build overlap graph

Apply transitive-edge
reduction [Myers 2005](b) all-vs-allalignments

23

Output nonredundant containedreadsIndian Institute of Science

Benchmark datasets• Simulated error-free long reads; length distribution matches real data
• Human genomes: CHM13 (haploid), HG002 (diploid)
Data set

24

Count of reads

N50 length

Max length

HAPLOID-20x-ONT-1

3.7M

40K

570K

HAPLOID-20x-ONT-2

3.7M

40K

540K

HAPLOID-20x-HiFi-1

2.9M

21K

49K

HAPLOID-20x-HiFi-2

2.9M

21K

49K

DIPLOID-30x-ONT-1

5.3M

40K

540K

DIPLOID-30x-ONT-2

5.3M

40K

570K

DIPLOID-30x-HiFi-1

4.2M

21K

49K

DIPLOID-30x-HiFi-2

4.2M

21K

49K

Cheng et al. "Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm" (2021)

Nurk et al. "The complete sequence of a human genome" (2022)

A double repeat of
A doublelength `repeat
0 inof glength
is a substring` Indian
0 inInstitutexg 2is⌃asubstrnec
of Scienceappearing at distinct
appearing
positionsat idistinct
in g; i.e.,i1g[iand+ `g; i.e.,chag[
1 and i2positions
1 :ii21 in
1] = g[i1]1]==g[ix.2Similarly,: i2 +`
1]a triple= x. Similarly,
repeatof lengtha triple` repeatis ando
2 : i2 +`Coverage
gaps
observedbyremoving
contained
reads
asubstringxthatappearsa
substring
at
three
x
that
distinctappears
locations
at
three
in
distincts(possibly
locationsintec• Step 1: Identify contained reads by all-vs-all read alignments
overlapping); i.e.,overlapping);
g[i1 : i1 + ` i.e.,1] =g[ig[i+ ` 1] 1]==g[ig[i`
1 :2i1: i+
2`
2 :3i2: + sam• Step 2: map non-contained reads to genome
i3 + ` 1] = x fori3 +distinct` 1]i=
anddistinct
i3 (modulo
i1 , i2G,
andgiven
i3 (modulo
the andG
1 , ix2 for
Coverage-gaps
circularData
assumption
circularong).
assumption
If
x
is
a
doubleong).
repeat
If
x
is
but
a
doublenotatriplerepeatbuCount of contained reads
rec
Count
Maximum length
repeat, we say that
repeat,it is weprecisely
say thata double
it is precisely
repeat. aA double
double repeawh
HAPLOID-20x-ONT-1
3.2M
0
repeatxismaximalrepeat
if
it
x
is
is
not
maximalasubstring
if
it
is
of
not
any
a
strictly
substring
longerofanyst
HAPLOID-20x-ONT-2
3.2M
0
inf
double repeat. Finally,
doubleif 1.9M
repeat.
x = g [iFinally,[i[il]l]and= g [igiv
HAPLOID-20x-HiFi-1
1 : i1 +0ifl] x==g g2 1: i:-2i1++
2
0
00
0
0
HAPLOID-20x-HiFi-2
1.9M
0
y = g [j1 : j1 + ly] == gg[j[j12 :: jj12 ++ll ]] =forg some[j2 : j2i1+, i2l ,]l,for
j1 j2some,libe1,
where x, y are maximal
where x,andy are
i1 <maximal
j1 < i2and< ji21, <then
j1 x< andi2 <y j2 , thcontained
Reads
2
7
form an interleaved
form
repeat.
aninterleaved
Examplesrepeat.
areshown
Examples
in
Fig.
are
1.
shown
in
3
5
8
1

4

6

Haploid genome
25

Identical
maximal double repeat repeats maximal double
repeatrepeat
doubleLappdouble repeatA double repeat of
A doublelength `repeat
0 inof glength
is a substring` Indian
0 inInstitutexg 2is⌃asubstrnec
of Scienceappearing at distinct
appearing
positionsat idistinct
in g; i.e.,i1g[iand+ `g; i.e.,chag[
1 and i2positions
1 :ii21 in
1] = g[i1]1]==g[ix.2Similarly,: i2 +`
1]a triple= x. Similarly,
repeatof lengtha triple` repeatis ando
2 : i2 +`Coverage
gaps
observedbyremoving
contained
reads
asubstringxthatappearsa
substring
at
three
x
that
distinctappears
locations
at
three
in
distincts(possibly
locationsintec• Step 1: Identify contained reads by all-vs-all read alignments
overlapping); i.e.,overlapping);
g[i1 : i1 + ` i.e.,1] =g[ig[i+ ` 1] 1]==g[ig[i`
1 :2i1: i+
2`
2 :3i2: + sam• Step 2: map non-contained reads to genome
i3 + ` 1] = x fori3 +distinct` 1]i=
anddistinct
i3 (modulo
i1 , i2G,
andgiven
i3 (modulo
the andG
1 , ix2 for
Coverage-gaps
circularData
assumption
circularong).
assumption
If
x
is
a
doubleong).
repeat
If
x
is
but
a
doublenotatriplerepeatbuCount of contained reads
rec
Count
Maximum length
repeat, we say that
repeat,it is weprecisely
say thata double
it is precisely
repeat. aA double
double repeawh
HAPLOID-20x-ONT-1
3.2M
0
repeatxismaximalrepeat
if
it
x
is
is
not
maximalasubstring
if
it
is
of
not
any
a
strictly
substring
longerofanyst
HAPLOID-20x-ONT-2
3.2M
0
inf
double repeat. Finally,
doubleif 1.9M
repeat.
x = g [iFinally,[i[il]l]and= g [igiv
HAPLOID-20x-HiFi-1
1 : i1 +0ifl] x==g g2 1: i:-2i1++
2
0
00
0
0
HAPLOID-20x-HiFi-2
1.9M
0
y = g [j1 : j1 + ly] == gg[j[j12 :: jj12 ++ll ]] =forg some[j2 : j2i1+, i2l ,]l,for
j1 j2some,libe1,
where x, y are maximal
where x,andy are
i1 <maximal
j1 < i2and< ji21, <then
j1 x< andi2 <y j2 , th
Reads
2
7
form aninterleaved
form
repeat.
an interleaved
Examplesrepeat.
areshown
Examples
in
Fig.
are
1.
shown
in
5
8
1

gap

4

6

Haploid genome
26

Identical
maximal double repeat repeats maximal double
repeatrepeat
doubleLappdouble repeat

Indian Institute of Science

Coverage gaps observed by removing contained reads• Step 1: Identify contained reads by all-vs-all read alignments• Step 2: map non-contained reads to genome
Data

Count of contained reads

Count

Maximum length

DIPLOID-30x-ONT-1

4.6M

46

53K

DIPLOID-30x-ONT-2

4.6M

54

101K

DIPLOID-30x-HiFi-1

2.5M

1

2K

DIPLOID-30x-HiFi-2

2.5M

1

0.2K

Hap1 reads

Diploid genome

1

2

C

Hap2

A
6

7

4

3

Hap1

Hap2 reads
27

Coverage-gaps

8

5

10
9

containedIndian Institute of Science

Coverage gaps observed by removing contained reads• Step 1: Identify contained reads by all-vs-all read alignments• Step 2: map non-contained reads to genome
Data

Count

Maximum length

DIPLOID-30x-ONT-1

4.6M

46

53K

DIPLOID-30x-ONT-2

4.6M

54

101K

DIPLOID-30x-HiFi-1

2.5M

1

2K

DIPLOID-30x-HiFi-2

2.5M

1

0.2K

Hap1 reads

Diploid genome

1

2

C

Hap2

A
6

7

4

3

Hap1

Hap2 reads
28

Count of contained reads

Coverage-gaps

8

gap

5

10

Indian Institute of Science

Read length distributions
Count ofreads

Count ofreads
Read length

6e-05

Read length

0.4M

Read length

Read length

1.0M

0.2M
Frequency

Frequency

0.00012

3e-05

0.5M

6e-05

0

40 kbp40000

29

160 kbp
160000Count: 406014 Mean: 20518.7 Stdev: 24250.2 Min: 1000.0 Max: 320722.0 Precision: 1 Bins: 70

280 kbp280000

0

5 kbp
5000

25 kbp
25000

Count: 499734 Mean: 21316.3 Stdev: 3510.4 Min: 1000.0 Max: 47444.0 Precision: 1 Bins: 70

Read length

Read length

Oxford Nanopore (ONT)PacBio HiFi

45 kbp45000

Indian Institute of Science

Evaluation of proposed heuristics

Data

DIPLOID-30x-ONT-1

DIPLOID-30x-HiFi-1

MethodCount of contained Count of junction Gaps introduced inreads retained
verticesthe genome

ContainX

ContainX

LOWER ISBETTER
30

LOWER ISBETTER

Indian Institute of Science

Existing solutions besides ContainX• Other solutions to identify “useful” contained reads• [Hui et al. ISIT 2016]• Contained read is removed if it has an inconsistent pair of parent reads
CA
contained readinconsistent• Hi asm [Cheng et al. 2021]• Recovers contained reads which join a broken haplotype walkfi

31

Indian Institute of Science

Existing solutions besides ContainX• Other solutions to identify “useful” contained reads• [Hui et al. ISIT 2016]• Contained read is removed if it has an inconsistent pair of parent reads
CA
contained read

• Hi asm Hybrid: PacBio HiFi + ultra-long ONT

inconsistent[Cheng et al. 2023]• Identi es useful contained reads by aligning ultra-long nanopore reads to graphvertices from contained reads
fi

fi

32

Indian Institute of Science

Existing solutions besides ContainX• Other solutions to identify “useful” contained reads• [Hui et al. ISIT 2016]• Contained read is removed if it has an inconsistent pair of parent reads
CA
contained read

• Hi asm Hybrid: PacBio HiFi + ultra-long ONT

inconsistent[Cheng et al. 2023]• Identi es useful contained reads by aligning ultra-long nanopore reads to graphvertices from contained reads
fi

fi

33

Indian Institute of Science

Existing solutions besides ContainX• Other solutions to identify “useful” contained reads• [Hui et al. ISIT 2016]• Contained read is removed if it has an inconsistent pair of parent reads
CA
contained read

• Hi asm Hybrid: PacBio HiFi + ultra-long ONT

inconsistent[Cheng et al. 2023]• Identi es useful contained reads by aligning ultra-long nanopore reads to graphvertices from contained reads
fi

fi

34

Indian Institute of Science

Evaluation of proposed heuristics

Data

MethodCount of contained Count of junction Gaps introduced inreads retained
verticesthe genome

Retain all
DIPLOID-30x-ONT-1

Hui-2016
ContainX
Remove all
Retain all

DIPLOID-30x-HiFi-1

Hui-2016
ContainX
Remove all
LOWER ISBETTER

35

LOWER ISBETTER

Indian Institute of Science

Evaluation of proposed heuristics

Data

MethodRetain all

DIPLOID-30x-ONT-1

DIPLOID-30x-HiFi-1

Count of contained Count of junction Gaps introduced inreads retained
verticesthe genome
2.8M

2.5M

0

Remove all

0

38.9K

46

Retain all

2.5M

3.4M

0

0

158.4K1

Hui-2016
ContainX

Hui-2016
ContainX
Remove all

LOWER ISBETTER
36

LOWER ISBETTER

Indian Institute of Science

Evaluation of proposed heuristics

Data

DIPLOID-30x-ONT-1

DIPLOID-30x-HiFi-1

MethodCount of contained Count of junction Gaps introduced inreads retained
verticesthe genome

Retain all

2.8M

2.5M

0

Hui-2016

2.5M

2.3M

0

ContainX

28.5K

53.9K

2

Remove all

0

38.9K

46

Retain all

2.5M

3.4M

0

Hui-2016

2.5M

3.3M

0

ContainX

39.8K

184.1K

0

Remove all0

158.4K

1

LOWER ISBETTER
37

LOWER ISBETTER

Indian Institute of Science

Evaluation of proposed heuristics

Data

DIPLOID-30x-ONT-1

DIPLOID-30x-HiFi-1

Retain all

2.8M

2.5M

0

Hui-2016

2.5M

2.3M

0

ContainX

28.5K

53.9K

2

Hi asm

4.0K

1.7K

33

Remove all

0

38.9K

46

Retain all

2.5M

3.4M

0

Hui-2016

2.5M

3.3M

0

ContainX

39.8K

184.1K0

Hi asm

164

36.9K

0

Remove all0

158.4K1* Hi asm is an end-to-end genome assembler,uses multiple graph sparsi cation heuristics
fi

fi

fi

fi

38

MethodCount of contained Count of junction Gaps introduced inreads retained
verticesthe genome

LOWER ISBETTER

LOWER ISBETTER

Indian Institute of Science

Conclusions• Provably-good graph models will be useful for reliable and accurate
human genome reconstruction• String graph model is used commonly, but
it violates the ‘safety’ guarantee, both in theory and practice.

• Optimal sparsi cation of overlap graphs remains unsolved. We
proposed safe rules and promising heuristics.
chirag@iisc.ac.in
github.com/at-cg/ContainX

fi

39

Link to publication:
https://doi.org/10.1093/bioinformatics/btad124

Indian Institute of Science

Acknowledgement
Mehak Bindra (Project sta ), Sudhanva Shyam Kamath (Ph.D. student)
Valuable feedback from:• Prof. Sunil Chandran, Prof. Debnath Pal [IISc]• Haowen Zhang [Georgia Tech]• Mile Sikic, Robert Vaser [Genome Institute of Singapore]• Brian Walenz, Sergey Nurk [NHGRI]• Haoyu Cheng [Dana Farber Cancer Institute]ff

40