Predictive model

The functions in this tutorial are a simplified interface to the KDM library. To reproduce the results reported in the paper, use the scripts in CODE/R-SCRIPTS of the GitHub repository, where cross-validation folds are fixed within the datasets. Numbers obtained with this tutorial may therefore differ slightly from those in the paper.

KDM-LRLM turns the discovered motifs into a binding predictor using Lasso-regularized logistic regression. The Lasso penalty drives most motif coefficients to exactly zero, so the fitted model is also a selection of the motifs that matter.

KDM-LRLM is implemented in three distinct functions: KDM.LRLM.train for training the model, KDM.LRLM.performance to evaluate predictive performances, and KDM.LRLM.predict to predict bound sequences.

Model identification

Function: `KDM.LRLM.train`

Parameter	Meaning
`motifs`	a `kdmotifs` object containing the motifs
`pseq`	a set of "positive" sequences (i.e. bound)
`nseq`	a set of "negative" sequences (i.e. unbound)
`win.mod`	a data frame containing the definition of the regions the sequences are split into, in order to obtain the features on which the model is identified. It must contain `start`, the starting position of the region; `end`, the ending position of the region; and `group`, the group the regions belong to — features on regions of the same group are collapsed into a single one (useful for ChIP-seq derived sequences in which distal regions are equivalent)
`nfold`	the number of folds used internally by `cv.glmnet` for cross validation
`seed`	the seed used for cross validation by `cv.glmnet` for reproducibility (default: `100`)
Output	a list containing `motifs`, an object of class `kdmotifs` with the motifs the model retained; `coefficients`, a matrix with the model coefficients organized by columns, each corresponding to a group as defined in `win.mod`; `intercept`, the intercept of the model; and `win.mod`, the input parameter reported back

To use the function we must first define the model windows. For transcription factors we used three windows: one centered around the midpoint of the sequence (which, in our dataset, corresponds to the summit of the ChIP-seq peaks), and two distal whose features are merged together. So we can define:

model_windows<-data.frame(start=c(1,91,181),end=c(121,211,301),group=c("distal","central","distal"))

At this point we can call the training function using the training sequences:

model<-KDM.LRLM.train(motifs=motifs,pseq=NRF1$train$pseq,nseq=NRF1$train$nseq,win.mod=model_windows,nfold=10)

We have:

model$motifs

kdm motifs matrix
 #N motifs: 58
 Kmer length: 8
 N# gaps: 2
 Double strand: TRUE

head(model$coefficients)

       central   distal
M-1   4.9847624 0.000000
M-3   1.4497146 0.000000
M-12  0.0000000 8.198660
M-15  5.3623925 3.851842
M-16 -0.6610866 0.000000
M-19  1.5381432 0.000000
. . .

model$intercept

[1] -4.825192

Model evaluation

Predictive performance measures are obtained through the KDM.LRLM.performance function:

Function: `KDM.LRLM.performance`

Parameter	Meaning
`model`	the model returned by `KDM.LRLM.train`
`pseq`	a set of "positive" sequences (i.e. bound)
`nseq`	a set of "negative" sequences (i.e. unbound)
Output	a named numeric array: `AUROC`, area under the ROC curve; `AUPRC`, area under the Precision-Recall curve; and `ACCURACY`, the highest accuracy over all possible score thresholds

In our example we use the test set sequences for evaluating the predictive performance:

performances<-KDM.LRLM.performance(model,NRF1$test$pseq,NRF1$test$nseq)
performances

   AUROC     AUPRC  ACCURACY
0.9920480 0.9812351 0.9780000

Model predictions

🚧 To be completed. Add the worked prediction example — take a long sequence (e.g. a promoter), split it into 301-bp sub-sequences, and show the KDM.LRLM.predict call and how to read the resulting profile.

KDM — source on GitHub · preprint