short read genome assembly Sorin Istrail CSCI1820 Short-read genome - - PowerPoint PPT Presentation

short read genome assembly

Sorin Istrail CSCI1820 Short-read genome assembly algorithms 3/6/2014

1

Genomathica Assembler

• Mathematica notebook for genome assembly simulation
• Assembler can be found at:

http://cs.brown.edu/courses/csci1820/software/minimal_asse mbler.nb

• Sample FASTA genome phix174.fasta can be found in HW5

Biology: http://cs.brown.edu/courses/csci1820/software/phix174.fasta

• Remember to

– Change the input genome to your FASTA file’s location – Evaluate each cell initially, then you only need to evaluate the last two cells to re-run the assembly, and display the results respectively – Mathematica can be downloaded here: http://www.brown.edu/information-technology/software/

2

coverage = 1

• Sequence reads

are in black

• Contiguous

strings of assembled DNA (contigs) are in red

coverage = 2

• Sequence reads

are in black

• Contiguous

strings of assembled DNA (contigs) are in red

coverage = 3

• Sequence reads

are in black

• Contiguous

strings of assembled DNA (contigs) are in red

coverage = 4

• Sequence reads

are in black

• Contiguous

strings of assembled DNA (contigs) are in red

coverage = 5

• Sequence reads

are in black

• Contiguous

strings of assembled DNA (contigs) are in red

coverage = 2, paired ends

Raw Sequence Reads Sample prep

Sequence data

• wet-lab experimental methods to isolate, prepare, and

sequence the DNA

• results in a number of large FASTQ files
• FASTQC can be used to check basic statistics of the files

–http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

• many tools available for QC

–e.g. http://hannonlab.cshl.edu/fastx_toolkit/

Wetterstrand KA. DNA Sequencing Costs: Data from the NHGRI Genome Sequencing Program (GSP) Available at: www.genome.gov/sequencingcosts. Accessed April 2013.

http://www.ncbi.nlm.nih.gov/Traces/sra/

Genome Assembly Software

• Overlap-layout-consensus
• Celera: http://wgs-assembler.sourceforge.net/
• K-mer based
• Velvet: http://www.ebi.ac.uk/~zerbino/velvet/
• SOAP-denovo: http://soap.genomics.org.cn/soapdenovo.html
• ALLPATHS-LG:

http://www.broadinstitute.org/software/allpaths-lg/blog/

• IDBA-UD: http://i.cs.hku.hk/~alse/hkubrg/projects/idba_ud/

Two graph models

• A first graph model

– Nodes (vertices) are contiguous sequences of k characters (k-mer) – Directed edge from vi to vj if vi[2..k]=vj[1..k-1]

A C G T T C ACG CGT GTT TTC

Two graph models

• De-bruijn graph

– Nodes (vertices) are contiguous sequences of k-1 characters (k-1-mer) – Directed edge from vi to vj if vi[1..k-1]+vj[k-1] are a valid k-mer

A C G T T C AC CG GT TT TC ACG CGT GTT TTC

Compeau et al. (2011) How to apply de Bruijn graphs to genome assembly

Note edges that are not reflected in the input!

Genome Assembly

• Building the k-mer graph

– nodes as k-mers, edges (k-1) overlap

17

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT

GACG ACGT

k=4 k=3

GAC ACG CGT CGTA GTA

1 1 1 1 1

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT

GACG ACGT

k=4 k=3

GAC ACG CGT CGTA GTA GTAC TACG TAC

1 1 1 1 2 1 1 1 1

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT

GACG ACGT

k=4 k=3

GAC ACG CGT CGTA GTA GTAC TACG TAC CGTT GTT

1 1 1 2 2 1 1 1 2 1 1

Genome Assembly

• Building the k-mer graph

– nodes as k-mers, edges (k-1) overlap – nodes as (k-1)-mers, edges form k-mers

21

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT k=4 k=3

GA AC CG GT TA GAC ACG CGT GTA

1 1 1 1 1 1 1

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT k=4 k=3

GA AC CG GT TA GAC ACG CGT GTA TAC

1 1 1 3 2 2 1 2 1 1

Genome assembly

Genome GACGTACGTT Reads GACGTA CGTACG TACGTT k=4 k=3

GA AC CG GT TA GT GAC ACG CGT GTA TAC GTT TT

1 2 1 4 2 2 1 2 2 1 1 2 1

Genome Assembly

• Building the k-mer graph

– G(k): nodes as k-mers, edges (k-1) overlap – H(k): nodes as (k-1)-mers, edges form k-mers

• H(k)=G(k-1)

– So it really does not matter which you choose to implement

• Where does the complexity come from?

– Sequencing errors, repeats, uneven coverage, contamination from other organisms, ploidy, unsequenced regions

25

Popping bubbles

Error occurs in the middle of a read and is propagated to many k-mers.

Trimming tips

Error creates an erroneous ending k-mer

Compeau et al. (2011) How to apply de Bruijn graphs to genome assembly

Chimeric extensions

Errors connect two nodes in the graph which do not correspond to a valid extension in the genome sequence

Repetitive regions

• Satellites, SINEs, LINEs
• Homologous Genes

– Ortholog: descended from the same ancestral sequence and separated by speciation – Paralog: genes created by a duplication event

29

30 Compeau et al. (2011) How to apply de Bruijn graphs to genome assembly

Compeau et al. (2011) How to apply de Bruijn graphs to genome assembly

Velvet assembler

• Four stages

– Hashing reads into k-mers – Constructing the de Bruijn graph (not all 4^k k- mers, only those that exist in input) – Correct errors – Resolve repeats

• But what after?

– Paper gives very little information on this...

32

The Chinese postman problem (CPP)

• Compute a closed tour of minimum length

that visits each edge at least once

– Similar to what we want except we may want to visit edges more than once due to repeats

• How do we deal with repeats?

– Also, the starting and ending vertices are distinct in genome assembly

• How can we convert the closed tour to an open one?

33

Your homework

• You are not required to implement section 4
• f

http://web.eecs.umich.edu/~pettie/matching /Edmonds-Johnson-chinese-postman.pdf

• You are not even required to model genome

assembly as CPP

• But you do have to build the k-mer graph,

correct errors, resolve repeats, and compute a CPP or Eulerian-like tour.

34

Evaluating assembly

• The Assemblathon2 study lists 102 measures for evaluating

assembly quality.

• Bradnam et al. (2013) Assemblathon 2: evaluating de novo

methods fo genome assembly in three vertebrate species

• 1. NG50 scaffold length: a length x where all scaffolds of length x or

longer consists of at least 50% of the genome size

• 2. NG50 contig length: a length x where all contigs of length x or

longer consists of at least 50% of the genome size

• 3. Amount of gene-sized scaffolds (>25 kbp). Useful for gene

finding.

• 4. CEGMA: Number of 458 core genes mapped

• 5. Fosmid coverage: How many validated fosmid regions were

captured in assembly

• 6. Fosmid validity: Percentage of assembly validated by validated

fosmid regions

• 7. Validated fosmid region tag scaffold summary score: number of

validated fosmid region tag pairs that match the same scaffold multiplied by the percentage of uniquely mapping tag pairs that map with correct distance. Rewards short-range accuracy.

• 8. and 9. Using local and global alignments of optimal map data,

how well the assembly is ordered.

• 10. REAPR summary score: a tool that evalutes accuracy of assembly

using paired reads

Evaluating assembly