参与者：佐藤健太 (@bicycle1885)
导师：丹尼尔·C·琼斯 (@dcjones)
感谢戈登和贝蒂·摩尔基金会的资助，我享受了由 NumFOCUS 管理的 2015 年 Julia 夏日编程，并前往波士顿参加了 2015 年 JuliaCon。在此期间，我创建了几个关于序列分析数据结构和算法的包，主要针对生物信息学。尽管 Julia 拥有大量用于浮点数数值计算的实用包，但它缺乏生物信息学中至关重要的有效且紧凑的数据结构。
高通量 DNA 测序仪的最新发展使得能够在一天内从生物样本中测序大量 DNA 片段（称为读段）。序列分析的第一步是确定这些片段在其他长参考序列中的位置，然后我们可以根据结果检测遗传变异或基因表达。此步骤称为序列比对或比对，由于参考序列最常见的是基因组规模（人类约为 32 亿长度），因此使用全文检索索引来加速此比对过程。这种全文检索索引在许多生物信息学工具中实现，最显著的是 bowtie2 和 BWA，它们的论文被引用了数千次。
我的项目的主要重点是创建一种在 Julia 中易于使用且在实际应用中高效的全文检索索引。为了实现这一目标，我创建了几个可用作其他数据结构构建块的包。我将在本文中向大家介绍这些包。
IntArrays.jl
IntArrays.jl 是一个用于无符号整数数组的包。那么，它有用吗？有用！这是因为在这个包中实现的 IntArray 类型可以尽可能小地存储整数。IntArray 类型有一个类型参数 w，它表示对数组中元素进行编码所需的位数。例如，如果每个元素都是 0 到 3 之间的整数，你只需要使用两位来编码它，并且 w 可以设置为 2 或更大。这些 2 位整数被打包到一个缓冲区中，因此数组消耗的空间仅为普通数组的四分之一。以下是 [0x01, 0x03, 0x02, 0x00] 字节序列的示例
index: 1 2 3 4
byte sequence (hex): 0x01 0x03 0x02 0x00
byte sequence (bin): 0b00000001 0b00000011 0b00000010 0b00000000
packed sequence (w=2): 01 11 10 00
in-memory layout: 00101101
完整的类型定义是 IntArray{w,T,n}，其中 w 如我所述是每个元素的位数，T 是元素的类型，n 是数组的维度。这种类型是 AbstractArray{T,n} 的子类型，并且将表现得像熟悉的数组一样；支持分配、随机访问和更新。IntVector 和 IntMatrix 也分别定义为类似于 Vector 和 Matrix 的类型别名。
以下是一个示例
julia> IntArray{2,UInt8}(2, 3)
2x3 IntArrays.IntArray{2,UInt8,2}:
0x00 0x00 0x01
0x00 0x00 0x03
julia> array = IntVector{2,UInt8}(6)
6-element IntArrays.IntArray{2,UInt8,1}:
0x00
0x00
0x03
0x03
0x02
0x00
julia> array[1] = 0x02
0x02
julia> array
6-element IntArrays.IntArray{2,UInt8,1}:
0x02
0x00
0x03
0x03
0x02
0x00
julia> sort!(array)
6-element IntArrays.IntArray{2,UInt8,1}:
0x00
0x00
0x02
0x02
0x03
0x03并且 IntArray 的内存占用要小得多
julia> sizeof(IntVector{2,UInt8}(1_000_000))
250000
julia> sizeof(Vector{UInt8}(1_000_000))
1000000由于在缓冲区中打包和解包整数需要额外的操作，因此操作中会有开销，而 IntArray 通常比 Array 慢。我已尽力将这种差异保持在尽可能小的范围内，但 IntArray 在排序时速度大约慢 4-5 倍
julia> array = rand(0x00:0x03, 2^24);
julia> sort(array); @time sort(array);
0.488779 seconds (8 allocations: 16.000 MB)
julia> iarray = IntVector{2}(array);
julia> sort(iarray); @time sort(iarray);
2.290878 seconds (18 allocations: 4.001 MB)如果你有改善性能的好主意，请告诉我！
IndexableBitVectors.jl
下一个包是 IndexableBitVectors.jl。你一定熟悉标准库中的 BitVector 类型；我的包中定义的类型是它的静态但可索引的版本。这里的“可索引”意味着可以**在常数时间内**回答查询以询问任意范围内的位数。如果你已经熟悉 简洁数据结构，你可能知道这是其他简洁数据结构（如小波树、LOUDS 等）的重要构建块。
该包导出了两种此类位向量的变体：SucVector 和 RRR。SucVector 比 RRR 更简单、更快，但 RRR 可压缩，如果 0/1 位在位向量中局部化，则 RRR 将更小。两种类型都将位向量分成块，并缓存到该位置的位数。在 SucVector 中，额外空间大约为每位 1/4 位，因此它将比原始位向量大约 25%。
对这些数据结构最重要的查询操作是 rank1(bv, i) 查询，它计算 bv[1:i] 中 1 位的个数。由于缓存了位计数，因此我们可以在常数时间内完成排名操作
julia> using IndexableBitVectors
julia> bv = bitrand(2^30);
julia> function myrank1(bv, i) # count ones by loop
r = 0
for j in 1:i
r += bv[j]
end
return r
end
myrank1 (generic function with 1 method)
julia> myrank1(bv, 2^29); @time myrank1(bv, 2^29);
0.714866 seconds (6 allocations: 192 bytes)
julia> sbv = SucVector(bv);
julia> rank1(sbv, 2^29); @time rank1(sbv, 2^29); # much faster!
0.000003 seconds (6 allocations: 192 bytes)
julia> rrr = RRR(bv);
julia> rank1(rrr, 2^29); @time rank1(rrr, 2^29); # much faster, too!
0.000004 seconds (6 allocations: 192 bytes)select1(bv, j) 查询在许多情况下也很有用，它在位向量 bv 中定位第 j 个 1 位。例如，如果一组正整数在该位向量中表示，你可以有效地查询集合中第 j 个最小的成员。
让我们看看 SucVector 的内部表示来理解其中的奥秘。一个位向量被分成较大的块
type SucVector <: AbstractIndexableBitVector
blocks::Vector{Block}
len::Int
end每个大块包含 256 位，由四个分别包含 64 位的小块组成，一个大块存储到其起始位置为止的全局 1 的计数，一个小块存储从其父大块的起始位置开始的局部 1 的计数。位本身存储在与小块相对应的四个位块中
immutable Block
# large block
large::UInt32
# small blocks
# the first small block is used for 8-bit extension of the large block
# hence, 40 (= 32 + 8) bits are available in total
smalls::NTuple{4,UInt8}
# bit chunks (64bits × 4 = 256bits)
chunks::NTuple{4,UInt64}
end
由于第一个小块的位计数始终为零，因此我们可以利用此空间来扩展大块的缓存（红色框）。在运行 rank1(bv, i) 查询时，它首先选择第 i 个位所属的大块和小块对，然后将它们的缓存位计数加起来，最后动态地计算剩余的 1 位。
正如我所说，这种数据结构可以用作各种数据结构的构建块。我将介绍的下一个包就是其中之一。
WaveletMatrices.jl
你可能已经了解 小波树，它支持类似于 SucVector 和 RRR 的排名和选择查询，但元素不限于 0/1 位。实际上，排名和选择查询在任意无符号整数上都可用。小波树可以被认为是可索引位向量的泛化。我实现的不是众所周知的小波树，而是一种称为“小波矩阵”的变体。你可以在 WaveletMatrices.jl 中找到实现和论文链接。根据论文作者的说法，小波矩阵“比逐层小波树更易于构建、更易于查询，并且在实践中更快”。
WaveletMatrix 类型需要三个类型参数：w、T 和 B。w 和 T 与 IntArray{w,T,n} 中的类似，而 B 是可索引位向量的类型。
julia> using WaveletMatrices
julia> wm = WaveletMatrix{2}([0x00, 0x01, 0x02, 0x03])
4-element WaveletMatrices.WaveletMatrix{2,UInt8,IndexableBitVectors.SucVector}:
0x00
0x01
0x02
0x03
julia> wm[3]
0x02
julia> rank(0x02, wm, 2)
0
julia> rank(0x02, wm, 3)
1
julia> xs = rand(0x00:0x03, 2^16);
julia> wm = WaveletMatrix{2}(xs); # 2-bit encoding
julia> sum(xs[1:2^15] .== 0x03)
8171
julia> rank(0x03, wm, 2^15)
8171数据结构和算法的细节相对简单，但超出了本文的范围。对于对这种数据结构感兴趣的人来说，我上面提到的论文和我的实现会有所帮助。小波矩阵还可以有效地运行更多操作，这些操作将在未来添加。
FMIndexes.jl
生物信息学中 80% 的序列分析都是关于序列搜索的，包括模式搜索、同源基因搜索、基因组比较、短读比对等等。 FM-Index 通常被认为是最有效的全文搜索索引之一，我已在 FMIndexes.jl 包中实现了它。由于我迄今为止介绍的包，它的代码看起来非常简单。例如，统计文本中给定模式出现的次数可以写成如下（为了说明目的，稍作简化）
function count(query, index::FMIndex)
sp, ep = 1, length(index)
# backward search
i = length(query)
while sp ≤ ep && i ≥ 1
char = convert(UInt8, query[i])
c = index.count[char+1]
sp = c + rank(char, index.bwt, sp - 1) + 1
ep = c + rank(char, index.bwt, ep)
i -= 1
end
return length(sp:ep)
endFM-Index 的一个独特特性是索引本身只是原始文本字符的排列和其中包含的字符的计数。这种排列被称为 Burrows-Wheeler 变换（也称为 BWT），排列后的文本存储在一个小波矩阵（或小波树）中，以便有效地统计特定区域内字符的个数。因此，对文本建立索引所需的空间通常比其他全文索引更小（实际上，在实践中，有效地找到查询的位置还需要辅助数据）。此外，这种变换是 双射 的，因此原始文本可以从索引中恢复。
为全文搜索建立索引非常简单：只需将序列传递给构造函数即可
julia> using FMIndexes
julia> fmindex = FMIndex("abracadabra");FMIndex 类型支持两个主要查询：count 和 locate。count(query, index) 查询实际上统计了 query 字符串出现的次数，而 locate(query, index) 则定位 query 的起始位置。为了恢复原始文本，可以使用 restore 函数。以下是一个简单的用法
julia> count("a", fmindex)
5
julia> count("abra", fmindex)
2
julia> locate("a", fmindex) |> collect
5-element Array{Any,1}:
11
8
1
4
6
julia> locate("abra", fmindex) |> collect
2-element Array{Any,1}:
8
1
julia> bytestring(restore(fmindex))
"abracadabra"例如，对于生物信息学家，让我们尝试对染色体进行一些查询。你还需要安装 Bio.jl 包才能有效地解析 FASTA 文件。下面的脚本从 FASTA 文件中读取染色体，构建一个 FM-Index，然后将其序列化到文件中以备后用（我喜欢 Julia 的序列化程序，它们是免费提供的！）
index.jl
using Bio.Seq
using IntArrays
using FMIndexes
# encode a DNA sequence with 3-bit unsigned integers;
# this is because a reference genome has five nucleotides: A/C/G/T/N.
function encode(seq)
encoded = IntVector{3,UInt8}(length(seq))
for i in 1:endof(seq)
encoded[i] = convert(UInt8, seq[i])
end
return encoded
end
# read a chromosome from a FASTA file
filepath = ARGS[1]
record = first(open(filepath, FASTA))
println(record.name, ": ", length(record.seq), "bp")
# build an FM-Index
fmindex = FMIndex(encode(record.seq))
# save it in a file
open(string(filepath, ".index"), "w+") do io
serialize(io, fmindex)
end好的，那么为人类的 22 号染色体创建一个索引（你可以从 这里 下载它）
$ julia4 index.jl chr22.fa
chr22: 50818468bp
$ ls -lh chr22.fa.index
-rw-r--r--+ 1 kenta staff 74M 9 26 06:30 chr22.fa.index构建完成后（这将需要几分钟），在 REPL 中读取索引
julia> using FMIndexes
julia> fmindex = open(deserialize, "chr22.fa.index");现在你可以执行查询来搜索 DNA 片段
julia> using Bio.Seq
julia> count(dna"GACTTTCAC", fmindex) # this DNA fragment hits at 111 locations
111
julia> count(dna"GACTTTCACTTT", fmindex) # this hits at 3 locations
3
julia> locate(dna"GACTTTCACTTT", fmindex) |> collect # the loci of these hits
3-element Array{Any,1}:
36253071
47308573
34159872
julia> count(dna"GACTTTCACTTTCCC", fmindex) # found a unique hit!
1
julia> locate(dna"GACTTTCACTTTCCC", fmindex) |> collect
1-element Array{Any,1}:
36253071
julia> @time locate(dna"GACTTTCACTTTCCC", fmindex); # this can be located in 32 μs!
0.000032 seconds (5 allocations: 192 bytes)这个基因座，chr22:36253071，是APOL1 基因的起始位置。
应用
我创建这些包的目的是证明在 Julia 中实现用于生物信息学的高性能数据结构是可行的。我确信这是真的，但其他人可能持怀疑态度。因此，我将通过使用这些包编写有用且高效的应用程序来证明这一点。现在我正在开发 FMM.jl，它使用 FM-Index 和其他算法将大量 DNA 片段比对到基因组序列。这仍然是一个正在进行的工作，可能会有许多错误和特殊情况需要我注意，但它的性能与其他实现相比并不差。
BioJulia 项目也正在积极开发中。我制作的包旨在与 Bio.jl 包一起使用。如果你对 BioJulia 项目感兴趣，我们热烈欢迎你的贡献！