我们很高兴地宣布发布 Bio.jl 0.4,这是一个小版本,包含了重要的功能改进,正如我在 之前的博客文章 中承诺的那样。
自上次发布以来,添加了以下功能
在线序列搜索算法。
参考基因组的序列数据结构。
.2bit 文件格式的数据读取器和写入器。
SAM 和 BAM 文件格式的数据读取器和写入器。
序列拆分工具。
处理 BGZF 文件的包。
以及许多其他杂项性能和可用性改进!教程笔记本可在 https://github.com/BioJulia/BioTutorials 获取。在这里,我将简要地向您逐一介绍这些新功能。
序列搜索是序列分析中不可或缺的工具。自上次发布以来,我添加了精确、近似和正则表达式搜索算法。Bio.jl 的搜索接口模仿了 Julia 标准库的接口。
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
在 Bio.jl 中,DNA 序列默认使用每碱基 4 位进行编码,以便存储模糊核苷酸,这种编码在大多数情况下都能很好地工作。但是,一些生物序列(例如染色体序列)非常长,特别是对于真核生物,默认的 DNA 序列可能会导致内存空间浪费。ReferenceSequence
是 Bio.jl 中引入的一种新类型,它使用稀疏位向量压缩模糊核苷酸的位置。这种类型可以在常见的参考序列中实现几乎 2 位的编码,因为大多数模糊核苷酸都聚集在一个序列中,并且与其他明确的核苷酸相比,它们的数量很少。
# 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 是一种用于存储参考序列的二进制文件格式。这是一种 FASTA 的二进制对应格式,但专门用于 DNA 参考序列,以实现更小的文件大小和更快的加载速度。各种生物体的参考序列以这种文件格式从 UCSC 的下载页面 分发。2bit 的一个重要优势是序列按其名称索引,可以立即访问。
# 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}[]
SAM 和 BAM 文件格式旨在用于存储与参考序列比对的序列。SAM 是一种面向行的文本文件格式,易于使用 UNIX 命令行工具处理。BAM 是 SAM 的压缩二进制版本,适合用于将数据存储在磁盘中并使用专门构建的软件(如 samtools)进行处理。BAM 数据读取器经过精心调整,以便用户可以在大型文件的实际分析中使用它。还可以结合使用 BAM 读取器和 samtools view
命令来读取 CRAM 文件。
一个实验性功能是使用多个线程进行并行处理。Julia 0.5 中引入了多线程支持,我们使用它来并行化 BAM 文件的解压缩。这是一个简单的基准测试脚本,用于显示使用多个线程可以提高多少读取速度
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
环境变量控制工作线程的数量。下面的结果显示,使用两个线程后,经过的时间几乎减少了一半
~/.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(块式 Gzip 格式)是一种与 gzip 兼容的文件格式,通常用于生物信息学。BGZF 可以使用标准的 gzip 工具读取,但该格式的文件是逐块压缩的,并且添加了特殊的元数据来索引压缩文件以进行随机访问。BAM 文件以这种文件格式压缩,并且可以有效地检索特定基因组区域中的序列比对。BGZFStreams.jl 是一个用于处理 BGZF 文件(如普通 I/O 流)的新包,它建立在我们自己的 Libz.jl 包的基础之上。上面提到的并行解压缩是在此包层中实现的。
julia> using BGZFStreams
julia> stream = BGZFStream("/Users/kenta/.julia/v0.5/BGZFStreams/test/bar.bgz")
BGZFStreams.BGZFStream{IOStream}(<mode=read>)
julia> readstring(stream)
"bar"
序列拆分是一种使用人工附加的“条形码”序列来区分序列来源的技术。这通常在 多重测序(一种同时测序多个样本的常用技术)后的预处理阶段使用。然而,条形码序列可能会因测序错误而损坏,我们需要从条形码集中找到最佳匹配的条形码。Bio.jl 中实现的拆分器算法基于类似于 Trie 的数据结构,并有效地从给定 DNA 序列的前缀中找到最佳条形码。
# 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)
这仍然是我今年项目的上半部分。下一阶段将包含
支持更多文件格式,包括 GFF3、VCF 和 BCF。
与数据库集成。
与基因组浏览器集成。
当然,还要改进 Bio.jl 和其他包的现有功能。我们欢迎大家提出任何贡献和功能请求。 Gitter 聊天频道 是与开发人员和其他用户沟通的最佳场所。如果你喜欢 Julia 和/或生物学,还有什么理由不加入我们呢?
我衷心感谢摩尔基金会和 Julia 项目对 BioJulia 项目的支持。我还想感谢 Ben J. Ward 和 Kevin Murray 对我的程序代码的评论和其他贡献。