Clarification on GMM LOR calculation and cluster interpretation in Nanocompore

[nawa2179]

Hi,

Thank you for developing and maintaining Nanocompore.

I am currently using the tool and would appreciate some clarification regarding the GMM LOR calculation and the consistency of results across replicates.

I analyzed transcriptomic data from two biological replicates (WT vs KO for a given enzyme), generated using the ONT RNA004 kit, following your recommended pipeline. I performed three analyses:

both replicates combined
replicate 1 only
replicate 2 only

Using the following thresholds:

GMM_chi2_qvalue ≤ 0.1
|GMM_LOR| ≥ 0.5

I obtained highly inconsistent results:

2 sites (replicate 1)
626 sites (replicate 2)
15 sites (combined analysis)

To rule out potential data quality issues, I re-basecalled the data using Dorado and quantified m6A and pseU with Modkit. No major differences were observed between conditions.

I have two main questions:

1. GMM LOR definition
   From the documentation, it seems that GMM LOR is defined as:

ln((WT_C1 / WT_C2) / (KO_C1 / KO_C2))

However, could you clarify how this is actually computed internally?

2. Cluster interpretation
   Which GMM cluster corresponds to modified vs unmodified kmers?

From empirical testing, it appears that the formula may instead behave like:

ln((sample2_mod + 1 / sample2_unmod + 1) / (sample1_mod + 1 / sample1_unmod + 1))

Could you confirm whether:
clusters are assigned based purely on signal distribution or linked to modification status?
a pseudocount is systematically applied?
the LOR is computed from cluster proportions?
Any clarification or pointers to the relevant code sections would be greatly appreciated.

Thank you in advance for your help.

Best,

Nawad

[mzdravkov]

Hi, Nawad!

The difference between the two replicates indeed is rather big. Some things to consider:

1. You have written that your q-value threshold is 0.1 (not sure if it's a typo). That's quite permissive. We typically use something closer to 0.05 or 0.01. Of course the decision would depend on your specific use case and circumstances.
2. Are the sizes of the two libraries similar? Nanocompore would filter out positions with insufficient coverage. The default threshold is minimum 30 reads (for each condition). It may be a good idea to check if: (a) the number of tested positions (i.e. lines in the output file) is similar and (b) the coverage for the significant positions in the two replicates is similar. Differences in the coverage can mean that you're testing wildly different number of positions or that higher coverage of specific transcripts gives you more statistical power to detect things there.

Regarding your questions about how Nanocompore works:

When the GMM test is used, Nanocompore would first run GMM clustering on all data points for the position twice. Once to fit a single cluster to the data and another time to fit two clusters. Then it uses the Bayesian Information Criterion to decide if one or two clusters fit better. If it's one, then the position would not be significant. If it's two, it'll try to use a chi-square test between the condition labels and the cluster assignments (this gives the p-value). However, the clustering is completely blind to the condition labels and uses the signal data alone.

The LOR is calculated, as you guessed, using the formula ln((COND1_C1/COND1_C2) / (COND2_C1/COND2_C2)), although note that COND1 and COND2 is a bit arbitrary (based on the order in which you defined them in the configuration file). So the LOR probably should be used as an absolute value to get the effect size. You're also correct that we use pseudo counts (adding 1 to each of the four numbers in the LOR calculation) to avoid division by zero issues.

Nanocompore will infer which cluster corresponds to the modified state using the depleted_condition parameter in the configuration. It'll look at which cluster is more enriched for the depleted condition and will assume that the other one represents the modified state. Note, this doesn't affect the p-value and LOR outputs, but is important for the count columns (<SAMPLE>_mod and <SAMPLE>_unmod) and for the modification probabilities when using the single molecule resolution functionality.

If you're interested in taking a look at the code, the GMM test is implemented here: https://github.com/tleonardi/nanocompore/blob/master/nanocompore/comparisons.py#L357. In that function you'll also easily spot the ones for getting the cluster counts and calculating the LOR.

Best regards,
Mihail

[nawa2179]

Hi Mihail,

Thank you for the clarification.

As you suggested, I reviewed the library sizes and the number of positions tested. Both samples have a library size of approximately 9 million reads, and the number of positions tested in each replicate is relatively similar:

Replicate 1: 13,225,108 sites

Replicate 2: 12,274,990 sites

Total (replicates 1 + 2): 17,357,889 sites

It appears that the library size is not the source of the problem.

Thank you again for your advice.

Best regards

Nawad

[lmulroney]

Hi Naward,

If I understand your experimental design, you have 2 replicates of WT, and 2 replicates of the enzyme KO.
And you used Nanocompore in these three ways:

1. Both WT replicates compared against both KO replicates (found 15 sites)
2. WT replicate 1 compared against KO replicate 1 (found 2 sites)
3. WT replicate 2 compared against KO replicate 2 (found 626 sites)

Is this correct?

Also, what exactly do you mean when you say that there were no major differences between the replicates when using Dorado and modkit? Because the LoR will be a measure of effect size. And a detected site using modkit may find 1 or 2 reads with a modification out of 100+ that Nanocompore would remove because the proportion of modified reads is too low to be confident. And if you have incomplete KO of the enzyme, there may still be very low activity that can be detected with sensitive techniques like the dorado basecaller.

If you are willing to describe your experimental and analysis design in a little more detail, including assumptions and filtering criteria for both Nanocompore and Dorado, we may be able to offer more concrete advice.

Cheers,
Logan

[nawa2179]

Hi Logan,

Thank you for your feedback.

To provide some additional context, we performed direct RNA sequencing on two conditions (WT and enzymatic KO), each with two replicates, and we are using Nanocompore to identify target genes based on differential modifications.

Regarding our quality control, we used Modkit to quantify global m6A and pseudouridine sites. We observed no major differences in these modification levels between samples, confirming that our sequencing data are technically reliable and that the samples are broadly comparable.

Regarding KO efficiency: we verified inhibition by monitoring the expression of the targeted “writer” enzyme. In both KO replicates, we found very few mapped reads. Furthermore, IGV visualization confirmed that the majority of these remaining mapped reads contain micro-insertions/deletions generated by CRISPR-Cas9, which cause a frameshift and likely result in a complete loss of function.

I would also like to clarify the figures I provided earlier. The numbers I mentioned (13 million, 12 million, and 17 million sites) refer to the total number of positions analyzed by Nanocompore (raw results), not the number of significant sites identified after filtering. For the significant results, I applied the following filters:

GMM chi2 qvalue ≤ 0.1

|GMM_LOR| ≥ 0.5

Using these criteria, I obtained:

Replicate 1 (WT1 vs. KO1): 2 significant sites.

Replicate 2 (WT2 vs. KO2): 626 significant sites.

Combined (WT1+WT2 vs. KO1+KO2): 15 significant sites.

The fact that the total number of tested positions is consistent across replicates confirms that library size and coverage are not the issue. Given these results and our validation of the KO efficiency, I am unable to identify the source of this discrepancy.

Best,

Nawad

[mzdravkov]

Hi Nawad,

Thank you for the information!
I guess that you're knocking out a writer enzyme for a modification that is not supported by the Dorado models (or at least not m6A and pseudouridine) and this is why you're validating with Dorado that there's no change in modification levels between the two conditions. Is that correct?

I also assume that the configuration files for the runs are the same (excluding the sample definitions of course), but still you can double check that there aren't some differences by mistake.

Some things that come to my mind that you can try. In no particular order:

1. Cross-comparing your samples, i.e. compare WT1-KO2, WT2-KO1 and even WT1-WT2, KO1-KO2 could indicate if there's one sample that behaves unexpectedly.
2. You can use the plotting functions to inspect some of the 626 significant sites in the WT2-KO2 run one by one. Do the distributions from the two conditions look different or are they likely false positives.
3. If you haven't changed the default comparison_methods parameter in the configuration file, it would have executed both GMM and Kolmogorov-Smirnov tests. If you look at the KS_intensity_qvalue to select significant sites, is there still such a big difference in the results from the two replicates?
4. I'd take a look at the coverage for the significant positions in WT2-KO2 that are not in WT1-KO1 and compare them to the coverage for these positions in WT1-KO1. Even if the total sizes of the libraries are similar, there may be some redistribution between the transcripts. If that's the case, then the different statistical power could explain the different results.
5. Check the logs for errors. Specifically check if any errors are reported in WT1-KO1 for the transcripts on which the 626 significant sites are.

Best regards,
Mihail