Hello,

I am trying to understand the limma diffSplice function. When calculating the gene expression diffSplice weights the exon expression values by 1/stddev.unscaled^2.

#   Genewise betas
    u2 <- 1/exon.stdev.unscaled^2
    u2.rowsum <- rowsum(u2,geneid,reorder=FALSE)
    gene.betabar <- rowsum(exon.coefficients*u2,geneid,reorder=FALSE) / u2.rowsum

Why does it do that? And why not use the scaled standart deviation 1/(sigma^2+stdev.unscaled^2)? I guess I have also trouble understanding to what the stdev.unscaled vaules actually correspond to.

Further diffSplice scales the moderated t-value by some leverage.

exon.1mleverage <- 1 - (u2 / u2.rowsum[g,,drop=FALSE])
exon.coefficients <- exon.coefficients / exon.1mleverage
exon.t <- exon.t / sqrt(exon.1mleverage)

What is the reason behind that? I could not find anything in the documentation or in the limma userguide. Any pointers to literature or other explenations would be greatly appreciated.

Thanks,

Stefan

The definitions of sigma and stddev.unscaled are given in several places in the limma documentation, see help("MarrayLM-class) or Section 13.2 of the User's Guide.

The use of leverages is a advanced computational trick to convert a coefficient computed across all exons to a coefficient for all exons except one, for the purpose of comparing one exon to the average of all the rest (which is documentated to be the test conducted by the function).

The documentation is intended to explain the input and output from the user-level function diffSplice rather than to explain individual lines of code. The meaning of the t-statistics returned by the function is described in the Details section of the help page for that function.