Increment dependency compatibilities

Project.toml

@@ -17,9 +17,9 @@ TranscodingStreams = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 [compat]
 Automa = "0.7, 0.8"
 BGZFStreams = "0.3"
-BioAlignments = "2"
+BioAlignments = "3"
 BioGenerics = "0.1"
-BioSequences = "2.0.4"
+BioSequences = "3"
 GenomicFeatures = "2"
 Indexes = "0.1"
 TranscodingStreams = "0.6, 0.7, 0.8, 0.9"

src/bam/record.jl

@@ -91,7 +91,7 @@ function Base.empty!(record::Record)
     record.tlen       = 0
 
     #Note: data will be overwritten and indexed using data_size.
-    
+
     return record
 end
 
@@ -436,7 +436,8 @@ function alignment(record::Record)::BioAlignments.Alignment
     end
     seqpos = 0
     refpos = position(record) - 1
-    anchors = [BioAlignments.AlignmentAnchor(seqpos, refpos, BioAlignments.OP_START)]
+    alnpos = 0
+    anchors = [BioAlignments.AlignmentAnchor(seqpos, refpos, alnpos, BioAlignments.OP_START)]
     for (op, len) in zip(cigar_rle(record)...)
         if BioAlignments.ismatchop(op)
             seqpos += len
@@ -448,7 +449,8 @@ function alignment(record::Record)::BioAlignments.Alignment
         else
             error("operation $(op) is not supported")
         end
-        push!(anchors, BioAlignments.AlignmentAnchor(seqpos, refpos, op))
+        alnpos += len
+        push!(anchors, BioAlignments.AlignmentAnchor(seqpos, refpos, alnpos, op))
     end
     return BioAlignments.Alignment(anchors)
 end
@@ -504,12 +506,7 @@ function hastemplength(record::Record)
     return isfilled(record)
 end
 
-"""
-    sequence(record::Record)::BioSequences.LongDNASeq
-
-Get the segment sequence of `record`.
-"""
-function sequence(record::Record)
+function sequence(::Type{S}, record::Record) where S <: BioSequences.LongSequence #TODO: add optional part retrieval.
     checkfilled(record)
     seqlen = seqlength(record)
     if seqlen == 0
@@ -522,7 +519,16 @@ function sequence(record::Record)
         x = unsafe_load(src, i)
         data[i] = (x & 0x0f0f0f0f0f0f0f0f) << 4 | (x & 0xf0f0f0f0f0f0f0f0) >> 4
     end
-    return BioSequences.LongDNASeq(data, 1:seqlen, false)
+    return S(data, UInt64(seqlen))
+end
+
+"""
+    sequence(record::Record)::BioSequences.LongDNA{4}
+
+Get the segment sequence of `record`.
+"""
+function sequence(record::Record) #TODO: add optional part retrieval.
+    return sequence(BioSequences.LongDNA{4}, record)
 end
 
 function hassequence(record::Record)

src/sam/record.jl

@@ -387,34 +387,43 @@ function hastemplength(record::Record)
 end
 
 """
-    sequence(record::Record)::BioSequences.LongDNASeq
-
-Get the segment sequence of `record`.
-"""
-function sequence(record::Record)
-    checkfilled(record)
-    if ismissing(record, record.seq)
-        # missingerror(:sequence)
-        return nothing
-    end
-    seqlen = length(record.seq)
-    ret = BioSequences.LongDNASeq(seqlen)
-    copyto!(ret, 1, record.data, first(record.seq), seqlen)
-    return ret
+    sequence(::Type{S}, record::Record, [part::UnitRange{Int}])::S
+
+Get the sequence of `record`.
+`S` can be either a subtype of `BioSequences.BioSequence` or `String`.
+If `part` argument is given, it returns the specified part of the sequence.
+!!! note
+    This method makes a new sequence object every time.
+    If you have a sequence already and want to fill it with the sequence
+    data contained in a fasta record, you can use `Base.copyto!`.
+"""
+function sequence(::Type{S}, record::Record, part::UnitRange{Int}=1:lastindex(record.seq)) where S <: BioSequences.LongSequence
+     checkfilled(record)
+     if ismissing(record, record.seq)
+         # missingerror(:sequence)
+         return nothing
+     end
+     seqpart = record.seq[part]
+     return S(@view(record.data[seqpart]))
+ end
+
+ function sequence(::Type{String}, record::Record, part::UnitRange{Int}=1:lastindex(record.seq))
+     checkfilled(record)
+     return String(record.data[record.seq[part]])
+ end
+
+ """
+     sequence(record::Record)::BioSequences.LongDNA{4}
+
+ Get the segment sequence of `record`.
+ """
+function sequence(record::Record, part::UnitRange{Int}=1:lastindex(record.seq))
+    return sequence(BioSequences.LongDNA{4}, record, part)
 end
 
 function hassequence(record::Record)
     return isfilled(record) && !ismissing(record, record.seq)
 end
-"""
-    sequence(::Type{String}, record::Record)::String
-
-Get the segment sequence of `record` as `String`.
-"""
-function sequence(::Type{String}, record::Record)::String
-    checkfilled(record)
-    return String(record.data[record.seq])
-end
 
 """
     seqlength(record::Record)::Int

test/test_bam.jl

@@ -79,10 +79,10 @@
         @test BAM.flag(record) === UInt16(16)
         @test BAM.cigar(record) == "27M1D73M"
         @test BAM.alignment(record) == Alignment([
-            AlignmentAnchor(  0,   1, OP_START),
-            AlignmentAnchor( 27,  28, OP_MATCH),
-            AlignmentAnchor( 27,  29, OP_DELETE),
-            AlignmentAnchor(100, 102, OP_MATCH)])
+            AlignmentAnchor(  0,   1,   0, OP_START),
+            AlignmentAnchor( 27,  28,  27, OP_MATCH),
+            AlignmentAnchor( 27,  29,  28, OP_DELETE),
+            AlignmentAnchor(100, 102, 101, OP_MATCH)])
         @test record["XG"] == 1
         @test record["XM"] == 5
         @test record["XN"] == 0

test/test_sam.jl

@@ -104,10 +104,10 @@
         @test SAM.flag(record) == 16
         @test SAM.cigar(record) == "27M1D73M"
         @test SAM.alignment(record) == Alignment([
-            AlignmentAnchor(  0,   1, OP_START),
-            AlignmentAnchor( 27,  28, OP_MATCH),
-            AlignmentAnchor( 27,  29, OP_DELETE),
-            AlignmentAnchor(100, 102, OP_MATCH)])
+            AlignmentAnchor(  0,   1,   0, OP_START),
+            AlignmentAnchor( 27,  28,  27, OP_MATCH),
+            AlignmentAnchor( 27,  29,  28, OP_DELETE),
+            AlignmentAnchor(100, 102, 101, OP_MATCH)])
         @test record["XG"] == 1
         @test record["XM"] == 5
         @test record["XN"] == 0