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")
@@ -27,7 +27,7 @@ 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.
ismapped(record::Record)::Bool
Test if record is mapped.
isprimary(record::Record)::Bool
Test if record is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0.
refname(record::Record)::String
Get the reference sequence name of record.
position(record::Record)::Int
Get the 1-based leftmost mapping position of record.
rightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record.
isnextmapped(record::Record)::Bool
Test if the mate/next read of record is mapped.
nextrefname(record::Record)::String
Get the reference name of the mate/next read of record.
nextposition(record::Record)::Int
Get the position of the mate/next read of record.
mappingquality(record::Record)::UInt8
Get the mapping quality of record.
cigar(record::Record)::String
Get the CIGAR string of record.
alignment(record::Record)::BioAlignments.Alignment
Get the alignment of record.
alignlength(record::Record)::Int
Get the alignment length of record.
tempname(record::Record)::String
Get the query template name of record.
templength(record::Record)::Int
Get the template length of record.
sequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record.
sequence(::Type{String}, record::Record)::String
Get the segment sequence of record as String.
seqlength(record::Record)::Int
Get the sequence length of record.
quality(record::Record)::Vector{UInt8}
Get the Phred-scaled base quality of record.
quality(::Type{String}, record::Record)::String
Get the ASCII-encoded base quality of record.
auxdata(record::Record)::Dict{String,Any}
Get the auxiliary data (optional fields) of record.
The XAM package supports the following accessors for BAM.Record types.
flag(record::Record)::UInt16
Get the bitwise flag of record.
ismapped(record::Record)::Bool
Test if record is mapped.
isprimary(record::Record)::Bool
Test if record is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0.
refid(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
refname(record::Record)::String
Get the reference sequence name of record.
See also: BAM.refid
reflen(record::Record)::Int
Get the length of the reference sequence this record applies to.
position(record::Record)::Int
Get the 1-based leftmost mapping position of record.
rightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record.
isnextmapped(record::Record)::Bool
Test if the mate/next read of record is mapped.
nextrefid(record::Record)::Int
Get the next/mate reference sequence ID of record.
nextrefname(record::Record)::String
Get the reference name of the mate/next read of record.
nextposition(record::Record)::Int
Get the 1-based leftmost mapping position of the next/mate read of record.
mappingquality(record::Record)::UInt8
Get the mapping quality of record.
cigar(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.
alignment(record::Record)::BioAlignments.Alignment
Get the alignment of record.
alignlength(record::Record)::Int
Get the alignment length of record.
tempname(record::Record)::String
Get the query template name of record.
templength(record::Record)::Int
Get the template length of record.
sequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record.
seqlength(record::Record)::Int
Get the sequence length of record.
quality(record::Record)
Get the base quality of record.
auxdata(record::Record)::BAM.AuxData
Get the auxiliary data of record.
SAM 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.
ismapped(record::Record)::Bool
Test if record is mapped.
isprimary(record::Record)::Bool
Test if record is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0.
refname(record::Record)::String
Get the reference sequence name of record.
position(record::Record)::Int
Get the 1-based leftmost mapping position of record.
rightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record.
isnextmapped(record::Record)::Bool
Test if the mate/next read of record is mapped.
nextrefname(record::Record)::String
Get the reference name of the mate/next read of record.
nextposition(record::Record)::Int
Get the position of the mate/next read of record.
mappingquality(record::Record)::UInt8
Get the mapping quality of record.
cigar(record::Record)::String
Get the CIGAR string of record.
alignment(record::Record)::BioAlignments.Alignment
Get the alignment of record.
alignlength(record::Record)::Int
Get the alignment length of record.
tempname(record::Record)::String
Get the query template name of record.
templength(record::Record)::Int
Get the template length of record.
sequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record.
sequence(::Type{String}, record::Record)::String
Get the segment sequence of record as String.
seqlength(record::Record)::Int
Get the sequence length of record.
quality(record::Record)::Vector{UInt8}
Get the Phred-scaled base quality of record.
quality(::Type{String}, record::Record)::String
Get the ASCII-encoded base quality of record.
auxdata(record::Record)::Dict{String,Any}
Get the auxiliary data (optional fields) of record.
The XAM package supports the following accessors for BAM.Record types.
flag(record::Record)::UInt16
Get the bitwise flag of record.
ismapped(record::Record)::Bool
Test if record is mapped.
isprimary(record::Record)::Bool
Test if record is a primary line of the read.
This is equivalent to flag(record) & 0x900 == 0.
refid(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
refname(record::Record)::String
Get the reference sequence name of record.
See also: BAM.refid
reflen(record::Record)::Int
Get the length of the reference sequence this record applies to.
position(record::Record)::Int
Get the 1-based leftmost mapping position of record.
rightposition(record::Record)::Int
Get the 1-based rightmost mapping position of record.
isnextmapped(record::Record)::Bool
Test if the mate/next read of record is mapped.
nextrefid(record::Record)::Int
Get the next/mate reference sequence ID of record.
nextrefname(record::Record)::String
Get the reference name of the mate/next read of record.
nextposition(record::Record)::Int
Get the 1-based leftmost mapping position of the next/mate read of record.
mappingquality(record::Record)::UInt8
Get the mapping quality of record.
cigar(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.
alignment(record::Record)::BioAlignments.Alignment
Get the alignment of record.
alignlength(record::Record)::Int
Get the alignment length of record.
tempname(record::Record)::String
Get the query template name of record.
templength(record::Record)::Int
Get the template length of record.
sequence(record::Record)::BioSequences.LongDNASeq
Get the segment sequence of record.
seqlength(record::Record)::Int
Get the sequence length of record.
quality(record::Record)
Get the base quality of record.
auxdata(record::Record)::BAM.AuxData
Get the auxiliary data of record.
SAM 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")
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,7 +80,7 @@ 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
+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
@@ -104,4 +104,4 @@ SAM.Writer(IOStream(<file my-data.sam>))
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.Records or SAM.Records to file:
julia> write(bamw, rec) # Here rec is a `BAM.Record`
-330780
