questions regarding raw signal analysis

[jorisbalc]

Hi Marcus,

I've got a couple of questions regarding the api usage and raw signal analysis. Sorry in advance if some of the questions are super basic as I'm not too familiar with the remora api yet.

1. The region_basecalls = io_read.extract_ref_reg(ref_reg = ref_region, signal_type="norm") returns an object, sufficient for mapping raw signal data as it contains the seq_to_sig array. I'm wondering how these values are transformed as I've analzyed the pA values from sample to sample and my modified sample has roughly ~2pA higher values on average, comapred to non-modified sample (in non-modified sites, such as explicitly selected large kmers of 50-60 bases in length without any CG sites). Is it possible for such difference is caused by the samples being run at different times or even cells or would this strictly be caused by the DNA sample being modified with bulky mods at every CG site? I also wonder what's the best properly transform the norm signal values to represent something like the printout of remora analyze command.
2. I've encountered some reads, where the seq_to_sig array doesn't make any sense to me. The only logical conclusion I've come to is that somehow the strand "skips" and goes through way faster than intended? I'm posting a single example that gets returned by <io_read.extract_ref_reg>. In this particular case, I'm looking at 31-mers centered at CG sites. Would you perhaps help in understanding how such a read happens? Also, what exactly does sig_start denote?

read = io_read.extract_ref_reg()

ReadRefReg(
    read_id='784882a1-b632-4512-8322-6a7227f6183e',
    norm_signal=array([
        55827.85090193, 55757.59090376, 55944.95089888, 55804.43090254,
        55921.53089949, 55687.33090559, 55593.65090803, 55289.19091596,
        55289.19091596, 55476.55091108, 55804.43090254, 56272.83089034,
        55476.55091108, 54446.07093792, 55429.7109123 , 55218.93091779,
        55687.33090559, 55640.49090681, 55195.5109184 , 55008.15092328,
        55781.01090315, 55382.87091352, 55851.27090132, 56249.41089095,
        56296.25088973, 56553.87088302, 59504.79080615, 60160.55078907,
        60137.13078968, 59715.57080066, 59434.53080798, 59668.73080188,
        59738.99080005, 59715.57080066, 59879.51079639, 60020.03079273,
        59668.73080188, 58614.83082933, 57279.89086411, 56577.29088241,
        56577.29088241, 56670.97087997, 56319.67088912, 59434.53080798,
        59200.33081408, 59668.73080188, 59738.99080005, 59762.41079944,
        59481.37080676, 59668.73080188, 59106.65081652, 59528.21080554,
        59387.6908092 , 60160.55078907, 59457.95080737, 59973.19079395,
        59598.47080371, 59598.47080371, 59317.43081103, 59340.85081042,
        59247.17081286, 59926.35079517, 59856.090797  , 60020.03079273,
        59645.31080249, 59598.47080371, 59130.07081591, 59645.31080249,
        59996.61079334, 59949.77079456, 60301.07078541, 59926.35079517,
        59856.090797  , 60066.87079151, 59598.47080371, 60020.03079273,
        59996.61079334, 59528.21080554, 59575.05080432, 59692.15080127,
        59387.6908092 , 59645.31080249, 59528.21080554, 59902.93079578,
        59364.27080981, 59270.59081225, 57537.5108574 , 55429.7109123 ,
        55336.03091474, 55593.65090803, 56108.89089461, 56670.97087997,
        54422.65093853, 53579.53096049, 54399.23093914, 55172.09091901,
        55289.19091596, 55031.57092267, 55195.5109184 , 55218.93091779,
        55499.97091047, 54797.37092877, 55078.41092145, 55265.77091657,
        54961.3109245 , 54820.79092816, 55054.99092206, 55617.07090742,
        55593.65090803, 55734.17090437, 55921.53089949, 55663.9109062 ,
        55734.17090437, 55804.43090254, 55546.81090925, 55781.01090315,
        55382.87091352, 55406.29091291, 55359.45091413, 55944.95089888,
        55382.87091352, 55663.9109062 , 58029.33084459, 63579.87070001,
        63533.03070123
    ]),
    seq='AGTTAATGAGTGTGTCGAAAAACCTGAAATC',
    seq_to_sig_map=array([
        0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6,
        6, 7, 7, 8, 8, 9, 9, 10, 11, 23, 29, 41, 53,
        83, 89, 95, 113, 119, 125
    ]),
    ref_reg=RefRegion(ctg='JB0001.0', strand='+', start=17711, end=17742),
    sig_start=22687
)

3. We noticed some modified reads have 20-50% more <.norm_signal> entries. I assume this correlates to dwell time? In any case, is dwell time information stored in the <io_read> object?
4. Is the <ReadRefReg.seq> the basecalled sequence or the ref? I assume it's the basecalled seq.

As always, I appreciate any help on this

[marcus1487]

Thanks for the thoughtful questions. Not basic at all! You’re digging into important signal-level details, and it’s great to see this level of care in your analysis. I’ll go through each of your points in order:

1. pA Value Variation Between Runs and Samples

Yes, a ~2 pA shift between samples is entirely plausible due to run-to-run variability, including differences in flow cells, pores, temperature, and runtime conditions. This effect can be compounded when comparing different sequencing runs.

To mitigate this and allow signal-level comparisons across samples, we strongly recommend k-mer expected signal normalization. This approach uses canonical signal expectations for each k-mer to normalize out per-read and per-run shifts. You can find this implemented in the Remora notebooks, particularly in:

remora/notebooks/signal_mapping_refinement.ipynb

If you’re comparing groups of reads (e.g. mod vs control), rather than single reads, normalization becomes especially important. Also have a look at the metrics_api.ipynb notebook — it provides guidance for extracting signal metrics at the read or group level:

remora/notebooks/metrics_api.ipynb

2. Unexpected seq_to_sig_map Behavior / Signal Skipping

Without refinement using a k-mer expected signal model, the seq_to_sig_map array is directly inherited from the basecaller, which optimizes for decoding accuracy rather than strict per-base alignment to signal. This means:

• Some bases may be emitted “faster” in signal space if the basecaller had higher information at a portion of the read.
• Conversely, some regions may appear to stretch if information content was low.

This is why you may observe reads where signal indices assigned to bases seem uneven or implausible. To address this:

• Use signal mapping refinement, which aligns observed signal to expected k-mer levels to get a more faithful signal-to-sequence alignment.
• This process helps resolve issues like signal compression or expansion, and is crucial for dwell time or kinetic analysis.

If you can share a specific example read, I’d be happy to help interpret what might be going on.

3. Dwell Time and Signal Length Differences in Modified Reads

You’re exactly right. More norm_signal entries typically correspond to longer dwell times, i.e., more samples collected while a base was in the pore.

You can compute dwell time from the seq_to_sig_map using np.diff(seq_to_sig_map), which gives you the number of signal points associated with each base.

Again though, without refinement, these dwell estimates may be unreliable, especially near modifications where the basecaller behavior can differ.

4. Understanding ReadRefReg.seq vs Basecalls

The ReadRefReg object represents a segment of the read aligned to the reference, so .seq refers to the reference sequence, not the basecalled sequence.

If you want the basecalled view of the signal, you can:

• Use the full io_read object and access .basecall, .norm_signal, and .seq_to_sig_map directly.
• Or extract a ReadBasecallRegion via io_read.extract_basecall_region(...).

The region-based interfaces (like extract_ref_reg or extract_basecall_region) are just conveniences for plotting and extracting manageable segments, but they all reference the same underlying io.Read data.

Let me know if you’d like help walking through the refinement steps with a concrete read. Happy to help troubleshoot or clarify further.

[jorisbalc]

Thanks for the incredibly fast response Marcus! As I understand further analysis requires signal refinement. I will take a look at that next. As for the read, if the spoiler with the read info isn't working, I'm posting it with this comment. I except the problem with the read is exactly what you described.

ReadRefReg(
    read_id='784882a1-b632-4512-8322-6a7227f6183e',
    norm_signal=array([
        55827.85090193, 55757.59090376, 55944.95089888, 55804.43090254,
        55921.53089949, 55687.33090559, 55593.65090803, 55289.19091596,
        55289.19091596, 55476.55091108, 55804.43090254, 56272.83089034,
        55476.55091108, 54446.07093792, 55429.7109123 , 55218.93091779,
        55687.33090559, 55640.49090681, 55195.5109184 , 55008.15092328,
        55781.01090315, 55382.87091352, 55851.27090132, 56249.41089095,
        56296.25088973, 56553.87088302, 59504.79080615, 60160.55078907,
        60137.13078968, 59715.57080066, 59434.53080798, 59668.73080188,
        59738.99080005, 59715.57080066, 59879.51079639, 60020.03079273,
        59668.73080188, 58614.83082933, 57279.89086411, 56577.29088241,
        56577.29088241, 56670.97087997, 56319.67088912, 59434.53080798,
        59200.33081408, 59668.73080188, 59738.99080005, 59762.41079944,
        59481.37080676, 59668.73080188, 59106.65081652, 59528.21080554,
        59387.6908092 , 60160.55078907, 59457.95080737, 59973.19079395,
        59598.47080371, 59598.47080371, 59317.43081103, 59340.85081042,
        59247.17081286, 59926.35079517, 59856.090797  , 60020.03079273,
        59645.31080249, 59598.47080371, 59130.07081591, 59645.31080249,
        59996.61079334, 59949.77079456, 60301.07078541, 59926.35079517,
        59856.090797  , 60066.87079151, 59598.47080371, 60020.03079273,
        59996.61079334, 59528.21080554, 59575.05080432, 59692.15080127,
        59387.6908092 , 59645.31080249, 59528.21080554, 59902.93079578,
        59364.27080981, 59270.59081225, 57537.5108574 , 55429.7109123 ,
        55336.03091474, 55593.65090803, 56108.89089461, 56670.97087997,
        54422.65093853, 53579.53096049, 54399.23093914, 55172.09091901,
        55289.19091596, 55031.57092267, 55195.5109184 , 55218.93091779,
        55499.97091047, 54797.37092877, 55078.41092145, 55265.77091657,
        54961.3109245 , 54820.79092816, 55054.99092206, 55617.07090742,
        55593.65090803, 55734.17090437, 55921.53089949, 55663.9109062 ,
        55734.17090437, 55804.43090254, 55546.81090925, 55781.01090315,
        55382.87091352, 55406.29091291, 55359.45091413, 55944.95089888,
        55382.87091352, 55663.9109062 , 58029.33084459, 63579.87070001,
        63533.03070123
    ]),
    seq='AGTTAATGAGTGTGTCGAAAAACCTGAAATC',
    seq_to_sig_map=array([
        0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6,
        6, 7, 7, 8, 8, 9, 9, 10, 11, 23, 29, 41, 53,
        83, 89, 95, 113, 119, 125
    ]),
    ref_reg=RefRegion(ctg='JB0001.0', strand='+', start=17711, end=17742),
    sig_start=22687
)

[jorisbalc]

Hi marcus,

Thanks for the tips! There are a few more things I'm curious about. As far as I understand do_fix_gauge normalization is very strict and limits the metrics that I can pull out of the dataset. As an example, I'd want to see average |Δ| changes across multiple x-mers containing the modified CG sites, for the purpose of comparing how different mods influence the magnitude of the signal shift from control groups. Would it makes sense to not normalized the data here? Or is such a metric generally not really reliable due to the difference in in run parameters?

To generalize: I'm trying to look at multiple CG centered k-mers and calculate the mod vs control |delta| shift (and dwell time)

Also, when calculating base specific signal metrics, what sort of signal trim is recommended, if any? Is the default=1 on both start/end points?

[marcus1487]

Thanks again for your thoughtful questions. You’re asking exactly the right things for doing robust signal-level analysis. I’ll walk through your points one by one, clarifying what each normalization and parameter does, and offering some guidance for your analysis goals.

do_fix_gauge Normalization

The do_fix_gauge flag only applies to the input k-mer model table, and not to the processing of your observed signal data. It ensures that signal level expectations are on a fixed, consistent scale across contexts. This is almost always recommended for accurate mapping and comparison.

If your concern is about normalization affecting downstream signal metrics, the more relevant flag is likely do_rough_rescale, which controls whether observed signals are rescaled relative to expected values. More on that below.

Should I Normalize for Δ Signal Comparison?

Yes, normalization is highly recommended if you want to compare average signal shifts (|Δ|) between modified and control reads, especially across different runs, flow cells, or sequencing sessions. Even seemingly small differences in run parameters (temperature, voltage drift, pore batches) can lead to global shifts in current levels.

The do_rough_rescale=True option aligns each read’s signal to expected canonical k-mer levels, which helps cancel out these run-specific effects. This is robust and flexible enough to support meaningful mod-vs-control comparisons. You can always compute |Δ| directly from the rescaled values.

However, if you have a tightly controlled experimental setup and want to examine raw signal shifts without normalization, that can be done, just be aware that apparent shifts may reflect variance due to unwanted effects (run settings, applied voltage, flowcell variations, etc) rather than biological effects.

If you can share more details about the run parameter differences, I’m happy to help evaluate whether normalization is a must in your case.

Signal Trimming

Trimming isn’t always necessary, but can improve robustness for some analyses:

• Trimming is not performed by default in the production of an io.Read object.
• For highly sensitive comparisons, trimming 1–2 bases from each end of the region (especially after signal mapping refinement) helps reduce edge alignment artifacts.

So, I’d recommend evaluating this based on your observed consistency: if you’re seeing noisy or inconsistent terminal dwell/signal behavior, trimming is worth it.

Summary Recommendations for Your Use Case

To compute reliable |Δ| and dwell shifts across CG-centered kmers:
1. Use signal mapping refinement to align signal to expected k-mers.
2. Enable do_fix_gauge when using a k-mer table.
3. Apply do_rough_rescale=True to minimize inter-run variance.
4. Trim 1–2 bp from region ends if you’re seeing noisy edges after mapping.

Happy to walk through an example read or refine your analysis further.