Задача прогнозирования высоты AI openSNP
[изначально опубликовано ЗДЕСЬ]
CrowdAI - это ОЧЕНЬ крутой сайт, на котором проводятся соревнования по машинному обучению. Недавно они провели соревнование с использованием данных с openSNP, веб-сайта, который позволяет любому загружать свои генетические данные, делая их общедоступными. Около 3500 человек загрузили некоторое количество своих генетических данных (многие просто загрузили митохондриальную часть), и 921 человек также сообщили о своем росте. Соревнование CrowdAI было очень простым: можете ли вы предсказать рост в подмножестве этих людей, используя генетические данные от всех остальных для обучения модели?
А также? "Я выиграл!!!" Моя лучшая модель предсказывает 53,45% отклонения в росте, в то время как следующая лучшая модель предсказывает 48,62%. Но есть одна загвоздка: хотя это было соревнование по машинному обучению, я на самом деле не использовал машинное обучение ... ой
Этот пост предназначен для того, чтобы объяснить, почему я участвовал и как я выиграл, в надежде, что другие, обладающие большим опытом в области машинного обучения, смогут извлечь выгоду из моих знаний в области генетики.
Почему я участвовал
Я аспирант, изучаю нейробиологию и генетику. Это именно та проблема, над которой я размышляю целыми днями. В моих исследованиях чаще всего больше: Предсказывает ли генетический риск данного расстройства, насколько часто люди используют определенную часть своего мозга во время психологической задачи?, Но лежащие в основе методы остаются теми же. Фактически, я недавно был одним из преподавателей летнего семинара по применению методов исследования поведенческой генетики в психологии и нейробиологии. На самом деле я впервые нашел конкурс CrowdAI, потому что искал демо-данные, которые можно было бы передать участникам семинара.
Я вспомнил о конкурсе через месяц после окончания семинара, и когда я вернулся к нему, я увидел, что было много заявок, которые предсказывали только ~ 40% дисперсии в подмножестве окончательного набора данных тестирования (14 человек). Мне показалось, что это очень близко к тому, сколько отклонений можно ожидать от биологического пола человека.
Итак, я решил представить свое предположение исследователя. НЕ, потому что я хотел отказаться от машинного обучения, а потому, что мне нужен был эталонный тест, с которым можно было бы сравнить отправленные материалы по ML. Я полностью ожидал, что кто-то превзойдет мою догадку. Но при машинном обучении важно убедиться, что базовая модель является чем-то разумным, о чем говорится в этой недавней статье. Затем я попытался улучшить это с помощью собственного машинного обучения, но безуспешно.
На самом деле, вся эта ситуация напоминает мне Соревнование СДВГ-200, которое должно было правильно классифицировать пациентов с СДВГ на основе их состояния покоя, связанного с фМРТ-сетями мозга. Оказывается, победителю вообще не нужно было использовать какие-либо данные мозга!
Не поймите меня неправильно, меня очень интересует машинное обучение. Я думаю, что у него есть большой потенциал, чтобы помочь в изучении генетики человека. Эта статья доктора Кристин Никодемус, в которой она использует случайный лес для выявления нелинейного взаимодействия между тремя мутациями, повышающими риск шизофрении, является одним из моих любимых исследований. Это также было одним из основных источников вдохновения для моей заявки на получение стипендии для аспирантов NSF (спасибо NSF за финансирование!), В которой предлагалось объединить машинное обучение с генетикой и нейровизуализацией. Доктор Никодимус тоже не отказалась от машинного обучения - вот слайды недавнего выступления на эту тему. Также см. Эту недавнюю статью, где приведен пример применения машинного обучения для изучения роста. Но если вы посмотрите на обсуждение в твиттере вокруг него, вы заметите, что тот факт, что они не включали подход, который я использовал для этого конкурса, считается красным флагом .
Я думаю, что этот опыт также может послужить уроком для организаторов соревнований по ML в целом. Используя только свои знания в области генетики, я смог сгенерировать модель, которая представляла собой просто описание набора данных - пол + 3 основных компонента - который работал так же хорошо, как и среднее представление машинного обучения. Это то, что организаторы могут предоставить участникам конкурса, особенно если экспертные знания в определенной области не являются обязательным условием.
Движение вперед
Частично моя мотивация для публикации этого сообщения состоит в том, что я хотел бы дать другим участникам конкурса (и всем, кто заинтересован) шанс превзойти этот базовый уровень. Вы можете найти файл с пятью переменными в моей последней регрессии - да, буквально только с пятью переменными в простой линейной регрессии, кто знает, сколько нейронных сетей - ЗДЕСЬ.
Анализ - Краткое содержание
Модель, которую я построил, которая в итоге закончилась победой, состояла из 5 основных шагов:
- Настроить данные для анализа
- Угадайте биологический пол людей (не было)
Секс может очень хорошо определять рост. Он находится в генах (Y&X Chromosome), но данные о конкурсе не включали его. К счастью, это довольно легко вменить.
- Удалите связанных людей из данных обучения
Как правило, для генетического анализа требуется, чтобы степень родства между участниками была относительно постоянной. Часто это просто то, что все одинаково не связаны. Как вариант, у вас могут быть семьи, где есть подгруппы людей, которые гораздо больше связаны друг с другом, чем с кем-либо еще. В случае тренировочных данных для этого соревнования, то же самое. Участники, как правило, совершенно не связаны между собой, есть небольшая горстка людей, являющихся прямыми родственниками. Это может быть проблемой, потому что это приведет к смещению вашей модели в сторону способности предсказывать только рост людей в этих семьях, вместо того, чтобы предсказывать рост кого-либо.
- Создать основные компоненты
Данные openSNP многонациональны. То есть большинство участников белые, но есть также немало людей африканского и азиатского происхождения. Это проблема, потому что это практически гарантирует увеличение количества ложных срабатываний. Анализ основных компонентов - это стандартный базовый способ: 1) определить, что ваша выборка многоэтническая (в генетике мы называем это стратификацией населения), и 2) произвести начальную, но не полную или полностью достаточную поправку.
- Создайте полигенную оценку риска (PRS) для роста
Полигенные оценки риска заслуживают отдельной публикации. Это стандартный метод, который исследователи поведенческой генетики применяют при попытке предсказать признак в независимой выборке. Шкалы полигенных рисков (PRS) - не новый подход. Я не совсем уверен, какая газета использовала их первой, но точно в начале 2000-х. Они чрезвычайно просты, поэтому замечательно, что они оказались столь полезными. PRS начинается с независимого исследования, GWAS, где каждый аллель в геноме независимо проверяется, чтобы увидеть, предсказывает ли он какой-либо признак. Затем строится PRS путем взятия новой выборки, и для каждого участника в этой выборке складываются все имеющиеся у них аллели риска, взвешенные по размеру эффекта аллеля в первой выборке. Для PRS вы чаще всего включаете в выборку обнаружения только те SNP, которые ниже определенного порогового значения p. Вы можете подумать, что достаточно жесткий порог, скажем, только те SNP, p-значение которых выдерживает коррекцию Бонферрони для нескольких тестов, будет работать лучше всего. Одно из наиболее удивительных наблюдений в результате анализа PRS состоит в том, что обычно очень порог, скажем, p
При шизофрении подход PRS теперь может охватывать более 20% риска заболевания. Вот GWAS, опубликованный в июле этого года, в котором PRS использовалась в проверочной выборке. Основная оговорка в отношении подхода PRS (и, честно говоря, со всеми исследованиями в области генетики человека) заключается в том, что он работает лучше всего, когда все делается в пределах одной и той же этнической выборки. Многоэтнические образцы или образцы совершенно иных этнических групп, чем те, которые использовались в открытии GWAS, будут иметь худшие / плохие результаты.
В этом соревновании каждой команде был предоставлен тщательно подобранный набор из ~ 900 аллелей, которые, как известно, связаны с ростом. Эти аллели были идентифицированы в 2014 GWAS of Height. Организаторы crowAI проделали тяжелую работу по выбору переменных для каждого случая заранее. Теперь результаты GWAS публично доступны здесь. Я просто загрузил их и использовал для создания PRS для роста.
Анализ в деталях
Примечание о программном обеспечении
Я в основном проводил свой анализ в plink и R. Почти весь приведенный ниже код является R-кодом - даже plink-код написан как системный вызов через R. R - мое предпочтительное программное обеспечение (когда-нибудь скоро я Возможно, я столкнусь с проблемой, для которой мне действительно нужно знать Python, но этого еще не произошло), и вызов plink через R удобен, потому что он упрощает создание сценариев и циклов. Кроме того, я пишу это с помощью RMarkdown.
1) Настройка данных для анализа
Первый шаг - преобразовать файлы .vcf, заархивированные с помощью g-zip, предоставленные crowAI, в двоичный формат plink. ничего проблемного.
system(paste ("plink --vcf genotyping_data_fullset_test.vcf --biallelic-only strict --make-bed --out OpenSNP_test",sep=" "))
system(paste ("plink --vcf genotyping_data_fullset_train.vcf --biallelic-only strict --make-bed --out OpenSNP_train",sep=" "))Файлы VCF могут вызывать затруднения из-за большого количества переменных с именами «.» и другие ошибки. Для соревнований я решил проигнорировать их все и просто выбросить все проблемные.
library(data.table) Train_bim<-fread(c("OpenSNP_train.bim")) Train_bim2<-filter(Train_bim, V2 != ".") Train_bim4<-Train_bim2[!duplicated(Train_bim2[c('V2') ]), ]write.table(x=Train_bim4$V2,file=c("SNPS_keep.txt"),sep="",quote = F,col.names = F,row.names = F) rm(Train_bim, Train_bim4, Train_bim2)system(paste ("plink --bfile OpenSNP_train --extract SNPS_keep.txt --make-bed --out OpenSNP_train2",sep=" ")) system(paste ("plink --bfile OpenSNP_test --extract SNPS_keep.txt --make-bed --out OpenSNP_test2",sep=" "))
Теперь мы объединяем их и отбрасываем все остальное, что проблематично (верно, это не лучшая практика для исследования, но она работает для целей конкурса).
system(paste ("plink --bfile OpenSNP_train2 --bmerge OpenSNP_test2 --make-bed --out OpenSNP_combined",sep=" "))system(paste ("plink --bfile OpenSNP_train2 --exclude OpenSNP_combined-merge.missnp --make-bed --out OpenSNP_train3",sep=" ")) system(paste ("plink --bfile OpenSNP_test2 --exclude OpenSNP_combined-merge.missnp --make-bed --out OpenSNP_test3",sep=" "))system(paste ("plink --bfile OpenSNP_train3 --bmerge OpenSNP_test3 --make-bed --out OpenSNP_combined2",sep=" "))
Наконец, перехватите остальные ошибки и удалите их тоже (NB, это не много SNP)
Errors<-readLines(c("OpenSNP_combined2.log")) Rem<-as.data.frame(matrix(NA,length(Errors),1)) for (i in 1: length(Errors)){ Mult<-regexpr('Warning: Multiple positions seen for variant',Errors[i]) if(Mult>0){Rem[i,1]<-substr(Errors[i],47,(nchar(Errors[i])-2)) } Dup<-regexpr('Warning: Variants ',Errors[i]) if(Dup>0){ And<-regexpr(pattern = 'and',Errors[i]) Rem[i,1]<-substr(Errors[i],20,(And[1]-3)) } }Rem<-na.omit(Rem)write.table(x=Rem$V1,file=c("SNPS_remove.txt"),sep="",quote = F,col.names = F,row.names = F) system(paste ("plink --bfile OpenSNP_combined2 --exclude SNPS_remove.txt --make-bed --out OpenSNP_combined3",sep=" "))
2) Угадай биологический секс
Теперь, когда у нас есть все в одном файле plink, мы можем безошибочно угадать их биологический пол. Это делается с помощью флага plink –impute-sex, который вычисляет коэффициенты инбридинга X-хромосомы. Сначала нам нужно отсечь коррелированные snps (см. Документацию plink).
system(paste ("plink --bfile OpenSNP_combined3 --indep-pairphase 20000 2000 0.5 --make-bed --out OpenSNP_combined4",sep=" "))
system(paste ("plink --bfile OpenSNP_combined3 --extract OpenSNP_combined4.prune.in --make-bed --out OpenSNP_combined5",sep=" "))
system(paste ("plink --bfile OpenSNP_combined5 --impute-sex ycount --make-bed --out OpenSNP_combined6",sep=" "))Давайте взглянем на этот файл
library(data.table)
Sexcheck<-fread("OpenSNP_combined6.sexcheck")
hist(Sexcheck$F)
table(Sexcheck$SNPSEX)## ## 0 1 2 ## 1 524 396
Мы видим хорошее разделение оценок F, указывающее на то, что это, вероятно, довольно точно. Однако есть один человек, в котором алгоритм не уверен. Они вероятно биологически женского пола (F = 0,2), но я все равно уберу их. Тогда мы сможем изготовить первую модель!
library(dplyr, quietly = T,warn.conflicts = F,verbose = F)Height<-fread('phenotypes.txt') #This is the training data Heights Sex2<-merge(Sexcheck,Height, by="IID") Sex2<-filter(Sex2, SNPSEX > 0) fit1<-lm(Height ~ SNPSEX, data = Sex2) summary(fit1)## ## Call: ## lm(formula = Height ~ SNPSEX, data = Sex2) ## ## Residuals: ## Min 1Q Median 3Q Max ## -26.4193 -5.0445 -0.0445 5.5807 23.5807 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 193.7941 0.8433 229.8 <2e-16 *** ## SNPSEX -14.3748 0.5571 -25.8 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 7.719 on 781 degrees of freedom ## Multiple R-squared: 0.4602, Adjusted R-squared: 0.4595 ## F-statistic: 665.7 on 1 and 781 DF, p-value: < 2.2e-16
Неудивительно, что биологический пол является хорошим предиктором роста, улавливая ~ 46% дисперсии - женщины ниже мужчин. Теперь мы делаем прогнозы.
Sex3<-as.data.frame(setdiff(Sexcheck$IID,Sex2$IID))
colnames(Sex3)<-c("IID")
Sex3<-merge(Sex3,Sexcheck, by = "IID")
Sex3<-filter(Sex3, SNPSEX > 0)
Sex3$pred2<-predict(fit1,Sex3)
write.table(Sex3$pred2,file="PredSex2.txt",sep = ",",col.names = F,row.names = F)Вот как складывается это предсказание о том, что рост каждого человека является средним для его пола. В независимом наборе тестов только пол предсказывает 42,08% дисперсии.
3) Удалить связанных людей
Итак, следующим шагом будет идентификация членов семьи и их удаление. Я использовал флаг plink's –genome, который вычисляет оценку идентичности по происхождению. IBD отражает родство двух людей:
1.00 = Тот же человек или монозиготные близнецы
0,50 = родитель / ребенок или братья и сестры
0,25 = двоюродные братья и сестры
и т.д…
Пороговое значение 0,1875 (на полпути между двоюродными братьями и сестрами) довольно типично.
Итак, я извлечу все значения IBD выше нижнего порога 0,05 и нанесу на график все значения, превышающие IBD = 0,1.
system(paste ("plink --bfile OpenSNP_combined5 --genome --min 0.05 --out OpenSNP_combined5_IBD",sep=" "))
Мы видим несколько пар участников, которые сильно связаны, а затем одну ОГРОМНУЮ группу людей, которые все связаны друг с другом. На самом деле все эти люди являются троюродными братьями, но эта расширенная семья намного больше, чем можно было бы ожидать случайно. Имейте в виду, что это все люди, которые добровольно загрузили данные о своих генотипах на общедоступный веб-сайт и сделали его полностью открытым. Люди делают это по многим причинам, хотя одна из них: если у них есть какая-то болезнь, которая присутствует в семье, они могут надеяться узнать о ней больше.
Я удалил по одному участнику из каждой пары с IBD ›0,1875 и всех, кроме одного, из огромной группы.
Train_Exclude<-c(567, 191, 849, 796, 569, 680, 461, 706, 688, 457, 746, 737, 263, 390, 61, 529, 55, 338, 492, 846, 392, 511, 901, 457)
Train_Exclude<-cbind(Train_Exclude,Train_Exclude)
write.table(x=Train_Exclude, file="TrainExclude.txt", sep=" ", quote = F, row.names = F, col.names = F)
system(paste ("plink --bfile OpenSNP_combined3 --remove Train_exclude.txt --make-bed --out OpenSNP_combined7",sep=" "))4) Генерация основных компонентов
Мы можем использовать plink для генерации 20 основных компонентов.
system(paste ("plink --bfile OpenSNP_combined7 --pca 20 --out PCs20_SNP7",sep=" "))
Мы видим, что у первого ПК очень большое собственное значение. Это признак расслоения населения.
Оранжевый - тренировочный, синий - тестовый. Это каноническая форма полиэтнической выборки (см. Здесь рисунок 5). В этом комке около начала координат находятся все белые люди, а хвосты имеют африканское и азиатское происхождение (NB, это каноническая форма образца, состоящего в основном из белых).
Поскольку и обучающая, и тестовая выборки являются многоэтническими, на этом этапе я никого не удалял. Оглядываясь назад, может быть, мне следовало. Если бы я удалил всех небелых людей из обучающих данных, это, вероятно, улучшило бы мой прогноз в большей части тестовой выборки, что улучшило бы общую производительность моей модели.
Но в моем фактическом представлении я просто поправил три основных компонента. Для исследования этого, вероятно, будет недостаточно. Вот как это выглядит:
Sexcheck<-fread("OpenSNP_combined6.sexcheck") Sex2<-merge(Sexcheck,Height, by="IID") Height<-fread('phenotypes.txt')Sex2<-merge(Sexcheck,Height, by="IID") Dat<-merge(PCs, Sexcheck, by = "IID") Test<-fread('OpenSNP_test.fam') colnames(Test)[1]<-"IID" Test<-merge(Test,Dat, by="IID")Train<-merge(Height, Dat, by="IID")Fit_all<-lm(as.formula(paste("Height ~ SNPSEX + ", paste("PC",c(1:3),sep="",collapse = "+") )),data=Train) summary(Fit_all)## ## Call: ## lm(formula = as.formula(paste("Height ~ SNPSEX + ", paste("PC", ## c(1:3), sep = "", collapse = "+"))), data = Train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -25.4227 -5.2219 -0.1535 5.0593 22.3754 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 193.9242 0.8474 228.838 < 2e-16 *** ## SNPSEX -14.4399 0.5630 -25.648 < 2e-16 *** ## PC1 20.3413 9.0280 2.253 0.024537 * ## PC2 -28.0329 8.3567 -3.355 0.000835 *** ## PC3 18.6254 8.1439 2.287 0.022470 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 7.651 on 756 degrees of freedom ## Multiple R-squared: 0.4672, Adjusted R-squared: 0.4644 ## F-statistic: 165.7 on 4 and 756 DF, p-value: < 2.2e-16
Основные компоненты улучшают производительность нашей модели в обучающих данных примерно на 0,5%.
Точно так же основные компоненты улучшают производительность в тестовых данных на ~ 1%. Теперь модель объясняет 43,28%.
5) Создайте рейтинг полигенного риска
У доктора Сары Медланд есть отличная онлайн-презентация о том, как создавать шкалы полигенного риска. Я не буду вдаваться в подробности. Я использую PRSice для создания своих PRS. На самом деле, я использую собственный R-скрипт, основанный на PRSice, потому что PRSice предназначен только для Linux и не дает результатов с нужными мне порогами, но результат тот же. Я сгенерировал PRS для высоты с несколькими разными пороговыми значениями и нанес на график величину дополнительной дисперсии, которую они объясняют в обучающей выборке ниже.
Мы видим, что PRS объясняет дополнительные 6,9% - 10,0% дисперсии, при этом порог p очень изучены). Низкая производительность, вероятно, связана с многонациональной выборкой.
Вот как эта модель сравнивается в тестовом наборе:
PRS объясняет дополнительные 10% дисперсии тестовых данных. Приятно видеть, что в тесте он работает примерно так же хорошо, как и в обучающих данных!
Вот как я это сделал. Никакого машинного обучения, только экспертные знания. Как я уже сказал выше, я делюсь этими результатами в надежде, что другие смогут их улучшить. Если да, пожалуйста, крикни мне !!
