In areyesq89/HumanTissuesDEU: Analysis of tissue-dependent exon usage in humans

knitr::opts_chunk$set(tidy = FALSE,
                      cache = FALSE,
                      dev = "png",
                      message = FALSE,
                      error = FALSE,
                      warning = TRUE)

Estimating REUCs

Below are the steps that were followed to estimate relative exon usage coefficients (REUCs). We load the data objects and extract the indexes of the exonic regions of the first gene.

library( HumanDEU )
library( DEXSeq )
data( "dxdObjects" )
dxd <- dxd1
gene1Index <-
    which( mcols( dxd )$groupID %in% "ENSG00000000003" )

We fit the model for each exonic region using a large prior for shrinkage and vector of fixed dispersion estimates. We also use a single core for the calculations, since we are estimating these for a single gene.

dispThis <- rep( 0.1, length.out = length( gene1Index ) )
dispOthers <- rep( 0.1, length.out = length( gene1Index ) )
names( dispThis ) <- rownames( dxd )[gene1Index]
names( dispOthers ) <- rownames( dxd )[gene1Index]
testableVector <- rep( TRUE, length( gene1Index ) )
names( testableVector ) <- rownames( dxd )[gene1Index]
bjp <- SerialParam()

allCoefsWeakShrinkage <-
    fitAllExonsParallel(
        dxd[gene1Index,],
        dispThis,
        dispOthers,
        priorsd=3,
        bjp )

We used these coefficients with weak shrinkage to get dispersion estimates for each exonic region.

rawDisps <-
    estimateDispersionsParallel(
        dxd[gene1Index,],
        allCoefsWeakShrinkage,
        bjp)

head( rawDisps )

Using these dispersion estimates, we run now the fit with a smaller prior.

dispThis <- pmax( pmin( rawDisps[,"dispThis"], 10 ), 1e-6 )
dispOthers <- pmax( pmin( rawDisps[,"dispOthers"], 10 ), 1e-6 )

allCoefsGoodShrinkage <- fitAllExonsParallel(
    dxd[gene1Index,],
    dispThis,
    dispOthers,
    priorsd=0.05, bjp )

allCoefsSub <-
    HumanDEU:::arrangeInto3DArray(
        dxd[gene1Index,],
        allCoefsGoodShrinkage )

We verify that the estimated REUCs are the same as in the REUCs in the pre-computed object.

data("crossCoefs1")
all( allCoefsSub == crossCoefs1[gene1Index,,] )

Estimating p-values

As for the REUCs and RSICs, the code below demonstrates the functions that were used to test for differential exon usage. Due to long running times, we demonstrate the p-value calculation for the exonic regions of the first gene.

dxd <- dxd1

formulaFull <- ~ sample + exon + sex:exon + individual:exon + tissue:exon
formulaNull <- ~ sample + exon + sex:exon + individual:exon

disps <- estimateTestDispersionsParallel(
    dxd[gene1Index,],
    formula=formulaFull,
    bjp )

disps <- data.frame(
    dispThis=disps,
    dispOthers=disps )
rownames(disps) <- gsub("\\d\\.", "", rownames(disps))

pvals <- testForDEUParallel(
    dxd[gene1Index,],
    formulaNull,
    formulaFull,
    disps,
    testableVector,
    bjp )

data("pvals1")

all( pvalsTissues1[gene1Index] == pvals )