Hi, 

I have a three groups with 20-30 samples each, and I'd like to generate predictive gene signatures for each. Do any of you have suggestions on R packages that you've used in the past and that I can look into to generate these types of signatures?

Kindly,
Martin

If you have already derived a panel of genes that are statistically significantly differentially expressed, then there are various tools to use.

# continuous outcome

If you are modelling a ***continuous*** trait (e.g. blood levels of inflammatory markers, metabolite levels, etc.), then build a linear model using `lm()`:

For example:

    lm(MyContinuousTrait ~ gene1 + gene2 + gene3 + ... + geneX, data=MyData)

You can then do stepwise regression to refine the model further based on the AIC or BIC using `stepAIC()` from *MASS* package.

If you have many genes (e.g. > 30), then first test each gene independently, followed by building a combined model (like above) with the statistically significant genes:

    lm(MyContinuousTrait ~ gene1, data=MyData) # p=0.8
    lm(MyContinuousTrait ~ gene2, data=MyData) # p=0.03
    lm(MyContinuousTrait ~ gene3, data=MyData) # p=0.67
    lm(MyContinuousTrait ~ gene4, data=MyData) # p=0.02

    MyModel <- lm(MyContinuousTrait ~ gene2 + gene4, data=MyData)

With this final model, you can still do stepwise regression to refine it even further.

For performing QC on this model::

    library(car)

    # check for outliers based on Cook's distances (2 different plots)
    plot(cooks.distance(MyModel))
    cutoff <- 4 / ((nrow(MyData) - length(MyModel$coefficients) - 1)) 
    plot(MyModel, which=4, cook.levels=cutoff)

    # Bonferonni p-value for most extreme observations
    outlierTest(MyModel)

    # check for non-normality via QQ plot for 'studentised' residuals
    qqPlot(MyModel, main="QQ Plot")

    # evaluate homoscedasticity
    # non-constant error variance test
    ncvTest(MyModel)

    # multi-collinearity via variance inflation factors 
    vif(MyModel)
    # a VIF > 2.5 indicates that the R^2 of a predictor to all other predictors is 0.6, i.e., collinear
    sqrt(vif(MyModel)) > 2.5

    # non-linearity via component + residual plot
    crPlots(MyModel)

    # test for Autocorrelated Errors
    durbinWatsonTest(MyModel)

    # global test of model assumptions
    library(gvlma)
    gvmodel <- gvlma(MyModel) 
    summary(gvmodel)

    # k-fold cross-validation (here 3-fold)
    library(DAAG)
    cv.lm(MyData[!is.na(MyData$MyOutcomeTrait),], MyModel, m=3)

    # R^2 shrinkage using 10-fold cross-validation
    library(bootstrap)
    theta.fit <- function(x, y) { lsfit(x,y) }
    theta.predict <- function(fit, x) { cbind(1,x) %*% fit$coef }
    x <- as.matrix(MyData[,which(colnames(MyData) %in% names(summary(MyModel)$coefficients[,1]))])
    y <- as.matrix(MyData[,c("MyOutcomeTrait")]) 
    results <- crossval(x, y, theta.fit, theta.predict, ngroup=10)
    # Raw R^2 
    cor(y, MyModel$fitted.values)**2
    # cross-validated R^2
    cor(y, results$cv.fit)**2

Finally, if you have built a linear model, you can make a nice plot with *ggplot2* using the following code:

    MyData$MyContinuousTrait.pred = predict(MyModel)
    ggplot(MyData, aes(x=MyContinuousTrait, y=MyContinuousTrait.pred)) +
      geom_point() + xlim(0,5) + ylim(0,5) +
      geom_smooth(method='lm',formula=y~x) +
      stat_smooth(method="lm", fullrange=TRUE) +
      ggtitle("")

<a href="https://imgbb.com/"><img src="https://image.ibb.co/f8aSw5/Captura_de_tela_de_2017_09_20_18_01_38.png" alt="Captura_de_tela_de_2017_09_20_18_01_38" border="0"></a>

 - X-axis:  real values for the outcome variable
 - Y-axis:  model predicted values

# categorical outcome

If your outcome is ***categorical***, i.e., a binary trait or has greater than 2 categories, then first build a binary or multinomial logistic regression model using `glm(..., family="binomial")` with your predictors (genes). You can then choose the best predictors again based on the AIC or BIC using stepwise regression.

The checks for a `glm()` differ from those for a `lm()`:

    # McFadden R^2 index used to assess the model fit
    library(pscl)
    pR2(MyModel)

    # perform Chi^2 ANOVA
    anova(MyModel, test="Chisq")

    # perform cross validation
    # the delta values should not greatly differ
    # K=number of samples, i.e., leave-one-out CV.
    library(boot)
    cv.glm(MyData, MyModel, K=nrow(MyData))$delta

---------
-

For both `lm()` and `glm()` - based models, after you have chosen your best predictors and defined a final model, you can derive odds ratios (with confidence intervals) with:

    # get ORs and confidence intervals
    exp(cbind(OR=coef(MyModel), confint(MyModel)))

You can also perform ROC analysis using the *pROC* package.

Other options that can be used instead of `lm()`/`glm()` are lasso-penalized, elastic-net, and ridge regression (*glmnet* package).
