# Loess regression on genomic data

I am struggling with R loess function in R. I have a dataframe on which I would locally weighted polynomial regression For each ‘Gene’ is associated a Count (log10 transformed) which gives information regarding the gene expression. For each Gene is associated an ‘Integrity’ measurement (span 0-100) which tells you the quality of the ‘Count’ measurement for each gene. As a general principle, higher is the ‘Integrity’, more reliable is the ‘Count’ for the specific Gene. Below is reported a chunk of the dataframe Sample dataframe:

Gene Integrity Count
ENSG00000198786.2 96.6937 3.55279
ENSG00000210194.1 96.68682 1.39794
ENSG00000212907.2 94.81709 2.396199
ENSG00000198886.2 93.87207 3.61595
ENSG00000198727.2 89.08319 3.238548
ENSG00000198804.2 88.82048 3.78326

I would like to use loess to predict the ‘true’ value of genes with low ‘Integrity’ values (since less reliable).

I) Should I pre-process my dataframe in order to correctly apply loess ? From a pletora of examples I observed sinusoidal distributions of points (A), while my dataset seem distributed in a ‘rollercoaster’-like fashion (B).

II) I have around 15000 observations and I observed that in dependence of the number of observations given, the regression has different a performance. As an example, the fitted values predicted on 100 genes ( C) are clearly more realistic than the fitted values obtained for the very same genes but performing the regression on the whole dataset ( D).

III) How should I run loess? I cannot understand how to run loess with the correct syntax to differentially weighted observations: -1 loess( Count ~ Integrity, weight=None) -2 loess( Count ~ 1:nrow(dataframe), weight=Integrity)

I performed several tests. Fig. E-F used loess (stats), Fig. G-H run weightedloess (limma). I used 2 different packages since, from the loess docs it is clear that prior weights are set based on x distance between points. weightedloess function allow the user to give priors in order to perform regression.

E) loess(Count ~ Integrity),degree=2,span=0.1)
F) loess(Count ~ 1:nrow(df)),weigths=’Integrity’,degree=2,span=0.1)
G) weightedLowess(x=1:nrow(df), y=Count, weigths=’Integrity’, span=0.1)
H) weightedLowess(x=1:nrow(df), y=order(Count), weigths=’Integrity’, span=0.1)


Please find enclosed images cited in the question.

I would like to use loess to predict the ‘true’ value of genes with low ‘Integrity’ values.

That's probably not a good idea. Look at your example F, which plots some log-gene-expression values against Integrity. If you use the smoothed plot to get a value of log-gene-expression corresponding to a particular Integrity, all of your values would be in the range of 2.5 to 3.0. That's throwing away a lot of detail. Even if individual values aren't equally reliable, the individual expression values (or their averages over multiple samples) are presumably the most reliable values that you have. There's no good reason to replace the observed values with what otherwise is close to the average expression across all genes.

You might, however, want to put less emphasis in your model on the values with the lowest Integrity. That's something that can be done by the limma package, using a similar approach in the opposite direction: modeling a type of Integrity as a function of gene-expression value.

When your model involves multiple genes measured in multiple samples, in an overall model you might weight each gene's value inversely with respect to the variance of its gene-expression estimate. With the sample sizes usually found in this type of study, however, those estimates can vary widely among genes even if the underlying distributions of the gene-expression values are similar. Relying on variance estimates from individual genes can be misleading.

Thus limma can pool information across genes and estimate the variance as a function of gene-expression level. It then uses the pooled estimate in its modeling. That's the opposite direction from what you propose. I haven't tried to use things like your Integrity values in that way to differentially weight genes in a limma model, but I suspect that it has been done and is documented in Bioconductor help threads.

In terms of your invocations of loess(), it's mostly a matter of how much wiggliness that your want to allow in the curve and whether you want all data points to be treated equally. For wiggliness there's no single best answer; the defaults often work well.

It's not clear to me that differential weighting according to Integrity will do much for you. In your plots with Integrity on the x-axis, that down-weights all points at any Integrity equally, increasingly so as you move from Integrity=100 down toward Integrity=0. But the local regression done by loess() is based on regions of similar values along the x-axis already, so the relative weighting by Integrity won't matter much locally. If you have a very wide span, it will tend to put more emphasis on the information from the right side of the plot.

If you instead model Integrity as a function of gene expression, then you might weight by how precisely the Integrity values themselves are estimated.