Исходные данные. в аннотированных энхансерах и промоторах в геноме вычислены уровни различных эпигенетических модификаций: модификации гистонов H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3, H3K9me3, ChIP-seq на РНК-полимеразу II (Pol2) и CTCF.

В таблице http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter.tab приведены уровни эпигенетических сигналов, а первый столбец Type указывает на то, является отрезок энхансером (1) или промотором (0).

Задача: научиться отличать по эпигенетическим модификациям промоторы от энхансеров.

tab2=read.table("http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter.tab", header=T, stringsAsFactors=F)

Визуализация данных

library(gplots)

#Корреляции между предикторами, а также между предикторами и outcome
heatmap.2(cor( tab2 ), col="bluered", symbreaks=T)

heatmap: возьмем 100 случайных промоторов и 100 случайных энхансеров:

heatmap.2(data.matrix(tab2[c(sample(which(tab2$Type==0), 100), sample(which(tab2$Type==1), 100)),2:ncol(tab2)]), trace="none", col="bluered", scale="column", breaks= (-20:20)/10, Rowv=F, RowSideColors=c(rep("red", 100), rep("orange", 100)))

#более детально:
#вектор номеров строк, первые 100 указывают на строки с промоторами, последние 100 - на строки с энхансерами
#c(sample(which(tab2$Type==0), 100), sample(which(tab2$Type==1), 100))
#Не перемешиваем по строкам (уже упорядочены так, как нам надо): Rowv=F, 
#Цвет сбоку от heatmap-а указывает, является элемент промотором (красный) или энхансером (оранжевый):
#RowSideColors=c(rep("red", 100), rep("orange", 100))

Построение и сравнение обобщенных линейных моделей (логистическая регрессия)

Для начала, проведем логистическую регрессию от одной переменной, например, H3k4me3. По тем коэффициентам регрессии, которые получаются после обучения, модель строит линейную комбинацию предикторов - число от -Inf до Inf. Затем данное число нужно превратить в вероятность, что Type==1 (то есть объект - энхансер, а не промотер). Такие вероятности хранятся в векторе gl$fitted (для тех точек, на которых модель была обучена).

gl=glm(Type ~ H3k4me3, family=binomial, data=tab2)
subsamp=sample(1:nrow(tab2), 200)
plot(tab2[subsamp,]$H3k4me3, tab2[subsamp,]$Type, pch=19, col='blue', xlab="H3K4me3", ylab="Probability of being an enhancer")
points(tab2[subsamp,]$H3k4me3, gl$fitted[subsamp], pch=19, col='red')
legend("topright", legend=c("Real (0 or 1)", "Predicted by model"), fill=c("blue", "red"))

Модель с двумя предикторами (H3k4me1 + H3k4me3)

gl=glm(Type ~ H3k4me1 + H3k4me3, family=binomial, data=tab2)
anova(gl2, gl1, test="Chisq")

Сравнение моделей

gl1=glm(Type ~ ., family=binomial, data=tab2)
gl2=glm(Type ~ H3k4me1 + H3k4me3, family=binomial, data=tab2)
anova(gl2, gl1)
anova(gl2, gl1, test="Chisq")

Кросс-валидация.

cost=function( y, ypred ) {
  mean( abs(y-ypred)>0.5 )
}

library(boot)
cv1=cv.glm(tab2, gl, cost, K=5)
cv1$delta

Примените кросс-валидацию к другим возможным моделям: например, к моделям с одним предиктором или к модели, учитывающей все предикторы (Type ~ .).

Можно ли без ущерба оставить только два предиктора?

(*) Автоматический выбор параметров

library(glmulti)
gl=glmulti(Type ~ ., family=binomial, data=tab2, level=1)

To be continued (SVM и другие классификаторы).