用julia语言计算测序数据的Insert Size?

Julia 读取 BAM 的库

想要计算Insert size，需要提供一个基因组比对后的文件，sam也好，bam也罢。那么，使用 julia 语言计算该值的第一步便是了解如何读取和解析BAM文件格式。

我们使用BioJulia提供的XAM包来读取 BAM 文件。所以我们需要首先安装该包。
打开 julia，输入]进入 Pkg 模式

add XAM

使用 XAM 读取和解析 BAM 文件的一般格式

reader = open(BAM.Reader, "data.bam")
record = BAM.Record()
while !eof(reader)
    empty!(record)
    read!(reader, record)
    # 做一些事情，例如解析和计算
end

上面的流程总体上就是

使用 BAM.Reader 读取文件
定义一个空的 Record 对象
从头至尾循环整个 BAMfile
将每行 bam 读入 Record

这样的操作方式可以节省内存，避免循环很大的 bam 文件爆内存。

什么是 Insert Size？

通俗的讲，Insert 长度就是指双端序列比对后，模版的长度。所以我们要计算需要保证如下条件

reads 是成对，最好是 Proper paired
只需要计算 read1 即可，不然就算重了
即使 Proper paired 也会有负的 Insert，需要移除

代码如下

using XAM
using Statistics

function insert_size_dist(Reader::XAM.BAM.Reader)
    insert_length = Int64[]
    record = BAM.Record()
    while !eof(Reader)
        empty!(record)
        read!(Reader, record)
        if BAM.flag(record) & 0x2 != 0 ## paired
            if BAM.flag(record) & 0x40 != 0  ## first in pair
                t_len = BAM.templength(record)
                if t_len > 0
                    push!(insert_length, t_len)
                end
            end
        end
    end
    (mean(insert_length), std(insert_length))
end