Попробуем применить различные классификаторы на уже знакомом нам датасете про эпигенетические модификации в промоторах и энхансерах.

Загрузим данные об эпигенетических модификациях в энхансерах и промоторах. Для простоты создано две таблицы - обучающая выборка (train) и тестовая (test) выборка. Данные для удобства уже прологарифмированы.

train=read.table("http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter_train.tab", header=T)
test=read.table("http://makarich.fbb.msu.ru/artemov/R/enhancer_promoter_test.tab", header=T)

(Шаг 0). Посмотрим, как устроены данные, методом PCA:

#все колонки, кроме y (столбец Type)
train2=train[,2:ncol(train)]
pc=prcomp(train2)
#покрасим по типа элемента (энхансер и промотер)
plot(pc$x, pch=19, col=train$Type)

Сразу нарисовать графики со всеми главными осями

pairs(pc$x, col=train$Type)
#или
pairs(pc$x[,1:4], col=train$Type)

Какие переменные x вносят наибольший вклад в главные оси?

plot(pc$rotation, pch=19)
text(pc$rotation[,1], pc$rotation[,2], rownames(pc$rotation))

barplot(pc$rotation[,1], las=2, main="PC1")

Как строить разные классификаторы:

SVM

install.packages("e1071")
library(e1071)
model=svm(Type ~ . , data=train)

Пример нелинейного классификатора + рисунок (можно нарисовать для двухмерного случая - два предиктора)

train1=train[,c(1,8,9)]
model=svm(Type ~ H3k4me1 + H3k4me3, data=train1)
plot(model, train1)

Ошибка модели на тестовой выборке

test1<-test[,c(1,8,9)]
mean(test1$Type!=predict(model, test1))

(*) Поэкспериментируйте с выбором ядра (параметр kernel):

model=svm(Type ~ H3k4me1 + H3k4me3, data=train1, kernel="polynomial")
model=svm(Type ~ H3k4me1 + H3k4me3, data=train1, kernel="radial")
model=svm(Type ~ H3k4me1 + H3k4me3, data=train1, kernel="linear")

Попробуем подобрать параметры

tuned <- tune.svm(Type ~., data = train1, gamma = 10^(-4:-1), cost = 10^(0:2)) # tune
summary(tuned)

Random Forest

Чтобы понять, как работает метод, построим одно дерево принятия решений

library(tree)
tr=tree(Type ~ . , data=train)
plot(tr)
text(tr, use.n=TRUE, all=TRUE, cex=.8)

Random Forest строит много деревьев принятия решений, предсказание получается "голосованием". При построении каждого из деревьев из обучающей матрицы с данными выкидываются фрагменты, поэтому деревья получаются разными.

library(randomForest)
model = randomForest(Type ~ ., train)
#важность переменных для классификатора
importance(model)

#вытащим одно из деревьев леса
getTree(model, k=10, labelVar=T)

Нейронные сети

install.packages("neuralnet")
library(neuralnet)
net.promot<- neuralnet(as.numeric(Type=="P")~Control+Ctcf+H2az+H3k27ac+H3k27me3+H3k36me3+H3k4me1+H3k4me3+H3k9me3+Pol2, data=train, hidden=5, threshold=0.01)
plot(net.promot)
net.results <- compute(net.promot, test[,2:11])
ls(net.results)
print(net.results$net.result)

Кросс-валидация (для любого метода)

pr=predict(model, test)
#сравним pr - предсказанный по модели класс и test$Type - настоящий класс элемента генома
table(predicted=pr, real=test$Type)
accuracy=mean(test$Type!=predict(model,test))