Determining a p-value for a test statistic that depends on other test statistics

Question

Sorry for the confusing wording of the title. If some has any better way to word it, please feel free to change it.

Background

For those unfamiliar with bioinformatics data, I have data from a microRNA-seq experiment. In this experiment, small pieces of our DNA, called reads, are extracted, aligned to the genome, and each part (microRNA) of the genome that gets a read aligned to it, gets a count. In the end, we get counts for each microRNA in each sample. The data looks something like this

  Mature.miRNA.ID DRX003170_Lung DRX003171_Lung DRX005946_Saliva
1 hsa-let-7a-2-3p            0.0            0.0              0.0
2   hsa-let-7a-3p            0.9            1.9              0.0
3   hsa-let-7a-3p            0.9            1.9              0.0
4   hsa-let-7a-5p        71231.0        81214.2           1534.5
5   hsa-let-7a-5p        71234.1        81219.4           1448.0

Note: count data is supposed to be integer, but these values have been normalized using a technique called RPM. The normalized count for microRNA i in sample j is $$10^6\frac{c_{ij}}{\sum_{k=1}^{n} c_{kj}}$$ Basically, the normalized count for a microRNA in a sample reflects what fraction of the reads mapped to that microRNA in that sample (multiplied by a million)

Problem

It might happen that even if some microRNA is not present in the cell, a read can still align to it, giving it a count (some microRNAs are similar so misalignments can happen). So if a microRNA has a count of, say 4.2, we have to ask, is this because it was present in the cell, or were the reads actually supposed to go to another microRNA and got mapped to this one accidentally? You can see here how the p-value can get complicated. If microRNA A and microRNA B are very similar and A has a count of 532.1 while B has a count of 0, we should be able to safely conclude that A was present in the cell and B wasn't. But if A has a count of 532.1 and B has a count of 3.4, you can't confidently say that B was present in the cell. But if A had a count of 4.1 and B had a count of 3.4, you might be a bit more confident that B was present in the cell. So we can see how the p-value for presence of B depends on the count of A. This can get more complicated when you have more than two related microRNAs. Say A, B, C, D, E are all similar. Then the p-value for A's presence depends on the counts of B, C, D, and E.

So I'm wondering what kind of test I can use to determine if A is present given the counts of its related microRNAs. I hope the question makes sense.

Answer 1

what kind of test I can use to determine if A is present given the counts of its related microRNAs?

The probability that a miRNA is actually absent from a sample is based not only on the probability of finding false counts when it is absent but also on the probability of finding a count when it is present. For example, your claim

If microRNA A and microRNA B are very similar and A has a count of 532.1 while B has a count of 0, we should be able to safely conclude that A was present in the cell and B wasn't.

isn't necessarily true. It might just be that B is present at a very low level. As you don't sequence and map all of the RNA copies in a biological sample, you might just have missed getting any mappable copies of B in the random subsample of RNA that was submitted for sequencing.

There's also a question about just what you mean by "related" miRNAs. Functionally related miRNAs might often be expressed together, so that the absence/presence of one might be taken to provide evidence for absence/presence of the other.

You, however, seem to have in mind an alternative--sequence-related miRNAs, which would be used together in an opposite way. If the sequences are close enough that a base-read error from A leads to an erroneous mapping to B, the presence of high levels of A might be used as evidence against the presence of B. For that type of "related" miRNAs it might be possible to use base-calling error rates to estimate the probability of mis-mapping, but I'm not familiar with that in miRNA-seq.

I'm more familiar presence/absence calls for gene DNA mutations, estimated for example by MutSigCV software. In that context:

A critical component of MutSigCV is the background model for mutations, the probability that a base is mutated by chance.

In your application, the corresponding issue is whether you find an erroneous miRNA mapping by chance. Mapping errors in miRNA-seq, unfortunately, seem to depend heavily on the particular methods used, so there might be no one-size-fits-all model for miRNA reads that combines read-mapping errors with the probability of true detection to get a presence/absence call.

Others reading this thread should know that this type of presence/absence call is not the usual use case of miRNA-seq data. More typical is evaluating differential expression of miRNAs under a set of experimental conditions. Differential expression analysis does share read-count information among samples and genes to get refined estimates. See for example the DESeq2, edgeR, and limma packages in Bioconductor. The type of negative binomial modeling of original counts used in the first two of those packages might provide some hints for how to approach the presence/absence calling addressed in this question.