# Julia package to calculate coverage and depth from sam/bam files
**URL:** https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888
**Category:** Biology, Health, and Medicine
**Tags:** question, package, biology, io
**Created:** [November 22, 2021, 9:04am UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888 "2021-11-22T09:04:33Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Lamma](https://discourse.julialang.org/u/Lamma)
#### Post date: [November 22, 2021, 9:04am UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888/1 "2021-11-22T09:04:33Z")
Hi there,
I am wanting to generate coverage and depth statistic for some sam/bam files I have and was wondering if Julia has a package that can do this for me?
---
### Author:  [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)
#### Post date: [November 22, 2021, 9:26am UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888/2 "2021-11-22T09:26:05Z")
Disclaimer: I have no idea what “sam/bam files” are.
A google search for “sam/bam files jl” led me directly to [https://biojulia.net/BioAlignments.jl/latest/hts-files.html](https://biojulia.net/BioAlignments.jl/latest/hts-files.html) which seems to suggest that [BioAlignments.jl](https://github.com/BioJulia/BioAlignments.jl) might be what you’re looking for 🙂
Alternatively, [https://github.com/BioJulia/XAM.jl](https://github.com/BioJulia/XAM.jl) might be worth a look.
More generally, you should probably check out [https://biojulia.net/](https://biojulia.net/).
---
### Author:  [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)
#### Post date: [November 22, 2021, 11:08am UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888/3 "2021-11-22T11:08:22Z")
I think there’s a function that computes coverage from an array of alignements somewhere, but here’s how you can do it from scratch, it’s a little bit complicated but it gives you complete control on how the coverage is computed (you can filter out reads, softclips, etc) :
```julia
using BioAlignments, XAM
import XAM.BAM
function get_coverage(bam_file, chr)
# use dictionnaries to store coverage
coverage = Dict(chr => Dict{Int,Int}() for chr in chr)
reader = open(BAM.Reader, bam_file)
record = BAM.Record()
while !eof(reader)
read!(reader, record)
!BAM.ismapped(record) && continue
chr = BAM.refname(record)
aln = BAM.alignment(record)
# anchor is one base before sequence, so +1
seqrange = (first(aln.anchors).seqpos + 1):last(aln.anchors).seqpos
#loop through the sequence coordinates
for i in seqrange
# get corresponding position in the reference genome
refpos, OP = seq2ref(aln, i)
# count only matches (skip softclip, deletions, ..)
if OP == OP_MATCH
coverage[chr][refpos] = get(coverage[chr], refpos, 0) + 1
end
end
end
close(reader)
coverage
end
bam_file = "..."
chr = vcat(string.(1:22), "X", "Y")
cov = get_coverage(bam_file, chr)
sum(values(cov["20"])) # total coverage on chr. 20
```
This helps with alignments and anchors : [https://biojulia.net/BioAlignments.jl/latest/alignments.html#Overview-1](https://biojulia.net/BioAlignments.jl/latest/alignments.html#Overview-1)
---
### Author:  [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)
#### Post date: [November 22, 2021, 8:40pm UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888/4 "2021-11-22T20:40:17Z")
Just note that for paired-end reads that might counts alignments twice, but not 100% sure.
---
### Author:  [@Lamma](https://discourse.julialang.org/u/Lamma)
#### Post date: [November 23, 2021, 7:42am UTC](https://discourse.julialang.org/t/julia-package-to-calculate-coverage-and-depth-from-sam-bam-files/71888/5 "2021-11-23T07:42:27Z")
> [@jonathanBieler](#):
>
> hink there’s a function that computes coverage from an array of alignements somewhere, but here’s how you can do it from scratch, it’s a little bit complicated but it gives you complete control on how the coverage is computed (you can filter out reads, softclips, etc) :
Many thanks, I will give this a go and see how it compares to other tools!