Занятие 4 if, for, apply, merge тесты ассоциации 30 сентября 2015. - презентация

1 Занятие 4 if, for, apply, merge тесты ассоциации 30 сентября 2015

2 Биологическая задача Найти аллельные варианты, которые отличают популяцию европейцев и африканцев Для каждого человека можем узнать: К какой расе принадлежит? Присутствует ли данный аллельный вариант* в данном месте последовательности ДНК? A G CA TA присутствие аллельного варианта («мутация») *в гетеро- или гомозиготном состоянии 2

3 Биологическая задача Задача: какие генетические варианты (мутации) отличают популяцию европейцев и африканцев Исходные данные: > variant_counts chrposEA_refEA_altAA_refAA_alt chr chr Хромосома Координата Европейцы Африканцы Кол-во людей без альт. варианта Кол-во людей с альтернативным вариантом … variant_counts.csv 3

4 Биологическая задача Второй шаг – понять, в каких генах эти расположены наиболее значимые мутации A G CA TA gene1 chrposgenes chr DUSP22 chr DUSP22 Хромосома КоординатаГены … >variant_genes variant_genes.csv 4

5 Что нужно для решения такой задачи? Какая нужна статистика? Проверка гипотезы об ассоциации. Точный тест Фишера, Хи^2 Как посчитать что-то для всех строк таблицы? Условия и циклы (if, for) Функции, применение их к элементам вектора, рядам и колонкам матрицы Как объединить две таблицы? SQL-подобные операции с таблицами: merge, aggregate 5

6 ПРОВЕРКА ГИПОТЕЗЫ ОБ АССОЦИАЦИИ Точный тест Фишера Хи^2 6

7 Тесты ассоциации У некоторых объектов есть два свойства, про которые, в простейшем случае, можно сказать «да» или «нет», например, больной – здоровый, носитель мутации – не носитель, курит – не курит. Связаны ли эти свойства? Примеры задач: Ассоциирована ли определенная мутация с болезнью Верно ли, что на биофаке значимо большая доля девушек, чем на ВМК? 7

rownames(gwas)=c("Healthy", "Diseased") > gwas refe" title="Тесты ассоциации. Мутация, ассоциированная с болезнью > gwas=matrix(c(3324,1896,2676,2104), nrow=2, ncol=2) > gwas [,1] [,2] [1,] 3324 2676 [2,] 1896 2104 > colnames(gwas)=c("reference", "mutant") > rownames(gwas)=c("Healthy", "Diseased") > gwas refe" class="link_thumb"> 8 Тесты ассоциации. Мутация, ассоциированная с болезнью > gwas=matrix(c(3324,1896,2676,2104), nrow=2, ncol=2) > gwas [,1] [,2] [1,] [2,] > colnames(gwas)=c("reference", "mutant") > rownames(gwas)=c("Healthy", "Diseased") > gwas reference mutant Healthy Diseased Таблица сопряженности Носитель мутации* Не-носитель мутации * Пояснение «для биологов»: подразумеваем, например, что носитель мутации = носитель болезни (т.е. мутантный аллель - доминантный) Здоровый (контроль) Больной 8 rownames(gwas)=c("Healthy", "Diseased") > gwas refe"> rownames(gwas)=c("Healthy", "Diseased") > gwas reference mutant Healthy 3324 2676 Diseased 1896 2104 Таблица сопряженности Носитель мутации* Не-носитель мутации * Пояснение «для биологов»: подразумеваем, например, что носитель мутации = носитель болезни (т.е. мутантный аллель - доминантный) Здоровый (контроль) Больной 8"> rownames(gwas)=c("Healthy", "Diseased") > gwas refe" title="Тесты ассоциации. Мутация, ассоциированная с болезнью > gwas=matrix(c(3324,1896,2676,2104), nrow=2, ncol=2) > gwas [,1] [,2] [1,] 3324 2676 [2,] 1896 2104 > colnames(gwas)=c("reference", "mutant") > rownames(gwas)=c("Healthy", "Diseased") > gwas refe">

9 Тест Фишера и Хи^2 > fisher.test(gwas ) Fisher's Exact Test for Count Data data: gwas p-value = 4.874e-15 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: sample estimates: odds ratio > chisq.test(gwas) Pearson's Chi-squared test with Yates' continuity correction data: gwas X-squared = , df = 1, p-value = 5.055e-15 9

10 Тест Фишера и Хи^2 Другой (иногда более удобный) способ – передавать не четырехпольные таблицы, а два логических вектора, соответствующие двум свойствам, между которыми ищем ассоциацию > people=data.frame(mutant=c(T,F,F,F,T,T,F,T), diseased=c(T,F,F,F,F,T,T,T)) > fisher.test(people$mutant, people$diseased) Удобный способ получить 4-польную таблицу: > xtabs(~ mutant + diseased, data=people) diseased Mutant FALSE TRUE FALSE 3 1 TRUE

11 Точный тест Фишера не мутантмутантВСЕГО здоровыйaba+b больнойcdc+d ВСЕГОa+cb+dn Точный тест Фишера (Fisher exact test): Сколько способов разложить n независимых элементов по 4 ячейкам с данным количеством элементов? p-value – сумма вероятностей таблиц с не меньшим OR. Отношение рисков (odds ratio) 11

12 Тест Хи^2 (Chi-squared) Не рекомендуется использовать, если хотя бы в одной ячейке меньше 5-7 образцов H0: Предположение независимости. Например, p(больной & мутант) = p(больной)*p(мутант) Тест Хи^2: Посчитаем значение в каждой ячейке в предположении независимости, это будет ожидаемое значение. Далее по всем ячейкам суммируем 12 d.f.=(nrow-1)(ncol-1)

13 «ПРОГРАММИРОВАНИЕ» Условия и циклы Функции и apply 13

source(my_" title="Цикл for Пример из домашнего задания. N_iter=20 y=rep(NA, N_iter) for(i in 1:N_iter){ rands=rnorm(2^i) y[i]=max(rands) } plot(2^(1:N_iter), y, type="l") Совет: многострочные исходные коды удобно сохранять в файл, читать и выполнять из R: > source(my_" class="link_thumb"> 14 Цикл for Пример из домашнего задания. N_iter=20 y=rep(NA, N_iter) for(i in 1:N_iter){ rands=rnorm(2^i) y[i]=max(rands) } plot(2^(1:N_iter), y, type="l") Совет: многострочные исходные коды удобно сохранять в файл, читать и выполнять из R: > source(my_program.R") Вектор (в данном случае от 1 до N_iter) i проходит все элементы вектора Для каждого i выполняется код в { } 14 source(my_"> source(my_program.R") Вектор (в данном случае от 1 до N_iter) i проходит все элементы вектора Для каждого i выполняется код в { } 14"> source(my_" title="Цикл for Пример из домашнего задания. N_iter=20 y=rep(NA, N_iter) for(i in 1:N_iter){ rands=rnorm(2^i) y[i]=max(rands) } plot(2^(1:N_iter), y, type="l") Совет: многострочные исходные коды удобно сохранять в файл, читать и выполнять из R: > source(my_">

15 for. Пример для таблицы > head(grades) id write math science socst #посчитаем средний балл каждого студента и средний балл за каждый экзамен > matrix_grades=as.matrix(grades[,2:ncol(grades)]) > meangrades=c() > for(i in 1:nrow(grades)){ + meangrades=c(meangrades, mean(matrix_grades[i,])) + } > head(meangrades) [1]

16 apply. Пример для таблицы #средний балл каждого студента > meangrades=apply(grades[,2:ncol(grades)], 1, mean) > head(meangrades) [1] #средний балл по каждому экзамену > examgrades=apply(grades[,2:ncol(grades)], 2, mean) > head(examgrades) write math science socst Применить функцию к каждому ряду или колонке data frame или элементу вектора mean( c(52,41,47,57) ) mean( c(59,53,63,61) ) mean( c(33,54,58,31) ) Матрица или data frame Измерение (1-строки или 2-столбцы) Функция 16

17 Функции – использование Имя функции со скобками – её вызов, в скобках аргументы. Функция может возвращать значение. tt=t.test(a, b, paired=T) Имя функции без скобок можно использовать как имя переменной – посмотреть содержимое или передать как аргумент другой функции > apply function (X, MARGIN, FUN,...) { … Помимо стандартных функций (например, mean, all, any) можно передать apply свою функцию Именованный аргумент 17

18 Функции – создание Создадим функцию, которая берет на вход вектор из 4 чисел, создает четырехпольную таблицу, выполняет тест Фишера или Хи^2 association.test <- function(v){ m=matrix(v, nrow=2, ncol=2) #создаем матрицу res=fisher.test(m) return(res$p.value) ИЛИ res$p.value } обычное присваивание Всегда возвращается последнее вычисленное значение, даже если не писать return 18

19 Функции – создание. if и аргумент по умолчанию association.test <- function(v, test="fisher"){ m=matrix(v, nrow=2, ncol=2) #создаем матрицу if ( test=="fisher" ){ res=fisher.test(m) } else if( test== "chi2" ) { res=chisq.test(m) } else { return (NA) } return(res$p.value) ИЛИ res$p.value } Значение по умолчанию условие 19

20 Другие разновидности apply Если на входе список (list), то lapply, sapply > l=list(1, c(1,2,3), c(2,1)) > lapply(l, sum) [[1]] [1] 1 [[2]] [1] 6 [[3]] [1] 3 > sapply(l, sum) [1] #список, каждый элемент – результат работы функции #Пытается «упростить» список. В данном случае свёл к вектору 20 См. также unlist

21 Другие разновидности apply Более содержательные пример – а пусть функция выдает не один элемент, а вектор из двух > lapply(l, function(x){return(c(min(x), max(x)))}) [[1]] [1] 1 1 [[2]] [1] 1 3 [[3]] [1] 1 2 > sapply(l, function(x){return(c(min(x), max(x)))}) [,1] [,2] [,3] [1,] [2,] min max min max #В данном случае «упростил» до матрицы 21

22 Вернемся к мутациям Уже можем посчитать для каждой позиции (строки в таблице) p-value ассоциации соответствующего варианта с расой apply(___, ___, association.test, test="fisher") Подсказка: все аргументы, которые apply не распознала как свои, передаются вызываемой функции Осталось понять, каким генам соответствуют мутации (сделать из двух таблиц одну) 22

23 SQL-ПОДОБНЫЕ ОПЕРАЦИИ С ТАБЛИЦАМИ merge aggregate, tapply 23

24 merge Студент ФакультетКурс Пол Иванов Петров СтудентКР1КР2КР3 Иванов Сидоров Студент ФакультетКурс ПолКР1КР2КР3 Иванов 24

25 merge Пример: 2 таблицы id -> оценки, id -> метаданные студента Объединить таблицы в одну по идентификатору студента 25 > head(grades) id write math science socst > head(metadata) id female race schtyp prog

26 merge > head(grades) id write math science socst > head(metadata) id female race schtyp prog > head(grades) id write math science socst > head(metadata) id female race schtyp prog > meta_grades=merge(metadata, grades, by.x=c(1), by.y=c(1)) > head(meta_grades) id female race schtyp prog write math science socst > meta_grades=merge(metadata, grades, by=c("id")) > head(meta_grades) id female race schtyp prog write math science socst > meta_grades=merge(metadata, grades, by=c("id"), all.x=F, all.y=T) В итоговой табл – все эл-ты второй таблицы. Там, где не нашлось пары, заполняется NA by может содержать несколько элементов 26 head(meta_grades) id female race schtyp prog write math science socst 1 1 1 1 1 3 44 40 39 41 2 2 1 1 1 3 41 33 42 41 3 3 0 1 1 2 65 48 63 56 4 4 1 1 1 2 50 41 39 51 5 5 0 1 1 2 40 43 45 31 6 6 1 1 1 2 41 46 40 41 > meta_grades=merge(metadata, grades, by=c("id"), all.x=F, all.y=T) В итоговой табл – все эл-ты второй таблицы. Там, где не нашлось пары, заполняется NA by может содержать несколько элементов 26">

27 aggregate Объединение рядов с одинаковыми значениями определенных столбцов id студента ФакультетОценка 1ФББ ВМК Биофак Group Оценка ФББ4.5 ВМК4.67 Биофак 4.5 mean(c(5,5,4)) 27

28 aggregate > mean_race=aggregate(meta_grades[,6:9], by=list(meta_grades$race), FUN=mean) > mean_race Group.1 write math science socst > mean_ gender =aggregate(meta_grades[,6:9], by=list(meta_grades$female), FUN=mean) > mean_gender Group.1 write math science socst

29 aggregate > mean_gender_race=aggregate(meta_grades[,6:9], by=list(meta_grades$female, meta_grades$race), FUN=mean) > mean_gender_race Group.1 Group.2 write math science socst

30 tapply Похожая идея, но для векторов На входе – вектор (а не таблица) и список (длины N) из векторов факторовов той же длины На выходе – N-мерный вектор > mean_race=tapply(meta_grades$math, list(meta_grades$race), FUN=mean) > mean_race > mean_race=tapply(meta_grades$math, list(meta_grades$race, meta_grades$female), FUN=mean) > mean_race раса Раса Пол 30

31 Ещё раз про домашнее задание 31 chrposEA_refEA_altAA_refAA_alt chr chr … variant_counts.csv … P-value Fisher P-value χ2 Odds ratio Для каждой строки посчитать: добавить как столбцы к data frame chrposgenes chr DUSP22 chr DUSP22 … Объединить таблицы по 2 общим полям …gene…P-value DUSP221e-10 DUSP221e-5 HLA-A1e-13 HLA-A1e-10 …gene…P-value DUSP221e-10 HLA-A1e-13 Группировать по генам, взяв min p-value variant_genes.csv