Kodomo

Пользователь

Учебная страница курса биоинформатики,
год поступления 2015

Множественное тестирование

Генерим на каждом шаге 20 значений из нормального распределения, считаем pvalue тестом Стьюдента 

set.seed(12345)
pValues <- rep(NA,1000)
for(i in 1:1000)
{
        x <- rnorm(20)
        pValues[i] <- t.test(x)$p.value
}

Считаем количество значимых экспериментов без поправки и с поправками 

sum(pValues < 0.05) # ≈ 0.05*1000
sum(p.adjust(pValues, method="bonferroni") < 0.05)
sum(p.adjust(pValues, method="BH") < 0.05)

Генерим 500 выборок из стандартного нормального распределения и 500 из нормального со средним 1.5. Помечаем результаты, чтобы посчитать TP, FP и т.п. 

set.seed(12345)
pValues <- rep(NA,1000)
for(i in 1:500){pValues[i] <- t.test(rnorm(20))$p.value}
for(i in 501:1000){pValues[i] <- t.test(rnorm(20, mean=1.5))$p.value}
trueStatus <- rep(c("zero", "not zero"), each=500)

Анализируем результат 

table(pValues < 0.05, trueStatus)
table(p.adjust(pValues, method="bonferroni")<0.05, trueStatus)
table(p.adjust(pValues, method="BH")<0.05,trueStatus)

Задание

1) В качестве примера, García-Arenzana et al. (2014) проверили ассоциации из 25 диетических переменных с маммографической плотностью, важным фактором риска развития рака молочной железы, у испанских женщин. Результаты тут: dietary.csv

Какие переменные значимо ассоциированы с маммографической плотностью на уровне значимости 5%? Как изменится значимость после поправки на множественное тестирование?

Кластеризация

Работаем с расстояниями 

dat <- rbind(samp1=c(0,0,0), samp2=c(0,0,3), samp3=c(1,1,1), samp4=c(2,2,2))
dist(dat, method="euclidian")
as.matrix(dist(dat, method="euclidian"))

Иерархическая кластеризация 

y <- matrix(rnorm(50), 10, 5, dimnames=list(paste("g", 1:10, sep=""), paste("t", 1:5, sep=""))) # генерируем данные 
c <- cor(t(y), method="spearman") 
d <- as.dist(1-c)
hr <- hclust(d, method = "complete", members=NULL) 
plot(hr)
plot(as.dendrogram(hr), edgePar=list(col=3, lwd=4), horiz=T)
str(as.dendrogram(hr))
library(ctc); hc2Newick(hr)#скобочная структура

Частота слов в блогах: blogdata.csv 

blogdata<-read.csv("blogdata.csv")
idx <- sample(1:dim(blogdata)[1], 10)
blogdataSample <- blogdata[idx,]
blogdataSample$Blog <- NULL
hc <- hclust(dist(blogdataSample), method="single")
plot(hc, hang=-1, labels=blogdata$Blog[idx])

hc <- hclust(dist(blogdataSample), method="complete")

Обрезаем дерево 

mycl <- cutree(hc, k=3) 
print(mycl)
plot(hc)
abline(h=55, col="red")

Неиерархическая кластеризация. k-means

km <- kmeans(blogdataSample, 4)

DBSCAN

library(dbscan)
data(iris)
iris <- as.matrix(iris[,1:4])
kNNdistplot(iris, k = 5)
abline(h=.5, col = "red", lty=2)

res <- dbscan(iris, eps = .5, minPts = 5)
print(res)
pairs(iris, col = res$cluster + 1L)

Сгенерируем данные, на которых не сработают алгоритмы кластеризации, основанные на измерении расстояния 

r1=runif(200)
alpha1=runif(200, 0, 2*pi)
cl1=data.frame(x=r1*cos(alpha1), y=r1*sin(alpha1))
plot(cl1)


alpha2=runif(200, 0, 2*pi)
r2=rnorm(200, mean=5)
cl2=data.frame(x=r2*cos(alpha2), y=r2*sin(alpha2))
plot(cl2)

cl_both=rbind(cl1, cl2)
colors=c(rep("black", 200), rep("red", 200))

plot(cl_both, col=colors, pch=19)

Применим k-means clustering (k=2) 

km = kmeans(cl_both,    2)
plot(cl_both, col=1+km$cluster, pch=19)

Применими DBSCAN 

install.packages("dbscan")
library(dbscan)

res = dbscan(cl_both, eps=0.3)
plot(cl_both, col=1+res$cluster, pch=19)

Задание

Загрузим данные об эпигенетических модификациях в энхансерах и промоторах. Для простоты создано две таблицы - обучающая выборка (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)

1. Постройте иерархическую кластеризацию участков генома. (Проверьте, что вы не передаете на вход алгоритму кластеризации столбец с типом геномного участка).

2. Постойте k-means кластеризацию с k=2. Получилось ли увидеть кластеры энхансеров и промоторов.

3. Примените DBSCAN к тем же данным.