如何计算 read 的 covergae

Original JunJunLab 老俊俊的生信笔记 2022-08-17

收录于合集 #Ribo-seq 39个

流逝

1引言

我们在得到比对的 sam 或者 bam 文件后,可以直接加载到 IGV 里面进行可视化,查看比对情况。

或者使用 deeptools 转为 bigwig 文件进行可视化,在文章里比较常见的如下面这些 track 图:

基本上纵坐标是 read 的 coverage,也就是覆盖度。我们可以自己拿 sam/bam 文件自己计算。

read counts 和 coverage 的区别 (手绘):

2测试

思路:

主要思路为对 sam 文件比对上的 reads 每个碱基的位置进行计数,保存为字典,然后遇到相同的位置就叠加。
代码里还提供了 reads counts 计数的方法,可以选择。

定义函数:(julia语言编写)

# load package
using XAM

# define function
function CalculateRNACoverage(inputFile,outputFile;type="counts")
    # save in dict
    coverage_dict = Dict{String,Float64}()

    # open sam file
    reader = open(SAM.Reader,inputFile)
    record = SAM.Record()

    # loop
    while !eof(reader)
        empty!(record)
        read!(reader, record)
        # do something
        if SAM.ismapped(record)
            # tags
            refname = SAM.refname(record)
            align_pos = SAM.position(record)

            if type == "coverage" # get read coverage
                read_length = SAM.seqlength(record)
                # read right position
                End5 = align_pos + read_length - 1
                # ribo density
                for elem in range(align_pos,End5)
                    key = join([refname,elem],"|")
                    if !haskey(coverage_dict,key)
                        coverage_dict[key] = 1
                    else
                        coverage_dict[key] += 1
                    end
                end
            elseif type == "counts" # get read counts
                key = join([refname,align_pos],"|")
                if !haskey(coverage_dict,key)
                    coverage_dict[key] = 1
                else
                    coverage_dict[key] += 1
                end
            else
                println("error!")
                break
            end
        end
    end

    # output file
    outfile = open(outputFile,"w")

    # total densitys
    total_density = sum(values(coverage_dict))

    # sort keys
    for key in sort(collect(keys(coverage_dict)))
        id,align_pos = split(key,"|")
        raw_denisty = coverage_dict[key]

        # RPM normalization
        rpm = (raw_denisty/total_density)*1000000
        write(outfile,"$id\t$align_pos\t$raw_denisty\t$rpm\n")
    end
    close(outfile)
end

运行:

# calculate coverage
CalculateRNACoverage("test.sam","test.counts.txt",type="coverage")

# calculate read counts
CalculateRNACoverage("test.sam","test.counts.txt",type="counts")

输出结果为染色体,比对的位置,coverage,RPM:

10 100014676 1.0 0.007010635358178688
10 100014677 3.0 0.021031906074536068
10 100014678 4.0 0.028042541432714754
10 100014679 6.0 0.042063812149072136
10 100014680 8.0 0.05608508286542951
10 100014681 8.0 0.05608508286542951