Учебная страница курса биоинформатики,
год поступления 2015
Классификаторы
Попробуем применить различные классификаторы на уже знакомом нам датасете про эпигенетические модификации в промоторах и энхансерах.
Загрузим данные об эпигенетических модификациях в энхансерах и промоторах. Для простоты создано две таблицы - обучающая выборка (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)
Валидация
pr=predict(model, test) #сравним pr - предсказанный по модели класс и test$Type - настоящий класс элемента генома table(predicted=pr, real=test$Type) accuracy=mean(test$Type!=predict(model,test))
Вопрос: какое количество деревьев стоит по-умолчанию? Попробуйте разные варианты данного параметра.