By: JuliaLang - The Julia programming language
Re-posted from: https://julialang.org/blog/2016/09/biojulia2016-mid/
BioJulia 2016 – online sequence search, sequence demultiplexing, new readers and much more! | We are pleased to announce releasing…
Read more
By: Julia Developers
Re-posted from: http://feedproxy.google.com/~r/JuliaLang/~3/12oUw98qs-Q/biojulia2016-mid
We are pleased to announce releasing
Bio.jl 0.4, a minor release including
significant functionality improvements as I promised in the previous blog
post.
The following features are added since the post:
And many other miscellaneous performance and usability improvements! Tutorial
notebooks are available at https://github.com/BioJulia/BioTutorials. Here I
briefly introduce you to these new features one by one.
Sequence search is an indispensable tool in sequence analysis. Since the last
post, I have added exact, approximate and regex search algorithms. The search
interface of Bio.jl mimics that of Julia’s standard library.
julia> using Bio.Seq
julia> seq = dna"ACAGCGTAGCT"
11nt DNA Sequence:
ACAGCGTAGCT
# Exact search.
julia> search(seq, dna"AGCG")
3:6
# Approximate search with one error or less.
julia> approxsearch(seq, dna"AGGG", 1)
3:6
# Regular expression search.
julia> search(seq, biore"AGN*?G"d)
3:6
In Bio.jl DNA sequences are encoded using 4 bits per base by default in order to
store ambiguous nucleotides and this encoding does well in most cases. However,
some biological sequences such as chromosomal sequences are so long especially
for eukaryotic organisms and the default DNA sequences may result in a waste of
memory space. ReferenceSequence is a new type introduced in Bio.jl that
compresses positions of ambiguous nucleotides using a sparse bit vector. This
type can achieve almost 2-bit encoding in common reference sequences because
most of the ambiguous nucleotides are clustered in a sequence and the number of
them is small compared to other unambiguous nucleotides.
# Converting a DNASequence object to ReferenceSequence.
julia> ReferenceSequence(dna"ACGT"^10000)
40000nt Reference Sequence:
ACGTACGTACGTACGTACGTACGTACGTACGTACGTACG…CGTACGTACGTACGTACGTACGTACGTACGTACGTACGT
# Reading chromosome 1 of human from a FASTA file.
julia> open(first, FASTAReader{ReferenceSequence}, "hg38.fa")
Bio.Seq.SeqRecord{Bio.Seq.ReferenceSequence,Bio.Seq.FASTAMetadata}:
name: chr1
sequence: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN…NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
metadata: Bio.Seq.FASTAMetadata("")
julia> sequence(ans)
248956422nt Reference Sequence:
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN…NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
2bit is a binary file format to store reference sequences. This is a kind of
binary counterpart of FASTA but
specialized for DNA reference sequences to enable smaller file size and faster
loading. Reference sequences of various organisms are distributed from the
download page of UCSC in this
file format. An important advantage of 2bit is that sequences are indexed by its
name and can be accessed immediately.
# Opening a sequence file of yeast (S.cerevisiae).
julia> reader = open(TwoBitReader, "sacCer3.2bit");
# Loading a chromosome VI using random access index.
julia> reader["chrVI"]
Bio.Seq.SeqRecord{Bio.Seq.ReferenceSequence,Array{UnitRange{Int64},1}}:
name: chrVI
sequence: GATCTCGCAAGTGCATTCCTAGACTTAATTCATATCTGC…GTGTGGTGTGTGGGTGTGGTGTGTGGGTGTGGTGTGTGG
metadata: UnitRange{Int64}[]
The SAM and BAM file formats are designed for storing sequences aligned to
reference sequences. SAM is a line-oriented text file format and easy to handle
with UNIX command line tools. BAM is a compressed binary version of SAM and
suitable for storing data in disks and processing with purpose-built softwares
like samtools. The BAM data reader is carefully
tuned so that users can use it in real analysis with large files. It is also
feasible to read a CRAM file
combining the BAM reader and samtools view command.
An experimental feature is parallel processing using multiple threads.
Multi-threading support is introduced in Julia 0.5 and we use it to parallelize
decompression of BAM files. Here is a simple benchmark script to show how
much reading speed can be improved with multiple threads:
using Bio.Align
# Count the number of mapped records.
function countmapped(reader)
ret = 0
record = BAMRecord()
while !eof(reader)
# in-place reading
read!(reader, record)
if ismapped(record)
ret += 1
end
end
return ret
end
println(open(countmapped, BAMReader, ARGS[1]))
JULIA_NUM_THREADS environment variable controls the number of worker threads.
The result below shows that the elapsed time is almost halved using two threads:
~/.j/v/Bio $ time julia countmapped.jl SRR1238088.sort.bam
28040186
29.27 real 28.64 user 0.66 sys
~/.j/v/Bio $ env JULIA_NUM_THREADS=2 time julia countmapped.jl SRR1238088.sort.bam
28040186
17.40 real 32.31 user 0.63 sys
BGZF (Blocked GZip Format) is a gzip-compliant file format commonly used in
bioinformatics. BGZF can be read using standard gzip tools but files in the
format are compressed block by block and special metadata are added to index the
compressed files for random access. BAM files are compressed in this file format
and sequence alignments in a specific genomic region can be retrieved
efficiently. BGZFStreams.jl is a new package to handle BGZF files like usual
I/O streams and it is built on top of our
Libz.jl package. Parallel decompression
mentioned above is implemented in this package layer.
julia> using BGZFStreams
julia> stream = BGZFStream("/Users/kenta/.julia/v0.5/BGZFStreams/test/bar.bgz")
BGZFStreams.BGZFStream{IOStream}(<mode=read>)
julia> readstring(stream)
"bar"
https://github.com/BioJulia/BGZFStreams.jl
Sequence demultiplexing is a technique to distinguish the origin of a sequence
using its artificially-attached “barcode” sequence. This is often used at a
preprocessing phase after multiplexed
sequencing,
a common technique to sequence multiple samples simultaneously. A barcode
sequence, however, may be corrupted due to sequencing error, and we need to find
the best matching barcode from a barcode set. The demultiplexer algorithm
implemented in Bio.jl is based on a trie-like data structure, and efficiently
finds the optimal barcode from the prefix of a given DNA sequence.
# Set DNA barcode pool.
julia> barcodes = DNASequence["ATGG", "CAGA", "GGAA", "TACG"];
# Create a sequence demultiplexer that allows one errors at most.
julia> dplxr = Demultiplexer(barcodes, n_max_errors=1, distance=:hamming)
Bio.Seq.Demultiplexer{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}}}:
distance: hamming
number of barcodes: 4
number of correctable errors: 1
# Demultiplex a given sequence from its prefix.
julia> demultiplex(dplxr, dna"ATGGCGNT") # 1st barcode with no errors
(1,0)
julia> demultiplex(dplxr, dna"CAGGCGNT") # 2nd barcode with one error
(2,1)
This is still the first half of my project this year. The next term will come
with:
And, of course, improving existing features of Bio.jl and other packages. We
welcome any contributions and feature requests from you all. Gitter chat
channel is the best place to communicate
with developers and other users. If you love Julia and/or biology, any reason
not to join us?
I gratefully acknowledge the Moore Foundation and the Julia project for
supporting the BioJulia project. I also would like to thank Ben J.
Ward and Kevin
Murray for comments on my program code and other
contributions.
By: Julia Developers
Re-posted from: http://feedproxy.google.com/~r/JuliaLang/~3/-XJNyq8rvJc/GSoC2016-Graft
This blog post describes my work on Graft.jl, a general purpose graph analysis package for Julia. For those unfamiliar with graph algorithms, a quick introduction might help.
My proposal, titled ParallelGraphs, was to develop a parallelized/distributed graph algorithms
library. However, in the first month or so, we decided to work towards a more general framework that supports data analysis on
networks (graphs with attributes defined on vertices and edges).
Our change in direction was mainly motivated by:
The modified proposal could be summarized as the development of a package that supports:
LightGraphsParallelGraphs turned out to be a misnomer, since we were moving towards a more general purpose data analysis framework. So we chose the name Graft, a kind of abbreviation for Graph Toolkit. The following sections detail Graft's features:
Graphs are often representations of real world entities, and the relationships between them. Such entities (and their relationships), often have data attached to them.
While it is quite straightforward to store vertex data (a simple table will suffice), storing edges and their data is very tricky. The data should be structured on the
source and target vertices, should support random access and should be vectorized for queries.
At first we tried placing the edge data in a SparseMatrixCSC. This turned out to be a bad idea, because sparse matrices are designed for numeric storage.
A simpler solution is to store edge metadata in a DataFrame, and have a SparseMatrixCSC map edges onto indices for the DataFrame. This strategy needed a lot less
code, and the benchmarks were more promising. Mutations such as the addition or removal of vertices and edges become more complicated however.
Most graph libraries do not support vertex labelling. It can be very confusing to refer to a vertex by its (often long) integer identifier. It is also
computationally expensive to use non-integer labels in the implementation of the package (any such implementation would involve dictionaries). There is no reason, however,
for the user to have to use integer labels externally. Graft supports two modes of vertex labelling. By default, a vertex is identified by its internal identifier. A user
can assign labels of any arbitrary Julia type to identify vertices, overriding the internal identifiers. This strategy, we feel, makes a reasonable compromise between
user experience and performance.
If vertex labels were used in the internal implementation, the graph data structure would probably look like this:
Dict(
"Alice" => Dict(
"age" => 34,
"occupation" => "Doctor",
"adjacencies" => Dict("Bob" => Dict("relationship" => "follow")))
),
"Bob" => Dict(
"age" => 36,
"occupation" => "Software Engineer",
"adjacencies" => Dict("Charlie" => Dict("relationship" => "friend"))
),
"Charlie" => Dict(
"age" => 30,
"occupation" => "Lawyer",
"adjacencies" => Dict("David" => Dict("relationship" => "follow"))
),
"David" => Dict(
"age" => 29,
"occupation" => "Athlete",
"adjacencies" => Dict("Alice" => Dict("relationship" => "friend"))
)
)
Cleary, using labels internally is a very bad idea. Any sort of data access would set off multiple dictionary look-ups. Instead, if a bidirectional map
could be used to translate labels into vertex identifiers and back, the number of dictionary lookups could be reduced to one. The data would also be better structured for query processing.
# Label Map to resolve queries
LabelMap(
# Forward map : labels to vertex identifiers
Dict("Alice" => 1, "David" => 4", "Charlie" => 3, Bob" => 2),
# Reverse map : vertex identifiers to labels
String["Alice, "Bob", "Charlie", "David"]
)
# Vertex DataFrame
4×2 DataFrames.DataFrame
│ Row │ age │ occupation │
├─────┼─────┼─────────────────────┤
│ 1 │ 34 │ "Doctor" │
│ 2 │ 36 │ "Software Engineer" │
│ 3 │ 30 │ "Lawyer" │
│ 4 │ 29 │ "Athlete" │
# SparseMatrixCSC : maps edges onto indices into Edge DataFrame
4×4 sparse matrix with 4 Int64 nonzero entries:
[4, 1] = 1
[1, 2] = 2
[2, 3] = 3
[3, 4] = 4
# Edge DataFrame
4×1 DataFrames.DataFrame
│ Row │ relationship │
├─────┼──────────────┤
│ 1 │ "follow" │
│ 2 │ "friend" │
│ 3 │ "follow" │
│ 4 │ "friend" │
Graft’s query notation is borrowed from Jplyr. The @query macro is used to simplify the query syntax, and
accepts a pipeline of abstractions separated by the pipe operator |>. The stages are described through abstractions:
Accepts an expression, that is run over every vertex. Vertex properties can be expressed using the dot notation. Some reserved properties are v.id, v.label,
v.adj, v.indegree and v.outdegree.
Examples:
# Check if the user has overridden the default labels
julia> @query(g |> eachvertex(v.id == v.label)) |> all
# Kirchoff's law :P
julia> @query(g |> eachvertex(v.outdegree - v.indegree)) .== 0
Accepts an expression, that is run over every edge. The symbol s is used to denote
the source vertex, and t is used to denote the target vertex in the edge. The symbol e is used to denote
the edge itself. Edge properties can be expressed through the dot notation. Some reserved properties are e.source, e.target, e.mutualcount, and e.mutual.
Examples:
# Arithmetic expression on edge, source and target properties
julia> @query g |> eachedge(e.p1 - s.p1 - t.p1)
# Check if constituent vertices have the same outdegree
julia> @query g |> eachedge(s.outdegree == t.outdegree)
# Count the number of "mutual friends" between the source and target vertices in each edge
julia> @query g |> eachedge(e.mutualcount)
Accepts vertex or edge expressions and computes subgraphs with a subset of vertices, or a subset
of edges, or both.
Examples:
# Remove vertices where property p1 equals property p2
@query g |> filter(v.p1 != v.p2)
# Remove self loops from the graph
@query g |> filter(e.source != e.target)
Returns a subgraph with a subset of vertex properties, or a subset of edge properties or both.
Examples:
# Preserve vertex properties p1, p2 and nothing else
@query g |> select(v.p1, v.p2)
# Preserve vertex property p1 and edge property p2
@query g |> select(v.p1, e.p2)
The typical workflow we hope to support with Graft is:
– Load a graph from memory
– Use the query abstractions to construct new vertex/edge properties or obtain subgraphs.
– Run complex queries on the subgraphs, or export data to LightGraphs and run computationally expensive algorithms there.
– Bring the data back into Graft as a new property, or use it to modify the graphs structure.
The following examples should demonstrate this workflow:
More information can be found here
This work was carried out as part of the Google Summer of Code program, under the guidance of mentors: Viral B Shah and Shashi Gowda.