How do I retrieve the sequence indices of an aligned region from a BioAlignments.jl alignment?

user_7210206 0
Viewed 76

I would like to access the indices of the alignment data that is returned from the Julia pairalign function in BioAlignments.jl to know where the alignment occurred in context of the original sequences.

using BioAlignments
using BioSequences
scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=-1);
my_alignment = pairalign(LocalAlignment(),dna"ATATTAGGTATTGATTATTGTACGCGGCCCGGC" , dna"TTGATTATTGT", scoremodel)
alignment(my_alignment)

For example, a script like this will output an alignment object from which I can access the score via score() function. However I wish to know where, in the original sequences that I supplied as inputs, the alignment occurred and know how to call the variable storing this index. Could not seem to find this anywhere in the documentation.

julia bioinformatics
1 Answers

While I do not use those libraries one of the approach to such problems in Julia is to dump such object and in result one can find out the location of needed information.

julia> dump(alignment(my_alignment))
PairwiseAlignment{LongSequence{DNAAlphabet{4}},LongSequence{DNAAlphabet{4}}}
  a: AlignedSequence{LongSequence{DNAAlphabet{4}}}
    seq: LongSequence{DNAAlphabet{4}}
      data: Array{UInt64}((3,)) UInt64[0x8814881844188181, 0x4422244242184881, 0x0000000000000002]
      part: UnitRange{Int64}
        start: Int64 1
        stop: Int64 33
      shared: Bool false
    aln: Alignment
      anchors: Array{AlignmentAnchor}((2,))
        1: AlignmentAnchor
          seqpos: Int64 10
          refpos: Int64 0
          op: Operation OP_START
        2: AlignmentAnchor
          seqpos: Int64 21
          refpos: Int64 11
          op: Operation OP_SEQ_MATCH
      firstref: Int64 1
      lastref: Int64 11
  b: LongSequence{DNAAlphabet{4}}
    data: Array{UInt64}((1,)) UInt64[0x0000084881881488]
    part: UnitRange{Int64}
      start: Int64 1
      stop: Int64 11
    shared: Bool false

And now you can see where is the information that you need:

julia> alignment(my_alignment).a.aln.anchors
2-element Array{AlignmentAnchor,1}:
 AlignmentAnchor(10, 0, '0')
 AlignmentAnchor(21, 11, '=')

The disadvantage of this approach is that the data structures are usually not the part of the library API and could change over the time with new package realeases.

user_99577100
Related