Nanopore signal to sequence alignement.

Fields

• read::XAM.XAMRecord: Reference to the read.
• sig_index::Union{SignalIndex, Missing}: Signal interval to k-mer interval index.
• seq_index::Union{SequenceIndex, Missing}: K-mer interval to signal interval index.
• reverse::Bool: True if the molecule is sequenced from 3′ to 5′ end (i.e. in dRNA-seq).
• kmer_size::Int: Size of the pore k-mer
• kmer_center::Int: Which position to use as a representative for a k-mer

Implements: read, sequence, isreverse, seq_ends, sig_ends, steps

SignalAlignment(
    read::XAM.XAMRecord;
    kmer_size=missing,
    kmer_center=missing,
    pore=missing,
    index=:both
)

Get a SignalAlignment object for the given read. The function will try to parse the tags in the BAM record in order to find alignment data produced by a resquiggler.

Currently supported resquiggling formats:

• F5c style "si"/"ss" tags. These can be produced by F5c or Uncalled4.

The produced alignment maps k-mer intervals (a [start, end] pair of sequence positions) to raw signal intervals. The kmer_size should match the correct number for the pore used for sequencing. Since the library also provides for convenience functions that use positions instead of k-mer intervals, the kmer_center parameter is used for controlling which position of the k-mer is considered the "representative one". The parameters kmer_size and kmer_center can be set manually, or alternatively the pore parameter can be used to get preset values.

Arguments

• read::XAM.XAMRecord: a XAM.XAMRecord with the read for which to get the signal alignment.
• kmer_size: How many neighbouring nucleotides influence a given measurement.
• kmer_center: Which position of the k-mer exerts the largest influence. This will be used to convert k-mer ranges to single positions.
• pore: Can be used instead of setting kmer_size and kmer_center manually. In this case, preset values for kmer_size and kmer_center will be used, sourced from https://doi.org/10.1038/s41592-025-02631-4. Supported values: :r9_4, :r10_4, :RNA001, :RNA002, RNA004.
• index: :sequence, :signal, or :both (default: :both) Used to control what indices to construct for accessing the alignment. Using only one index can limit the available functionality or degrade performance, but may be useful to reduce memory footprint.

read(sig_aln::SignalAlignment)

Return the SAM/BAM record for the read.

Arguments

• sig_aln::SignalAlignment: signal alignment.

sequence(sig_aln::SignalAlignment)

Returns the nucleotide sequence of the read.

Arguments

• sig_aln::SignalAlignment: signal alignment.

isreverse(sig_aln::SignalAlignment)

Returns true if the molecule is sequenced from 3′ to 5′ end (i.e. in dRNA-seq).

Arguments

• sig_aln::SignalAlignment: signal alignment.

seq_ends(sig_aln::SignalAlignment)

Returns the (first, last) sequence positions in the alignment.

Arguments

• sig_aln::SignalAlignment: signal alignment.

sig_ends(sig_aln::SignalAlignment)

Returns the (first, last) signal positions in the alignment.

Arguments

• sig_aln::SignalAlignment: signal alignment.

steps(sig_aln::SignalAlignment)

Returns a vector of ((sig_start, sig_end), (seq_start, seq_end)) tuples representing the signal-to-sequence alignment with gaps in both signal and sequence.

Arguments

• sig_aln::SignalAlignment: signal alignment.

fast_alignment(read::XAM.XAMRecord)

Parse the SAM/BAM tags to get a minimalistic alignment in a fast manner.

Arguments

• read::XAM.XAMRecord: SAM/BAM read with si/ss tags.

Returns

Vector with ((sig_start, sig_end), (seq_start, seq_end)) pairs. Gaps in the signal or sequence will be represented with missing.

Example

julia> fast_alignment(read)
127-element Vector{Tuple{Union{Missing, Tuple{Integer, Integer}}, Union{Missing, Tuple{Integer, Integer}}}}:
 ((10618, 10622), (424, 424))
 ((10623, 10697), (423, 423))
 ((10698, 10715), (422, 422))
 ((10716, 10737), (421, 421))
 ((10738, 10755), (420, 420))
 ⋮
 ((15575, 15618), (19, 19))
 ((15619, 15659), (18, 18))
 ((15660, 15666), (17, 17))
 ((15667, 15676), (16, 16))
 ((15677, 15696), (15, 15))

slice(
    sig_aln::SignalAlignment;
    from::Integer=1,
    to::Integer=typemax(Int64),
    by::Symbol=:sequence
)

Slice the SignalAlignment to a part of the alignment. Slicing can be done by :sequence or :signal:

• :sequence - obtain a part of the alignment within a specific interval of sequence positions.
• :signal - obtain a part of the alignment within a specific interval of the raw signal.

Note: All k-mer intervals or signals overlapping the given interval will be included. E.g. an interval of 1 base pair can return 9 k-mers that include that position in RNA004, where the kmer_size is 9.

Arguments

• sig_aln::SignalAlignment: signal alignment.
• from::Integer: start index of the slice (inclusive).
• to::Integer: end index of the slice (inclusive).
• by::Symbol: slice by :sequence or :signal.

kmer_to_sig(
    sig_aln::SignalAlignment;
    from::Integer=1,
    to::Integer=typemax(Int64),
    by::Symbol=:sequence
)

Returns a generator for iterating over sequence k-mers to signal interval pairs. Gaps in the sequence will be skipped. The function allows limiting the results to specific sequence position or signal intervals. All intervals used for limiting the results are closed (inclusive) and any overlapping k-mer or signal interval will be included. Results will be sorted in increasing order of sequence positions.

Arguments

• sig_aln::SignalAlignment: signal alignment.
• from::Integer: start index of the slice (inclusive).
• to::Integer: end index of the slice (inclusive).
• by::Symbol: slice by :sequence or :signal.

Examples

These are based on the examples from F5c's documentation:

• DNA: https://hasindu2008.github.io/f5c/docs/output#dna-example
• dRNA: https://hasindu2008.github.io/f5c/docs/output#direct-rna-example

(also shown in the Introduction above) Note: f5c uses 0-based indices and shows the k-mer start. We use 1-based indices and use the full k-mer interval. For this example, let's assume kmer_size=6 and kmer_center=3 (the third base of the 6-mer influences the signal the most).

DNA examples:

julia> kmer_to_sig(dna_aln; from=2, to=4) |> collect
4-element Vector{Pair{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (1, 6) => (3, 4)
 (2, 7) => (5, 7)
 (3, 8) => (8, 9)
 (4, 9) => (12, 12)


julia> kmer_to_sig(dna_aln; from=2, to=5, by=:signal) |> collect
2-element Vector{Pair{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (1, 6) => (3, 4)
 (2, 7) => (5, 7)

dRNA examples:

julia> kmer_to_sig(rna_aln; from=2, to=4) |> collect
3-element Vector{Pair{Tuple{Int64, Int64}}}:
 (1, 6) => (15, 17)
 (2, 7) => (13, 14)
 (3, 9) => missing

julia> kmer_to_sig(rna_aln; from=2, to=5, by=:signal) |> collect
2-element Vector{Pair{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (7, 12) => (5, 7)
 (8, 13) => (3, 4)

See also pos_to_sig for position-to-signal pairs and steps for full alignment (with sequence and signal gaps).

pos_to_sig(
    sig_aln::SignalAlignment;
    from::Integer=1,
    to::Integer=typemax(Int64),
    by::Symbol=:sequence
)

Returns a generator for iterating over sequence position to signal interval pairs. K-mers are converted to a single position using the kmer_center parameter (e.g. k-mer (41, 49) will map to position 45 if kmer_center is 5). Gaps in the sequence will be skipped. The function allows limiting the results to specific sequence position or signal intervals. All intervals used for limiting the results are closed (inclusive) and any overlapping position or signal interval will be included. Results will be sorted in increasing order of sequence positions.

Arguments

• sig_aln::SignalAlignment: signal alignment.
• from::Integer: start index of the slice (inclusive).
• to::Integer: end index of the slice (inclusive).
• by::Symbol: slice by :sequence or :signal.

Examples

These are based on the examples from F5c's documentation:

• DNA: https://hasindu2008.github.io/f5c/docs/output#dna-example
• dRNA: https://hasindu2008.github.io/f5c/docs/output#direct-rna-example

(also shown in the Introduction above) Note: f5c uses 0-based indices and shows the k-mer start. We use 1-based indices and use the full k-mer interval. For this example, let's assume kmer_size=6 and kmer_center=3 (the third base of the 6-mer influences the signal the most).

DNA examples:

julia> pos_to_sig(dna_aln; from=2, to=4) |> collect
2-element Vector{Pair{Int64, Tuple{Int64, Int64}}}:
 3 => (3, 4)
 4 => (5, 7)

julia> pos_to_sig(dna_aln; from=2, to=5, by=:signal) |> collect
2-element Vector{Pair{Int64, Tuple{Int64, Int64}}}:
 3 => (3, 4)
 4 => (5, 7)

dRNA examples:

julia> pos_to_sig(rna_aln; from=2, to=4) |> collect
2-element Vector{Pair{Int64, Tuple{Int64, Int64}}}:
 3 => (15, 17)
 4 => (13, 14)

julia> pos_to_sig(rna_aln; from=2, to=5, by=:signal) |> collect
2-element Vector{Pair{Int64, Tuple{Int64, Int64}}}:
  9 => (5, 7)
 10 => (3, 4)

See also kmer_to_sig for kmer-to-sequence pairs and steps for full alignment (with sequence and signal gaps).

sig_to_kmer(
    sig_aln::SignalAlignment;
    from::Integer=1,
    to::Integer=typemax(Int64),
    by::Symbol=:signal
)

Returns a generator for iterating over signal interval to k-mer pairs. Gaps in the signal will be skipped. The function allows limiting the results to specific signal or sequence position intervals. All intervals used for limiting the results are closed (inclusive) and any overlapping k-mer or signal interval will be included. Results will be sorted in increasing order of signal intervals.

Arguments

• sig_aln::SignalAlignment: signal alignment.
• from::Integer: start index of the slice (inclusive)
• to::Integer: end index of the slice (inclusive).
• by::Symbol: slice by :sequence or :signal.

Examples

These are based on the examples from F5c's documentation:

• DNA: https://hasindu2008.github.io/f5c/docs/output#dna-example
• dRNA: https://hasindu2008.github.io/f5c/docs/output#direct-rna-example

(also shown in the Introduction above) Note: f5c uses 0-based indices and shows the k-mer start. We use 1-based indices and use the full k-mer interval. For this example, let's assume kmer_size=6 and kmer_center=3 (the third base of the 6-mer influences the signal the most).

DNA example

julia> sig_to_kmer(dna_aln; from=2, to=12) |> collect
5-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => (1, 6)
   (5, 7) => (2, 7)
   (8, 9) => (3, 8)
 (10, 11) => missing
 (12, 12) => (4, 9)

julia> sig_to_kmer(dna_aln; from=2, to=5, by=:sequence) |> collect
5-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => (1, 6)
   (5, 7) => (2, 7)
   (8, 9) => (3, 8)
 (10, 11) => missing
 (12, 12) => (4, 9)

dRNA example

julia> sig_to_kmer(rna_aln; from=2, to=12) |> collect
5-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => (8, 13)
   (5, 7) => (7, 12)
   (8, 9) => (6, 11)
 (10, 11) => missing
 (12, 12) => (5, 10)

julia> sig_to_kmer(rna_aln; from=2, to=5, by=:sequence) |> collect
3-element Vector{Pair{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (12, 12) => (5, 10)
 (13, 14) => (2, 7)
 (15, 17) => (1, 6)

See also sig_to_pos for signal-to-pos pairs and steps for full alignment (with sequence and signal gaps).

sig_to_pos(
    sig_aln::SignalAlignment;
    from::Integer=1,
    to::Integer=typemax(Int64),
    by::Symbol=:signal
)

Returns a generator for iterating over signal interval to sequence position pairs. K-mers are converted to a single position using the kmer_center parameter (e.g. k-mer (41, 49) will map to position 45 if kmer_center is 5). Gaps in the signal will be skipped. The function allows limiting the results to specific signal or sequence position intervals. All intervals used for limiting the results are closed (inclusive) and any overlapping position or signal interval will be included. Results will be sorted in increasing order of signal intervals.

Arguments

• sig_aln::SignalAlignment: signal alignment.
• from::Integer: start index of the slice (inclusive)
• to::Integer: end index of the slice (inclusive).
• by::Symbol: slice by :sequence or :signal.

Examples

These are based on the examples from F5c's documentation:

• DNA: https://hasindu2008.github.io/f5c/docs/output#dna-example
• dRNA: https://hasindu2008.github.io/f5c/docs/output#direct-rna-example

(also shown in the Introduction above) Note: f5c uses 0-based indices and shows the k-mer start. We use 1-based indices and use the full k-mer interval. For this example, let's assume kmer_size=6 and kmer_center=3 (the third base of the 6-mer influences the signal the most).

DNA example

julia> sig_to_pos(dna_aln; from=2, to=12) |> collect
5-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => 3
   (5, 7) => 4
   (8, 9) => 5
 (10, 11) => missing
 (12, 12) => 6

julia> sig_to_pos(dna_aln; from=2, to=12, by=:sequence) |> collect
7-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => 3
   (5, 7) => 4
   (8, 9) => 5
 (10, 11) => missing
 (12, 12) => 6
 (13, 14) => 9
 (15, 17) => 10

dRNA example

julia> sig_to_pos(rna_aln; from=2, to=12) |> collect
5-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => 10
   (5, 7) => 9
   (8, 9) => 8
 (10, 11) => missing
 (12, 12) => 7

julia> sig_to_pos(rna_aln; from=2, to=12, by=:sequence) |> collect
7-element Vector{Pair{Tuple{Int64, Int64}}}:
   (3, 4) => 10
   (5, 7) => 9
   (8, 9) => 8
 (10, 11) => missing
 (12, 12) => 7
 (13, 14) => 4
 (15, 17) => 3

See also sig_to_kmer for signal-to-kmer pairs and steps for full alignment (with sequence and signal gaps).

aligned_signal(
    sig_aln::SignalAlignment,
    slow5::Slow5.Slow5File;
    to_pa::Bool=false,
    read_id::Union{String, Missing}=missing,
    read_id_tag::String="pi"
)

Returns the raw signal segments, aligned to the sequence positions, e.g.

[
    (1, [899, 886, 887, 903, 890]),
    (2, [1210, 1199, 1213, 1202])
]

Arguments

• sig_aln::SignalAlignment: SignalAlignment instance.
• slow5::Slow5.Slow5File: A handle to the SLOW5/BLOW5 file with the raw signal.
• to_pa::Bool: If true, will convert the signal to picoampers. Default: false.
• read_id::Union{String, Missing}: The read id in the SLOW5/BLOW5 file. If not provided, it'll try to extract it from the SAM/BAM using the read_id_tag parameter to query the SAM/BAM tags.
• read_id_tag::String: SAM/BAM tag with the id of the raw record in the SLOW5/BLOW5 file. A single raw signal sometimes may contain multiple molecules, which are going to be split during base calling. Dorado uses the pi tag to store the parent read. Default: "pi".

aligned_signal(
    sig_aln::SignalAlignment,
    slow5::Slow5.Read;
    to_pa::Bool=false
)

Returns the raw signal segments, aligned to the sequence positions, e.g.

[
    (1, [899, 886, 887, 903, 890]),
    (2, [1210, 1199, 1213, 1202])
]

Arguments

• sig_aln::SignalAlignment: SignalAlignment instance.
• slow5_record::Slow5.Read: A record from the SLOW5/BLOW5 with the raw data for the read.
• to_pa::Bool: If true, will convert the signal to picoampers. Default: false.

aligned_signal_summary(
    sig_aln::SignalAlignment,
    slow5::Slow5.Slow5File;
    summary_functions=[mean],
    to_pa=false,
    read_id::Union{String, Missing}=missing,
    read_id_tag::String="pi"
)

Executes summary functions on the segments of the signal aligned to each position and returns a DataFrame with the results. The DataFrame will have the following columns:

• pos: sequence position.
• kmer_start: first position of the k-mer, aligned to the segment.
• kmer_end: last position of the k-mer, aligned to the segment.
• dwell: duration that the k-mer spent in the pore, in seconds.
• <summaries>: one column per given summary function. By default it will be only the mean level of the event.

Note: The data frame will be according to the sequencing direction of the signal and the position and kmer columns will have missing values for sequence deletions (regions of the signal that didn't align to any sequence).

Arguments

• sig_aln::SignalAlignment: SignalAlignment instance.
• slow5::Slow5.Slow5File: A handle to the SLOW5/BLOW5 file with the raw signal.
• summary_functions: Functions that will be evaluated on each event in the signal. Can be a vector of functions or :name => fn pairs.
• to_pa::Bool: If true, will convert the signal to picoampers. Default: false.
• read_id::Union{String, Missing}: The read id in the SLOW5/BLOW5 file. If not provided, it'll try to extract it from the SAM/BAM using the read_id_tag parameter to query the SAM/BAM tags.
• read_id_tag::String: SAM/BAM tag with the id of the raw record in the SLOW5/BLOW5 file. A single raw signal sometimes may contain multiple molecules, which are going to be split during base calling. Dorado uses the pi tag to store the parent read. Default: "pi".

Examples

E.g. for dRNA sequencing which goes in 3′->5′ direction:

julia> aligned_signal_summary(
    aln,
    blow5;
    read_id=BAM.tempname(reads[1]),
    to_pa=true,
    summary_functions=[mean, :geometric_mean => xs -> reduce(*, xs)^(1/length(xs))])
125×6 DataFrame
 Row │  pos     kmer_start  kmer_end  dwell    mean      geometric_mean
     │  Int64?  Int64?      Int64?    Float64  Float64   Float64
─────┼─────────────────────────────────────────────────────────────────
   1 │     429         424       432  0.00125   69.2437         69.1985
   2 │     428         423       431  0.01875   91.5471         91.2353
   3 │     427         422       430  0.0045    102.751         102.495
   4 │     426         421       429  0.0055    90.1499         89.9085
   5 │     425         420       428  0.0045    68.6686         68.3237
  ⋮  │    ⋮         ⋮          ⋮         ⋮        ⋮            ⋮
 121 │      24          19        27  0.011     67.7238         67.3204
 122 │ missing     missing   missing  0.01025   58.6041         58.4725
 123 │      22          17        25  0.00175   63.3077         63.3012
 124 │      21          16        24  0.0025    67.4599         67.3375
 125 │      20          15        23  0.005     90.1732         89.9791

get_pa(blow5_record::Slow5.Read)

Returns the read's signal converted to picoampers.

pA = (raw_signal + offset) * range / digitisation))

Arguments

• blow5_record::Slow5.Read: a SLOW5/BLOW5 record.

scale_signal_f5c(pa_signal::Vector{Float64}, read::XAM.XAMRecord)

Get the signal scaled to the model using scaling coefficients from f5c. Expects a signal in picoampers and the BAM record for the read populated with f5c tags. The formula is:

scaled pA current values = (pA - sh) / sc

where pA is the signal in picoampers and sh and sc are f5c scaling parameters written as tags to the SAM/BAM.

Arguments

• pa_signal::Vector{Float64}: signal in picoampers.
• read::XAM.XAMRecord: SAM/BAM record that has f5c tags (sh and sc).

get_raw_read(
    read::XAM.XAMRecord,
    slow5::Slow5.Slow5File;
    read_id::Union{String, Missing}=missing,
    read_id_tag::String="pi"
)

Returns the SLOW5/BLOW5 record associated with the given SAM/BAM read.

If read_id is provided it'll look for record with that id. Otherwise, it'll try to find the SLOW5/BLOW5 record id by looking at the read_id_tag tag in the SAM/BAM record. If the read_id_tag is missing from the record, it'll try to look for a SLOW5/BLOW5 record with id equal to the SAM/BAM QNAME.

Arguments

• read::XAM.XAMRecord: SAM/BAM record for which to find the SLOW5/BLOW5 record.
• slow5::Slow5.Slow5File: A handle to the SLOW5/BLOW5 file with the raw signal.
• read_id::Union{String, Missing}: The read id in the SLOW5/BLOW5 file. If not provided, it'll try to extract it from the SAM/BAM using the read_id_tag parameter to query the SAM/BAM tags.
• read_id_tag::String: SAM/BAM tag with the id of the raw record in the Slow5/Blow5 file. A single raw signal sometimes may contain multiple molecules, which are going to be split during base calling. Dorado uses the pi tag to store the parent read. Default: "pi".

segment_signal(
    signal::Vector{<:Real};
    window::Integer=10,
    pvalue_threshold::AbstractFloat=0.05,
    clip::Integer=0
)

A simple signal segmentation strategy that slides two windows of the given size over the signal and performs a t-test between them at each position to find breaking points. The function will then run peak calling to find the positions of the most significant p-values (discarding the ones above pvalue_threshold).

Note: the default window size of 10 was selected in tests with RNA004. For DNA (R10.4) a velue of 5 was found to work better.

Arguments

• signal::Vector{<:Real}: A vector with the signal measurements.
• window: Size of the window. Two adjacent windows of the given size will be used.
• pvalue_threshold: Less significant positions will be ignored from the peak calling. I.e. positions with a t-test p-value > pvalue_threshold.
• clip: Clip the values in each window to ignore likely "transition values".

Returns

A vector of positions along the signal, representing the boundaries between the segments.

polyA_range(
    aln::SignalAlignment,
    seq::String,
    signal::Vector{<:Real};
    min_len::Integer=200,
    homoA_level_low::Number=70,
    homoA_level_high::Number=100,
    block::Integer=100
)

A quick and simple heuristic algorithm for finding the start and end of the poly(A) region in a signal. The method essentially will start from the signal region aligned to the last non-poly(A) k-mers of the sequence and look back until it finds a region resembling a poly(A) stretch. It'll then expand in both directions until it reaches levels outside the expected bounds.

Note: the function also has a method that does not require the signal-to-sequence alignment and relies solely on the raw signal.

Arguments

• aln::SignalAlignment: the signal alignment.
• seq::String: the sequence of the read.
• signal::Vector{<:Real}: the raw signal for the read in picoampers.
• min_len::Integer: minimum length of the poly(A) region (in sample points).
• homoA_level_low::Number: Minimum level for the poly(A) signal intensity.
• homoA_level_high::Number: Maximum level for the poly(A) signal intensity.
• block::Integer: Number of sample points that will be used to test if the signal deviates from the expected level.

Examples

julia> blow5 = slow5_open("raw_reads.blow5", "r")
julia> slow5_idx_load(blow5)
julia> read = first(open(BAM.Reader, "reads.bam"))
julia> aln = SignalAlignment(read; pore=:RNA004)
julia> raw_read = get_raw_read(read, blow5)
julia> signal = get_pa(raw_read)
julia> polyA_start, polyA_end = polyA_range(aln, string(BAM.sequence(read)), signal)
(3286, 4953)

Returns

A tuple of (from, to) positions.

polyA_range(
    signal::Vector{<:Real};
    min_len::Integer=200,
    homoA_level_low::Number=70,
    homoA_level_high::Number=100,
    block::Integer=10
)

A quick and simple heuristic algorithm for finding the start and end of the poly(A) region in a signal with no additional information. The method will essentially start from the beginning of the signal (assuming 3′→5′ sequencing) and look for an extended region between a given lower and upper picoamper bounds.

Note: the function also has a method that uses the signal-to-sequence alignment, which should provide more reliable results in most cases.

Warning: this function is experimental.

Arguments

• signal::Vector{<:Real}: the raw signal for the read in picoampers.
• min_len::Integer: minimum length of the poly(A) region (in sample points).
• min_start::Integer: minimum starting position of the poly(A) region (in sample points).
• homoA_level_low::Number: Minimum level for the poly(A) signal intensity.
• homoA_level_high::Number: Maximum level for the poly(A) signal intensity.
• block::Integer: Number of sample points that will be used to test if the signal deviates from the expected level.

Examples

julia> blow5 = slow5_open("raw_reads.blow5", "r")
julia> slow5_idx_load(blow5)
julia> read = first(open(BAM.Reader, "reads.bam"))
julia> raw_read = get_raw_read(read, blow5)
julia> signal = get_pa(raw_read)
julia> polyA_start, polyA_end = polyA_range(signal)
(3286, 4953)

Returns

A tuple of (from, to) positions.

polyA_length(
    aln::SignalAlignment,
    read::XAM.XAMRecord,
    blow5_record::Slow5.Read,
    emp_coef::Number=1.7
)

A heuristic method for estimating the length of the poly(A) tail. The method will try to find the poly(A) region in the signal and will calculate its duration. It'll estimate the sequencing speed for the read from the read sequence and signal length and will multiply the poly(A) duration by the sequencing speed weighted by an empirically determined coefficient.

Arguments

• aln::SignalAlignment: signal alignment.
• read::XAM.XAMRecord: SAM/BAM record for the read.
• blow5_record::Slow5.Read: The Slow5/Blow5 record for the read.
• emp_coef: an empirically determined coefficient used to scale the estimation.

The result of searching for a sequence in a raw signal.

Fields

• query::AbstractString: The query sequence.
• query_sig::Vector{Float64}: The expected signal for the query sequence. This is what the function actually searched for in the target signal.
• breaks::Vector{Int64}: The segmentation boundaries found in the target signal.
• start_event::Int64: The position of the first event of the best match, found in the segmented target signal. of the best match.
• start_signal::Int64: The corresponding start position of the best match in the unsegmented target signal.
• cost::Float64: The cost calculated for the best match.
• probability::Float64: The estimated probability (if calculated).

search_sequence(
    seq::AbstractString,
    signal::Vector{<:Real},
    kmer_model::Dict{String, <:Number};
    kmer_size::Union{Integer, Missing}=missing,
    is_reverse::Union{Bool, Missing}=missing,
    pore::Union{Symbol, Missing}=missing,
    segmentation_window::Integer=10,
    segmentation_clip::Integer=0,
    radius::Integer=6,
    estimate_proba::Union{Nothing, Symbol}=nothing,
    num_samples::Integer=40,
    seed::Union{Integer, Missing}=missing
)

Search for a sequence in the raw signal. The function will construct the expected signal for the given sequence using the kmer_model and use dynamic time-warping to look for a match in the raw signal.

If estimate_proba=nothing, the function will just return information about the best match, including the calculated cost, but without a probability estimation. The function can also attempt to estimate the probability that the match is a true occurance of the query sequence. It supports two methods to do so. When estimate_proba=:sample, it will generate num_samples random sequences and use those as a null distribution to calculate the probability. This is time consuming. When estimate_proba=:linear it will use a built-in linear classification model to estimate the probability.

Arguments

• seq::AbstractString: Query sequence.
• signal::Vector{<:Real}: Signal to search in.
• kmer_model::Dict{String, <:Number}: The k-mer model. Should be a dictionary of k-mer => level pairs.
• kmer_size::Union{Integer, Missing}: Size of the k-mer for the pore. Either kmer_size or pore should be set.
• is_reverse::Union{Bool, Missing}: Must be true when the sequencing is in 3′→5′ direction. Either is_reverse or pore should be set.
• pore::Union{Symbol, Missing}: The type of pore used in the sequencing. Supported values: :r9_4, :r10_4, :RNA001, :RNA002, RNA004. One must set either pore or kmer_size and is_reverse.
• segmentation_window::Integer: Size of the window used to segment the signal. Passed on to segment_signal.
• segmentation_clip::Integer: Number of most extreme sample measurements that will be ignored (lowest and highest) when segmenting. Passed on to segment_signal.
• radius::Integer: Radius for the dynamic time-warping.
• estimate_proba::Union{Nothing, Symbol}: Method for estimating the probability for a match. Can be: nothing, :sample, :linear. See above what each method does.
• num_samples::Integer: How many samples to use if estimate_proba=:sample.
• seed::Union{Integer, Missing}: Seed to use when generating the random sequences when estimate_proba=:sample.

seq_to_sig(
    seq::AbstractString,
    kmer_size::Integer,
    kmer_model::Dict{<:AbstractString, <:Number},
    is_reverse::Bool
)

Generate expected signal levels for a sequence, based on a k-mer model. Note: This function won't simulate a realistic signal, as it produces a single data point per k-mer.

Arguments

• seq::AbstractString: Sequence of nucleotides.
• kmer_size::Integer: The k-mer size of the model.
• kmer_model::Dict{<:AbstractString, <:Number}: The k-mer model, as a dictionary of k-mer => current level pairs.
• is_reverse::Bool: Should be true if the sequencing is in 3′→5′ direction (e.g. RNA004).

DNA Example (R10.4)

import DataFrames
import CSV

kmer_model = DataFrame(CSV.File("dna_r10.4.1_e8.2_400bps_9mer_model.txt", header=[:kmer, :level]))
kmer_model = Dict(String.(d.kmer) .=> d.level)

seq_to_sig("ACGTGCTGCTGACATGACAAGCTTCGCACAGGTACGTA", 9, kmer_model, false)

RNA Example (RNA004)

import DataFrames
import CSV

kmer_model = DataFrame(CSV.File("9mer_levels_v1.txt", header=[:kmer, :level]))
kmer_model = Dict(String.(d.kmer) .=> d.level)

seq_to_sig("ACGTGCTGCTGACATGACAAGCTTCGCACAGGTACGTA", 9, kmer_model, true)