windows 也能快速处理 BAM/SAM 大型文件

测试开头

没关注？伸出手指点这里—

1引言

面对高通量测序数据,我们会经常碰到处理 BAM/SAM 这些大型文件,python 则有 pysam 库支持处理,但是只能在 linux 系统里安装,windows 里基本上装不了,感觉还是不太有好的。

如果是 sam 文件则可以使用 for 循环逐行处理应付, 但是 bam 就不行了。Rsamtools 包则可以帮助我们在 R 里面处理 bam/sam 文件,但是遇到 大文件,多文件 则效率明显降低。

Julia 社区在逐渐强大起来,其中开发了 XAM package 来处理 bam/sam 文件,感觉速度还是挺快的,在 windows 系统里可以使用。

2python 读取

这里有 1.8G 的 sam 文件, 我们先使用 python 的 for 循环逐行读取,对 reads 的长度进行计数保存到字典:

len_dict = {}
with open("../5.map-data/SRR7471228.sam",'r') as input:
    for line in input:
        if line.startswith('@'):
            continue
        fileds = line.split()
        if fileds[2] == '*':
            continue
        length = len(fileds[9])
        if length not in len_dict:
            len_dict[length] = 1
        else:
            len_dict[length] += 1

# sort dict
len_dict_sorted = list(len_dict.items())
len_dict_sorted.sort()

时间为 21.3 秒。

3Julia 读取

实现的代码比较类似:

notin 可以打出∉符号。

len_dict = Dict()
for line in eachline("../5.map-data/SRR7471228.sam")
    if startswith(line,"@")
        continue
    end
    fileds = split(line)
    if fileds[2] == '*'
        continue
    end
    # julia is 1-based
    seq_length = length(fileds[10])
    # println(seq_length)
    # break
    if seq_length ∉ keys(len_dict)
        len_dict[seq_length] = 1
    else
        len_dict[seq_length] += 1
    end
end

# sorted dict by read length
len_dictSorted = sort(len_dict)

可以看到运行了 55.9s,居然比 python 慢了两倍多!给我整不会了。

4XAM 读取

使用 XAM 包读取,依然是 Julia 语言:

using XAM

len_dict = Dict()

# Open a SAM file.
reader = open(SAM.Reader, "../5.map-data/SRR7471228.sam")

# Iterate over SAM records.
for record in reader
    # `record` is a SAM.Record object.
    if SAM.ismapped(record)
        # Print the mapped position.
        read_length = SAM.seqlength(record)
        if read_length ∉ keys(len_dict)
            len_dict[read_length] = 1
        else
            len_dict[read_length] += 1
        end
    end
end

# Close the SAM file.
close(reader)

# sorted dict by read length
len_dictSorted = sort(len_dict)

只花了 11.6s,比 python 快了 1 倍。

另外一种更加快速的读取方式:

using XAM

len_dict = Dict()

reader = open(SAM.Reader, "../5.map-data/SRR7471228.sam")
record = SAM.Record()
while !eof(reader)
    empty!(record)
    read!(reader, record)
    # do something
    if SAM.ismapped(record)
        # Print the mapped position.
        read_length = SAM.seqlength(record)
        if read_length ∉ keys(len_dict)
            len_dict[read_length] = 1
        else
            len_dict[read_length] += 1
        end
    end
end

# sorted dict by read length
len_dictSorted = sort(len_dict)