ANN - Towards a (Bio)Julia powered Genome Graphs framework

Was that in a fresh new project? Can I see what happens upon Pkg add? GenomeGraphs currently doesen’t depend on BioAlignments so I’m a bit confused without seeing the console msges.

It works in a fresh environment yes, and it’s indeed BioAlignments that is causing an incompatibility (I wanted to get reads from a bam file).

Hello lovely fellow Julians!

I haven’t left an update since January so I thought I had better fill you in on progress on the genome graph framework.

ReadDatastores has gotten an update and a new release pending which allows you to choose to create stores for LongSequence{DNAAlphabet{2}}s or LongSequence{DNAAlphabet{4}}s now, rather than forcing you to use only DNAAlphabet{4} - so if you know your reads dont have IUPAC ambiguities, this will likely save you a lot more disk-space and may be faster (haven’t benchmarked yet).

A bug fix has also made it’s way into ReadDatastores which should reduce the number of buffer refreshes that occur, when streaming sequences from a buffered datastore.

Now ReadDatastores is in a state I like a bit more, I’m going to first update GenomeGraphs to use the new ReadDatastores and begin breaking out some of the Kmer counter data-structures, as the ability to project kmer counts and coverage over sequences in general - not just graphs - is - I’m finding - an insanely useful thing to do!

4 Likes

Hey guys, some progress and updates on the genome assembly and genome graph analysis framework:

I took the time to start breaking out the kmer counting functionality - barebones as it was - into its own package MerCounting.jl.

Currently the serial and all in RAM method of counting kmers that is used for GenomeGraphs.jl’s deBruijn graph builder made it in, but I wanted a more serious kmer counter algorithm, or rather a more serious choice.

So I’m thinking of 3 counter methods to begin with: 1). The serial-in-RAM one (called serial_mem) that already exists, for when the problem is trivial or memory isn’t and issue - I have access to a cluster configured for large memory requirements rather than speed, and sometimes its quicker and easier to just throw a big vector at the thing rather than spend time on something more clever. 2). A distributed-in-RAM one that uses “Distributed” (called dist_mem) - again not so concerned with memory, but speeds up serial_mem by cutting the ReadDatastores into equal-ish chunks and processes them according to serial_mem and then combines the result. I guess if your julia processes were spread across multiple machines in a cluster or such, this could help with memory. 3). A distributed counter that uses disk batches. This one does a pre-analysis of a ReadDatastore to look at the distribution of minimizers, and will allocate N “bins” that the minimizers are distributed across as evenly as possible. Then multiple processes will divide up the reads, kmerize them, and throw the kmers into RemoteChannels (bins) based on their minimizer. These channels are read by a series of writer processes that send the kmer to a disk file.
These bin files can then be processes individually and in parallel, and unlike the other two counting methods, which return a vector of sroted counts, this will probably return a KmerCountDatastore or some other disk backed store of the kmer counts (similar to the ReadDatastores I guess).

Anyway here’s a preview of testing that dist_mem works:

julia> using Distributed, BioSequences, ReadDatastores

julia> addprocs(4, exeflags="--project=.")
4-element Array{Int64,1}:
 2
 3
 4
 5

julia> @everywhere using MerCounting

julia> ds = buffer(@openreads "ecoli-pe.prseq")
Buffered Paired Read Datastore 'ecoli-pe': 20 reads (10 pairs)

julia> o = MerCounting.Counters.dist_mem(DNAMer{31}, ds, MerCounting.CANONICAL)
[ Info: Splitting all reads in ecoli-pe, across 5 processes and counting kmers
[ Info: [Process: 1] Collecting kmers from 4 reads in read datastore ecoli-pe
[ Info: [Process: 1] Collapsing collected kmers into counts
[ Info: [Process: 2] Collecting kmers from 4 reads in read datastore ecoli-pe
[ Info: [Process: 2] Collapsing collected kmers into counts
[ Info: [Process: 3] Collecting kmers from 4 reads in read datastore ecoli-pe
[ Info: [Process: 3] Collapsing collected kmers into counts
[ Info: [Process: 4] Collecting kmers from 4 reads in read datastore ecoli-pe
[ Info: [Process: 4] Collapsing collected kmers into counts
[ Info: [Process: 5] Collecting kmers from 4 reads in read datastore ecoli-pe
[ Info: [Process: 5] Collapsing collected kmers into counts
[ Info: Merging all counts from 5 processes
4916-element Array{MerCount{Mer{DNAAlphabet{2},31}},1}:
 AAAAAAAAGTAATCGTCGGCATGTCCGGCGG occurs 1 times
 AAAAAAAGTAATCGTCGGCATGTCCGGCGGT occurs 1 times
 AAAAAAATAAAATTTAGTGCTGTACAGAGCG occurs 1 times
 AAAAAAGTAATCGTCGGCATGTCCGGCGGTG occurs 1 times
 AAAAAATAAAATTTAGTGCTGTACAGAGCGC occurs 1 times
 AAAAAATCCACCAGCAGCAGGTTGATGTCAG occurs 1 times
 AAAAAATCGACCAGGAGCGGTTGGATGACAG occurs 1 times
 AAAAAATCGTTGCGCTTTTTCTCGCCGCCGG occurs 1 times
 AAAAACATCGCGGGTGTTGAAGAAAAATTAA occurs 1 times
 AAAAACCAACGTTTACCGAACGGGTTAATAA occurs 2 times
 AAAAACTGGTTACTGACACCTGGAATCTCCA occurs 1 times
 AAAAACTTCCTCAAACCAATGCTGATTGTGT occurs 1 times
 AAAAAGCACTATTTCATCCTTGTTGCTGAAG occurs 1 times
 AAAAAGCGCAACGATATTTTGCAAATGGCGG occurs 1 times
 AAAAAGTAATCGTCGGCATGTCCGGCGGTGT occurs 1 times
 AAAAATAAAATTTAGTGCTGTACAGAGCGCG occurs 1 times
 AAAAATCCACCAGCAGCAGGTTGATGTCAGC occurs 1 times
 AAAAATCGACCAGGAGCGGTTGGATGACAGC occurs 1 times
 AAAAATCGTTGCGCTTTTTCTCGCCGCCGGA occurs 1 times
 AAAAATTAACCTATACCGATACCTACGCGCA occurs 1 times
 AAAACAATCCCTTCTATCAATTATATCGGCT occurs 1 times
 AAAACACCGTTATTCATTGCTTCTACCCGTG occurs 1 times
 AAAACAGCGTATTGAAAGCCTTCATGTAAAA occurs 1 times
 AAAACATCGCGGGTGTTGAAGAAAAATTAAC occurs 1 times
 AAAACCAACGTTTACCGAACGGGTTAATAAA occurs 2 times
 ⋮                                             
 TGTTAAGGACTATCTTGTTAAATGCTCGGAA occurs 1 times
 TGTTTATTCGCACCGGCGACGTAGACGGCAA occurs 2 times
 TTAACCTATACCGATACCTACGCGCAGGAAA occurs 1 times
 TTAAGAAGTTGTTAAGGACTATCTTGTTAAA occurs 1 times
 TTACACCGGCGGTGACCAAAAACTTCCTCAA occurs 1 times
 TTACGACCTTGTGTAAACCAGTGGACGCTAA occurs 1 times
 TTACGTGGTGATCACCATTACCCCGAGCGAA occurs 1 times
 TTACTAAGGCCCGTTTTGATGCCGTCGCCAA occurs 1 times
 TTACTGACGTATCATTATCTGGTAGAAAAAA occurs 1 times
 TTAGGGGCCAACGTCGATTGTGACAGCACAA occurs 1 times
 TTATCTGTGTAACCCTTTTGCGTTGATAGAA occurs 1 times
 TTATTCCCCATGCTTCAGCAACAAGGATGAA occurs 1 times
 TTCAATGTCGTGAGTGATCCAATGTCTGAAA occurs 1 times
 TTCAATTTGCCTCCATTGGTGCAACCACCAA occurs 1 times
 TTCCTCCGTTTCTTCCTGGCTGTAGCAACAA occurs 1 times
 TTCGATTCTTCTTTGTCACCGCAGCCAGCAA occurs 1 times
 TTCGGCGTGCGACCGGCTTTATATTCGGCAA occurs 1 times
 TTCTGGCTGAAGAGGTGGTGGAAATTCCCAA occurs 1 times
 TTCTTGCTTAAGAGGTGGTGGAAATTCCCAA occurs 1 times
 TTGCCATTACCTATGTCTACAAAGGTGACAA occurs 1 times
 TTGCGAACAATTATTCCCCATGCTTCAGCAA occurs 1 times
 TTGCGCTTTTTCTCGCCGCCGGAAAAACCAA occurs 2 times
 TTGGCAGCATCTTCTTTGGTGGTTGCACCAA occurs 1 times
 TTGTCCCAACACATGTGCCAGCCGATAAAAA occurs 1 times
 TTTGATTTTCAGGATTTGATGGAAGAGAAAA occurs 1 times

The third more complex counter is on its way - I have code already for the minimizer analysis for deciding on the disk bins.

3 Likes

Hi All,

I’ve been fairly quite for a while so I thought an update might be appropriate.

One of the components that is super important for an assembler / contigger that is truly generally haplotype specific is the ability to project kmer coverage and counts into the graph. Doing this allows the boundaries of assembly sub-problems to be defined and solved iteratively to resolve the genome and simplify the graph.

But projecting kmer coverage is a generally useful tool as well, as are other analyses for example computing kmer spectra.

So KmerCounting.jl is now KmerAnalysis.jl, and will have the following to mark an initial release:

Kmer counter types, that count kmers in an input, and behave somewhat like functors. Here’s an example of the serial and in memory counter now:

struct SerialMemCounter{M<:AbstractMer,C<:CountMode} <: Function
    mer_buffer::Vector{M}
    count_buffer::Vector{MerCount{M}}
end

serial_mem(::Type{M}, count_mode::C) where {M<:AbstractMer,C<:CountMode} = SerialMemCounter{M,C}(Vector{M}(), Vector{MerCount{M}}())

function count!(smc::SerialMemCounter{M,C}, input, args...) where {M<:AbstractMer,C<:CountMode}
    _serial_mem!(smc.count_buffer, smc.mer_buffer, C(), input, args...)
end

function (smc::SerialMemCounter{M,C})(input, args...) where {M<:AbstractMer,C}
    count!(smc, input, args...)
    return copy(smc.count_buffer)
end

function _serial_mem!(count_buffer::Vector{MerCount{M}}, mer_buffer::Vector{M}, count_mode::CountMode, input, args...) where {M<:AbstractMer}
    collect_mers!(mer_buffer, count_mode, input, args...)
    collapse_into_counts!(count_buffer, mer_buffer)
end

I think defining kmer counters this way has a few advantages. First they can be passed around and treated like functions that customize the behaviour of an analysis. For example to customize how the counting is done during a kmer spectra analysis:

c = serial_mem(DNAMer{31}, CANONICAL)
speccy = spectra(c, my_reads)

Where the spectra method is something like:

function spectra(counter::C, input, min_count::Integer = 0) where {C<:AbstractKmerCounter}
    return spectra(counter(input), min_count)
end

These counter types could also edit the behaviour of code they are passed through simply by altering dispatch too, not just by being called as a function.

I think this is quite a nice API for this kind of thing, but I know it is possible to produce something the compiler finds hard to optimize when passing functions about so any reviews of the package when it’s released are most welcome.

As I hinted the package will also provide functionality to compute kmer spectra. These are useful to check the distribution of kmer frequencies in the reads before attempting a genome assembly. For example:

A plot of a kmer spectra from reads like this shows me the main sampling distribution and the errors in the 1 column are very distinct cutting off the kmer count at a threshold of 2 or three for DBG construction should retain most of the information, and discard most error.

You can also have 2D kmer spectra. They are useful for producing spectra copy number plots
which should be your first check that an assembly, at any stage is any good. Indeed we tell people at our institute “Spectra CN or your assembly didn’t happen”

spec = spectra(kc_reads, kc_genom)

# Basic recipe for cn plot
heights = cumsum(KmerAnalysis._find_plottable_subset(spec); dims = 2)
s = barplot([Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,7])], color = :orange)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,6])], color = :yellow)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,5])], color = :blue)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,4])], color = :green)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,3])], color = :purple)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,2])], color = :red)
barplot!(s, [Point2f0(i - 1, j) for (i, j) in enumerate(heights[:,1])], color = :black)
xlabel!(s, "K-mer multiplicity")
ylabel!(s, "Number of distinct k-mers")

It needs a legend but the colours basically correspond to how many times a kmer appears in the genome, (0 = black, red = 1, orange = 2 etc…) but this CN plot tells me that since the first (and only) peak (which corresponds to the first (and only) haplotype) is coloured red, the kmers from that peak appear once in the genome. Then the part of the spectra corresponding to sequencing errors is black - they appear 0 times in the genome, which is great because before you consider contiguity or N50 or whatever, you should first make sure you assembly has the correct motifs, the right number of times, and in the right order. The cn plot validates the first two criteria.

So the spectra and cn plots will be important for people working with the GenomeGraphs API to check as they go that they haven’t lost any genome content or even made some up which does happen!

Genomes with more than one haplotype, diploids, triploids and so on, get a bit more complex but the cn plot is simple enough to interpret - every peak corresponds to a level of coverage 1x, 2x and so on, and so a plot from a good haplotype-specific assembly will colour each of the peaks the correct colour.

Here’s an example from a decent diploid assembly where the assembler collapses bubbles (heterozyogous content) and basically chooses one of the two variants. Note the black in the first peak is the bubble halfs the assembler is discarding.

Whereas a correct haplotype specific assembler would look something like this:

So the ability to do these checks is provided in KmerAnalysis.jl. This will be the first package I’ve done where I’ve had to really try to do some nice visualisation, and I’m liking Makie so for, there are a few things regarding recipe’s I don’t get or can’t do - I’ve found I can edit x and y labels in plots, I don’t know if I can add colour legends yet, because recipes work on plots, not scenes.

The other datatype I’ve been working on is IndexedCounts.
So the idea with indexed counts is you first construct one from some sequences you want to use as an index (normally a genome assembly or other reference genome). Then you add many counts from different data sources - short reads, long reads, reads from different individuals and so on.

julia> KmerAnalysis.index!(ikc, genome) # TODO make this happen during construction.
[ Info: Indexing kmer counts of reference sequences
[ Info: Populating index with reference sequence mers
[ Info: Sorting the index
[ Info: Counting and collapsing mers in index
31-mer Indexed Counts Datastore 'ecoli': 1 stored 31-mer counts

julia> # Add a count called :pe_reads to the indexed counts collection.
       add_count!(ikc, :pe_reads, reads)

You can then project the counts from an IndexedCount onto some new query sequence.
Here I project the counts from the reads, back onto the first sequence of the genome:

lines(project_count(ikc, :pe_reads, genome[1]))

Most of the sequence is covered by the read kmers between 10 and 20 times, and there are some very high portions, overrepresented or very common motifs perhaps, and some regions clearly more covered than others. Projecting kmer count data over genome graph sequences in this way, and also - more importantly - it’s transformed KCI or Kmer Count Index equivalent (more on that later) is a central and critical tool to both do and describe and assess a haplotype-specific assembly graph.

Anyway KmerAnalysis.jl is comming given its utility in GenomeGraphs.jl and elsewhere I hope people end up taking a look and don’t hesitate to improve it or point out things to be improved.

6 Likes