NA in likelihood ratio test

[pimentel]

look into why we are getting NA in sleuth_lrt() and look into what p.adjust is doing in this situation (treating them as ones?)

[warrenmcg]

Hey @pimentel,

This issue was also raised in the forums here: link

I took a look at some of my datasets and found some NA values there as well. I did some digging and discovered that the NA values are from the shrink_df step of sleuth_fit. After creating the LOESS fit, the "shrink" value is then predicted from the fit. The target ids that generated NA values have mean_obs values outside of the range used by the LOESS fit. As stated in the ?predict.loess documentation:

When the fit was made using surface = "interpolate" (the default), predict.loess will not extrapolate – so points outside an axis-aligned hypercube enclosing the original data will have missing (NA) predictions and standard errors.

In my case, these NA values were only generated by the 'reduced' model reduced <- formula(~1).

There are a few different things you could do:

1. Just throw a warning that some values fell outside the range of the LOESS fit, and leave them as NAs [simplest at this point]
2. Revert to unshrunken sigma values [bad idea]
3. Attempt to extrapolate to these values [bad idea + I don't know how you would do it]
4. Rerun the LOESS using control = loess.control(surface = "direct") as an option. This would lead to an exact calculation, which could then be extrapolated. I have not tested how long this takes and whether the rest of the values significantly change when switching, but this is potentially worth looking into.

p.adjust ignores NA p-values (as if they didn't happen), so the "N" used to adjust them is just based on non-NA values. Quoting from a very old (2005) thread on making changes to p.adjust:

[T]he non-NA adjusted p-values should be the same as if there were no
NAs in 'p'. This is because NAs in 'p' virtually always indicate, not an
existing but unknown p-value, but a p-value which could not be computed
because there was insufficient data to do so. This means that there is
zero probability of rejecting the null hypothesis for these
cases. Therefore it is unnecessary to adjust the other p-values for these
cases.

[winni2k]

I think I might be running into a similar issue:

When I run sleuth_lrt on 1084 transcripts measured in eight samples across three groups (3, 3, and 2 samples per group), I get 855 NA values in the mean_obs column. I have checked the kallisto estimated counts and tpms for missing values using:

• sleuth_to_matrix(so, which_df='obs_raw', which_units='est_counts')
• sleuth_to_matrix(so, which_df='obs_norm', which_units='est_counts')
• sleuth_to_matrix(so, which_df='obs_norm', which_units='tpm')

, but I have found no missing values in those matrices. My transcripts are highly similar, so I would expect some unevenness in estimated transcript counts between samples in the same group, and I do observe this quite frequently.

My principal concern is that all of my 855 NA values are not counted in the Benjamini-Hochberg procedure, and that the quote above does provide a caveat on the reason for p-values being NA. The proportion of NA values is also quite high, much higher than anything I have seen anyone report here. On this basis is there a chance that that my FDR corrected p-values might be liberal?

@pimentel, have you had a chance to try the loess.control(surface = "direct") fix? Would you recommend I give the fix a shot?

PS Here is a tail of the results table that I get out of sleuth_lrt

     target_id test_stat pval qval      rss sigma_sq   tech_var mean_obs
1079      1077        NA   NA   NA       NA       NA         NA       NA
1080      1078        NA   NA   NA       NA       NA         NA       NA
1081      1079        NA   NA   NA       NA       NA         NA       NA
1082      1080        NA   NA   NA 72.61146 9.064751 0.01168158 4.579041
1083      1082        NA   NA   NA       NA       NA         NA       NA
1084      1083        NA   NA   NA       NA       NA         NA       NA
      var_obs sigma_sq_pmax smooth_sigma_sq final_sigma_sq degrees_free
1079       NA            NA              NA             NA           NA
1080       NA            NA              NA             NA           NA
1081       NA            NA              NA             NA           NA
1082 9.076432      9.064751       0.8138833       9.064751            2
1083       NA            NA              NA             NA           NA
1084       NA            NA              NA             NA           NA

[warrenmcg]

@winni2k What are you analyzing that you only have 1084 transcripts to analyze? That to me seems to be the reason why there are so many NAs coming out of your analysis. If this is the whole transcriptome under analysis, then I would definitely run the version of sleuth that has the second layer of using the direct/exact LOESS, and see what happens.

[winni2k]

Thanks @warrenmcg! I'll give that a try then. Is this functionality still only in the dev branch?

To answer your question: This is not the whole transcriptome. I am looking at a subset of transcripts from a transcriptome. Is this a mis-appropriation of kallisto and/or sleuth?

[warrenmcg]

Hi @winni2k, you should definitely do most of the pipeline on the whole transcriptome.

If this subsetting was done with kallisto, you are definitely distorting the results if there are fragments that could have mapped ambiguously to transcripts not included in the annotation.

If the subsetting was just done with sleuth, then it's the shrinkage estimation step that we're discussing now that will be distorted. Shrinkage estimation, which afaik was first used for RNA-Seq in DESeq and edgeR, is a way of borrowing from all of the transcripts present in the data to better estimate the variance for each individual transcript. You will get different estimates if you look only at a subset, and the way that sleuth does it in combination with so few transcripts is likely the reason why your estimation is turning up so many NAs.

If you're only looking at a subset of transcripts, it's better to focus the subsetting on the testing stage. In most cases, it's still not justified to only test a subset; instead you should just look at how your subset of transcripts perform when tested with the rest of the transcriptome.

However, if it is a justified situation to test just the subset, you can subset out your transcripts from the results table and repeat the p-value adjustment with just that subset. The particulars of your situation will depend on the details of your experiment, so if you want to discuss the approach and you're not comfortable talking about this on a public forum, feel free to email me to follow-up.

[winni2k]

Thank you for this extensive response @warrenmcg. I really appreciate it. I will need to do a bit more reading to figure out how to proceed.