The XAM
package offers high-performance tools for SAM and BAM file formats, which are the most popular file formats.
If you have questions about the SAM and BAM formats or any of the terminology used when discussing these formats, see the published specification, which is maintained by the samtools group.
A very very simple SAM file looks like the following:
@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r002 0 ref 9 30 3S6M1P1I4M * 0 0 AAAAGATAAGGATA *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
-r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1
Where the first two lines are part of the "header", and the following lines are "records". Each record describes how a read aligns to some reference sequence. Sometimes one record describes one read, but there are other cases like chimeric reads and split alignments, where multiple records apply to one read. In the example above, r003
is a chimeric read, and r004
is a split alignment, and r001
are mate pair reads. Again, we refer you to the official specification for more details.
A BAM file stores this same information but in a binary and compressible format that does not make for pretty printing here!
A typical script iterating over all records in a file looks like below:
using XAM
+r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1
Where the first two lines are part of the "header", and the following lines are "records". Each record describes how a read aligns to some reference sequence. Sometimes one record describes one read, but there are other cases like chimeric reads and split alignments, where multiple records apply to one read. In the example above, r003
is a chimeric read, and r004
is a split alignment, and r001
are mate pair reads. Again, we refer you to the official specification for more details.
A BAM file stores this same information but in a binary and compressible format that does not make for pretty printing here!
A typical script iterating over all records in a file looks like below:
using XAM
# Open a BAM file.
reader = open(BAM.Reader, "data.bam")
@@ -21,13 +21,13 @@ for record in reader
end
# Close the BAM file.
-close(reader)
The size of a BAM file is often extremely large. The iterator interface demonstrated above allocates an object for each record and that may be a bottleneck of reading data from a BAM file. In-place reading reuses a pre-allocated object for every record and less memory allocation happens in reading:
reader = open(BAM.Reader, "data.bam")
+close(reader)
The size of a BAM file is often extremely large. The iterator interface demonstrated above allocates an object for each record and that may be a bottleneck of reading data from a BAM file. In-place reading reuses a pre-allocated object for every record and less memory allocation happens in reading:
reader = open(BAM.Reader, "data.bam")
record = BAM.Record()
while !eof(reader)
empty!(record)
read!(reader, record)
# do something
-end
Both SAM.Reader
and BAM.Reader
implement the header
function, which returns a SAM.Header
object. To extract certain information out of the headers, you can use the find
method on the header to extract information according to SAM/BAM tag. Again we refer you to the specification for full details of all the different tags that can occur in headers, and what they mean.
Below is an example of extracting all the info about the reference sequences from the BAM header. In SAM/BAM, any description of a reference sequence is stored in the header, under a tag denoted SQ
(think reference SeQuence
!).
julia> reader = open(SAM.Reader, "data.sam");
+end
Both SAM.Reader
and BAM.Reader
implement the header
function, which returns a SAM.Header
object. To extract certain information out of the headers, you can use the find
method on the header to extract information according to SAM/BAM tag. Again we refer you to the specification for full details of all the different tags that can occur in headers, and what they mean.
Below is an example of extracting all the info about the reference sequences from the BAM header. In SAM/BAM, any description of a reference sequence is stored in the header, under a tag denoted SQ
(think reference SeQuence
!).
julia> reader = open(SAM.Reader, "data.sam");
julia> find(header(reader), "SQ")
7-element Array{Bio.Align.SAM.MetaInfo,1}:
@@ -52,17 +52,17 @@ julia> find(header(reader), "SQ")
Bio.Align.SAM.MetaInfo:
tag: SQ
value: SN=mitochondria LN=366924
-
In the above we can see there were 7 sequences in the reference: 5 chromosomes, one chloroplast sequence, and one mitochondrial sequence.
The XAM
package supports the following accessors for SAM.Record
types.
flag(record::Record)::UInt16
Get the bitwise flag of record
.
sourceismapped(record::Record)::Bool
Test if record
is mapped.
sourceisprimary(record::Record)::Bool
Test if record
is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0
.
sourcerefname(record::Record)::String
Get the reference sequence name of record
.
sourceposition(record::Record)::Int
Get the 1-based leftmost mapping position of record
.
sourcerightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record
.
sourceisnextmapped(record::Record)::Bool
Test if the mate/next read of record
is mapped.
sourcenextrefname(record::Record)::String
Get the reference name of the mate/next read of record
.
sourcenextposition(record::Record)::Int
Get the position of the mate/next read of record
.
sourcemappingquality(record::Record)::UInt8
Get the mapping quality of record
.
sourcecigar(record::Record)::String
Get the CIGAR string of record
.
sourcealignment(record::Record)::BioAlignments.Alignment
Get the alignment of record
.
sourcealignlength(record::Record)::Int
Get the alignment length of record
.
sourcetempname(record::Record)::String
Get the query template name of record
.
sourcetemplength(record::Record)::Int
Get the template length of record
.
sourcesequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record
.
sourcesequence(::Type{String}, record::Record)::String
Get the segment sequence of record
as String
.
sourceseqlength(record::Record)::Int
Get the sequence length of record
.
sourcequality(record::Record)::Vector{UInt8}
Get the Phred-scaled base quality of record
.
sourcequality(::Type{String}, record::Record)::String
Get the ASCII-encoded base quality of record
.
sourceauxdata(record::Record)::Dict{String,Any}
Get the auxiliary data (optional fields) of record
.
sourceThe XAM
package supports the following accessors for BAM.Record
types.
flag(record::Record)::UInt16
Get the bitwise flag of record
.
sourceismapped(record::Record)::Bool
Test if record
is mapped.
sourceisprimary(record::Record)::Bool
Test if record
is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0
.
sourcerefid(record::Record)::Int
Get the reference sequence ID of record
.
The ID is 1-based (i.e. the first sequence is 1) and is 0 for a record without a mapping position.
See also: BAM.rname
sourcerefname(record::Record)::String
Get the reference sequence name of record
.
See also: BAM.refid
sourcereflen(record::Record)::Int
Get the length of the reference sequence this record applies to.
sourceposition(record::Record)::Int
Get the 1-based leftmost mapping position of record
.
sourcerightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record
.
sourceisnextmapped(record::Record)::Bool
Test if the mate/next read of record
is mapped.
sourcenextrefid(record::Record)::Int
Get the next/mate reference sequence ID of record
.
sourcenextrefname(record::Record)::String
Get the reference name of the mate/next read of record
.
sourcenextposition(record::Record)::Int
Get the 1-based leftmost mapping position of the next/mate read of record
.
sourcemappingquality(record::Record)::UInt8
Get the mapping quality of record
.
sourcecigar(record::Record)::String
Get the CIGAR string of record
.
Note that in the BAM specification, the field called cigar
typically stores the cigar string of the record. However, this is not always true, sometimes the true cigar is very long, and due to some constraints of the BAM format, the actual cigar string is stored in an extra tag: CG:B,I
, and the cigar
field stores a pseudo-cigar string.
Calling this method with checkCG
set to true
(default) this method will always yield the true cigar string, because this is probably what you want the vast majority of the time.
If you have a record that stores the true cigar in a CG:B,I
tag, but you still want to access the pseudo-cigar that is stored in the cigar
field of the BAM record, then you can set checkCG to false
.
See also BAM.cigar_rle
.
sourcealignment(record::Record)::BioAlignments.Alignment
Get the alignment of record
.
sourcealignlength(record::Record)::Int
Get the alignment length of record
.
sourcetempname(record::Record)::String
Get the query template name of record
.
sourcetemplength(record::Record)::Int
Get the template length of record
.
sourcesequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record
.
sourceseqlength(record::Record)::Int
Get the sequence length of record
.
sourcequality(record::Record)
Get the base quality of record
.
sourceauxdata(record::Record)::BAM.AuxData
Get the auxiliary data of record
.
sourceSAM and BAM records support the storing of optional data fields associated with tags.
Tagged auxiliary data follows a format of TAG:TYPE:VALUE
. TAG
is a two-letter string, and each tag can only appear once per record. TYPE
is a single case-sensetive letter which defined the format of VALUE
.
Type | Description |
---|
'A' | Printable character |
'i' | Signed integer |
'f' | Single-precision floating number |
'Z' | Printable string, including space |
'H' | Byte array in Hex format |
'B' | Integer of numeric array |
For more information about these tags and their types we refer you to the SAM/BAM specification and the additional optional fields specification document.
There are some tags that are reserved, predefined standard tags, for specific uses.
To access optional fields stored in tags, you use getindex
indexing syntax on the record object. Note that accessing optional tag fields will result in type instability in Julia. This is because the type of the optional data is not known until run-time, as the tag is being read. This can have a significant impact on performance. To limit this, if the user knows the type of a value in advance, specifying it as a type annotation will alleviate the problem:
Below is an example of looping over records in a bam file and using indexing syntax to get the data stored in the "NM" tag. Note the UInt8
type assertion to alleviate type instability.
for record in open(BAM.Reader, "data.bam")
+
In the above we can see there were 7 sequences in the reference: 5 chromosomes, one chloroplast sequence, and one mitochondrial sequence.
The XAM
package supports the following accessors for SAM.Record
types.
flag(record::Record)::UInt16
Get the bitwise flag of record
.
sourceismapped(record::Record)::Bool
Test if record
is mapped.
sourceisprimary(record::Record)::Bool
Test if record
is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0
.
sourcerefname(record::Record)::String
Get the reference sequence name of record
.
sourceposition(record::Record)::Int
Get the 1-based leftmost mapping position of record
.
sourcerightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record
.
sourceisnextmapped(record::Record)::Bool
Test if the mate/next read of record
is mapped.
sourcenextrefname(record::Record)::String
Get the reference name of the mate/next read of record
.
sourcenextposition(record::Record)::Int
Get the position of the mate/next read of record
.
sourcemappingquality(record::Record)::UInt8
Get the mapping quality of record
.
sourcecigar(record::Record)::String
Get the CIGAR string of record
.
sourcealignment(record::Record)::BioAlignments.Alignment
Get the alignment of record
.
sourcealignlength(record::Record)::Int
Get the alignment length of record
.
sourcetempname(record::Record)::String
Get the query template name of record
.
sourcetemplength(record::Record)::Int
Get the template length of record
.
sourcesequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record
.
sourcesequence(::Type{String}, record::Record)::String
Get the segment sequence of record
as String
.
sourceseqlength(record::Record)::Int
Get the sequence length of record
.
sourcequality(record::Record)::Vector{UInt8}
Get the Phred-scaled base quality of record
.
sourcequality(::Type{String}, record::Record)::String
Get the ASCII-encoded base quality of record
.
sourceauxdata(record::Record)::Dict{String,Any}
Get the auxiliary data (optional fields) of record
.
sourceThe XAM
package supports the following accessors for BAM.Record
types.
flag(record::Record)::UInt16
Get the bitwise flag of record
.
sourceismapped(record::Record)::Bool
Test if record
is mapped.
sourceisprimary(record::Record)::Bool
Test if record
is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0
.
sourcerefid(record::Record)::Int
Get the reference sequence ID of record
.
The ID is 1-based (i.e. the first sequence is 1) and is 0 for a record without a mapping position.
See also: BAM.rname
sourcerefname(record::Record)::String
Get the reference sequence name of record
.
See also: BAM.refid
sourcereflen(record::Record)::Int
Get the length of the reference sequence this record applies to.
sourceposition(record::Record)::Int
Get the 1-based leftmost mapping position of record
.
sourcerightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record
.
sourceisnextmapped(record::Record)::Bool
Test if the mate/next read of record
is mapped.
sourcenextrefid(record::Record)::Int
Get the next/mate reference sequence ID of record
.
sourcenextrefname(record::Record)::String
Get the reference name of the mate/next read of record
.
sourcenextposition(record::Record)::Int
Get the 1-based leftmost mapping position of the next/mate read of record
.
sourcemappingquality(record::Record)::UInt8
Get the mapping quality of record
.
sourcecigar(record::Record)::String
Get the CIGAR string of record
.
Note that in the BAM specification, the field called cigar
typically stores the cigar string of the record. However, this is not always true, sometimes the true cigar is very long, and due to some constraints of the BAM format, the actual cigar string is stored in an extra tag: CG:B,I
, and the cigar
field stores a pseudo-cigar string.
Calling this method with checkCG
set to true
(default) this method will always yield the true cigar string, because this is probably what you want the vast majority of the time.
If you have a record that stores the true cigar in a CG:B,I
tag, but you still want to access the pseudo-cigar that is stored in the cigar
field of the BAM record, then you can set checkCG to false
.
See also BAM.cigar_rle
.
sourcealignment(record::Record)::BioAlignments.Alignment
Get the alignment of record
.
sourcealignlength(record::Record)::Int
Get the alignment length of record
.
sourcetempname(record::Record)::String
Get the query template name of record
.
sourcetemplength(record::Record)::Int
Get the template length of record
.
sourcesequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record
.
sourceseqlength(record::Record)::Int
Get the sequence length of record
.
sourcequality(record::Record)
Get the base quality of record
.
sourceauxdata(record::Record)::BAM.AuxData
Get the auxiliary data of record
.
sourceSAM and BAM records support the storing of optional data fields associated with tags.
Tagged auxiliary data follows a format of TAG:TYPE:VALUE
. TAG
is a two-letter string, and each tag can only appear once per record. TYPE
is a single case-sensetive letter which defined the format of VALUE
.
Type | Description |
---|
'A' | Printable character |
'i' | Signed integer |
'f' | Single-precision floating number |
'Z' | Printable string, including space |
'H' | Byte array in Hex format |
'B' | Integer of numeric array |
For more information about these tags and their types we refer you to the SAM/BAM specification and the additional optional fields specification document.
There are some tags that are reserved, predefined standard tags, for specific uses.
To access optional fields stored in tags, you use getindex
indexing syntax on the record object. Note that accessing optional tag fields will result in type instability in Julia. This is because the type of the optional data is not known until run-time, as the tag is being read. This can have a significant impact on performance. To limit this, if the user knows the type of a value in advance, specifying it as a type annotation will alleviate the problem:
Below is an example of looping over records in a bam file and using indexing syntax to get the data stored in the "NM" tag. Note the UInt8
type assertion to alleviate type instability.
for record in open(BAM.Reader, "data.bam")
nm = record["NM"]::UInt8
# do something
-end
The XAM
package supports the BAI index to fetch records in a specific range from a BAM file. Samtools provides index
subcommand to create an index file (.bai) from a sorted BAM file.
$ samtools index -b SRR1238088.sort.bam
+end
The XAM
package supports the BAI index to fetch records in a specific range from a BAM file. Samtools provides index
subcommand to create an index file (.bai) from a sorted BAM file.
$ samtools index -b SRR1238088.sort.bam
$ ls SRR1238088.sort.bam*
-SRR1238088.sort.bam SRR1238088.sort.bam.bai
The method eachoverlap(reader, chrom, range)
returns an iterator of BAM records overlapping the query interval:
reader = open(BAM.Reader, "SRR1238088.sort.bam", index="SRR1238088.sort.bam.bai")
+SRR1238088.sort.bam SRR1238088.sort.bam.bai
The method eachoverlap(reader, chrom, range)
returns an iterator of BAM records overlapping the query interval:
reader = open(BAM.Reader, "SRR1238088.sort.bam", index="SRR1238088.sort.bam.bai")
for record in eachoverlap(reader, "Chr2", 10000:11000)
# `record` is a BAM.Record object
# ...
end
-close(reader)
The eachoverlap
method also accepts the Interval
type defined in GenomicFeatures.jl.
This allows you to do things like first read in the genomic features from a GFF3 file, and then for each feature, iterate over all the BAM records that overlap with that feature.
using GenomicFeatures
+close(reader)
The eachoverlap
method also accepts the Interval
type defined in GenomicFeatures.jl.
This allows you to do things like first read in the genomic features from a GFF3 file, and then for each feature, iterate over all the BAM records that overlap with that feature.
using GenomicFeatures
using GFF3
using XAM
@@ -80,8 +80,8 @@ for feature in features
# ...
end
end
-close(reader)
In order to write a BAM or SAM file, you must first create a SAM.Header
.
A SAM.Header
is constructed from a vector of SAM.MetaInfo
objects.
For example, to create the following simple header:
@HD VN:1.6 SO:coordinate
-@SQ SN:ref LN:45
julia> a = SAM.MetaInfo("HD", ["VN" => 1.6, "SO" => "coordinate"])
+close(reader)
In order to write a BAM or SAM file, you must first create a SAM.Header
.
A SAM.Header
is constructed from a vector of SAM.MetaInfo
objects.
For example, to create the following simple header:
@HD VN:1.6 SO:coordinate
+@SQ SN:ref LN:45
julia> a = SAM.MetaInfo("HD", ["VN" => 1.6, "SO" => "coordinate"])
SAM.MetaInfo:
tag: HD
value: VN=1.6 SO=coordinate
@@ -97,11 +97,11 @@ SAM.Header(SAM.MetaInfo[SAM.MetaInfo:
value: VN=1.6 SO=coordinate, SAM.MetaInfo:
tag: SQ
value: SN=ref LN=45])
-
Then to create the writer for a SAM file, construct a SAM.Writer
using the header and an IO
type:
julia> samw = SAM.Writer(open("my-data.sam", "w"), h)
+
Then to create the writer for a SAM file, construct a SAM.Writer
using the header and an IO
type:
julia> samw = SAM.Writer(open("my-data.sam", "w"), h)
SAM.Writer(IOStream(<file my-data.sam>))
-
To make a BAM Writer is slightly different, as you need to use a specific stream type from the https://github.com/BioJulia/BGZFStreams.jl package:
julia> using BGZFStreams
+
To make a BAM Writer is slightly different, as you need to use a specific stream type from the https://github.com/BioJulia/BGZFStreams.jl package:
julia> using BGZFStreams
julia> bamw = BAM.Writer(BGZFStream(open("my-data.bam", "w"), "w"))
BAM.Writer(BGZFStreams.BGZFStream{IOStream}(<mode=write>))
-
Once you have a BAM or SAM writer, you can use the write
method to write BAM.Record
s or SAM.Record
s to file:
julia> write(bamw, rec) # Here rec is a `BAM.Record`
-330780