Генотип, ко-факторы, исход, как анализировать? |
Здравствуйте, гость ( Вход | Регистрация )
Генотип, ко-факторы, исход, как анализировать? |
15.06.2014 - 14:09
Сообщение
#1
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
Здравствуйте, коллеги! Прошу помощи в анализе данных.
Задача исследования - оценить связь между генотипом (15 SNP), "промежуточным фенотипом" (параметры биохимии, иммунологии и др.), исходом (ЗНО есть/нет). Существующий сервис "SNPstats" (http://bioinfo.iconcologia.net/SNPstats_web) выдает отношения шансов, "adjustet by фактор1+фактор2+...", используется при этом "logistic regression models" (то есть, логит-регрессию?). Хотелось бы поточнее узнать, что значит "adjusted by". Кроме того, в данной програме остается "за кадром", какой из факторов является ведущим. Возможно, есть какие то альтернативные методы анализа, позволяющие оценить вклад конкретных факторов? Посоветуйте, пожалуйста. Заранее благодарен. И есть ещё один вопрос: По разным SNP имеется разное количество генотипированных, как и разное количество известных значений по каждому из "промежуточных фенотипов" и исходов. То есть, грубо говоря, выборки по каждому из SNP перекрываются только отчасти. Нужно ли в этом случае рассматривать проблему множественных сравнений? Спасибо! |
|
15.06.2014 - 15:15
Сообщение
#2
|
||
Группа: Пользователи Сообщений: 116 Регистрация: 20.02.2011 Пользователь №: 23251 |
1) Что означает adjusted..Вы включаете в модель один фактор_1, рассчитываете отношения шансов для этого одного "снипа"..затем создаете новую модель, в которую включаете фактор_2, получаете отношения шансов для фактора_2. Теперь, если Вы включите фактор_1 и фактор_2 в одну модель, то получите скорректированные (adjusted) отношения шансов для фактора_1 и фактора_2, поскольку в модель включаются оба фактора, поэтому OR фактора_1 корректируются относительно присутствия фактора_2 и наоборот. То есть adjusted в данном случае - это просто расчет OR относительно количества предикторов, включенных в модель..
2) Как Вы знаете, снипы (SNPs) сами по себе обладают очень маленькими величинами эффектов, то есть сам по себе 1 полиморфизм не может определить предрасположенность к заболеванию (если Вы только не исследуете моногенные заболевания, в чем я сомневаюсь). Следовательно, риск развития наиболее распространенных полигенных заболеваний, таких как артериальная гипертензия, инфаркт миокарда, сахарный диабет и др., определяется совокупностью влияний снипов, а также факторами внешней среды. Поэтому, изучая вклад каждого снипа без учета взаимодействия (строя аддитивные модели) Вы никогда не докопаетесь до истины. Для этого необходимо изучать межгенные взаимодействия с целью определить, например, усиливает ли Снип1 влияние Снипа 2, какая связь между Снипом1 и Снипом2? Для решения подобных задач был разработан целый спектр подходов и алгоритмов, основные из них приведены на рисунке ниже. Из этого списка очень хорошие результаты показывают Multifactorial Dimensionality reduction и Random Forest с модификацией Catherine Strobl. Про первых подход можно ознакомиться на сайте MDR: http://www.multifactordimensionalityreduction.org/ Статьи про MDR, например: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1226028/ Вторую прикрепил здесь. multifactor1.pdf ( 648,61 килобайт ) Кол-во скачиваний: 1016 Метод хороший и работает как для подхода ген-кандидат (как в Вашем случае), так и для полногеномных исследований, когда количество снипов будет от 100 тысяч (здесь логистическая регрессия загнется в результате "curse of dimensionality"). 3) Какой из снипов делает наибольший вклад в риск развития заболевания позволит Вам определить, например, Random Forest с модификацией Catherine Strobl ("тетя" из Германии, которая всю жизнь занимается лишь случайными лесами, ей можно доверять ). Погуглите Strobl C. An Introduction to Recursive Partitioning: Rationale, Application and Characteristics of Classification and Regression Trees, Bagging and Random Forests. Должны быть статьи в свободном доступе, если не найдете, вышлю Вам. Также этот метод реализовал в скрипте R, если работаете в R (или пока не работаете) - тоже могу выслать. Для статьи в хорошем журнале или для кандидатской (докторской) хорошего уровня, расчет лишь OR c р-values - явно недостаточно. Любая хорошая конференция, посвященная проблемам биоинформатики и определения SNP's, включает вышеприведенные методы. Сообщение отредактировал TheThing - 15.06.2014 - 15:26 |
|
|
16.06.2014 - 00:05
Сообщение
#3
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Если данных достаточное количество, то стоило бы попробовать построить Bayesian networks http://www.bnlearn.com/
|
|
17.06.2014 - 09:09
Сообщение
#4
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
Благодарю, коллеги, за ответы. Буду разбираться...
|
|
25.09.2014 - 09:47
Сообщение
#5
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
Добрый день, коллеги!
Изучение вопроса привело меня к необходимости использовать random forests, а именно - модификацию, доступную в пакете party для R (cforest), где для построения леса используются "безошибочные" деревья, т.к. в этом случае адекватно оценивается важность переменных и нет преференций для количественных данных (как в "стандартном" random forests)[Strobl, 2007]. Кроме того, необходимо вычислять "кондиционную" значимость, которая адекватнее в случае коррелированных переменных [Strobl, 2009]. И тут возникает следующее препятствие: missing values, которые в моих данных для некоторых параметров составляют более 50%. Однако, и на этот счет у той же "тёти" Strobl нашлось решение [Hapfelmeier, 2013]. На данном этапе моего "мастерства" в R хватает только на импорт таблицы, вычисления значимости переменных и построения гистограммы значимости: read.table("gen2.txt", h=T)->gen2 set.seed(100) cfresults<-cforest(cr~.,data=gen2,controls=cforest_unbiased(mtry=10,ntree=100)) varimp(cfresults)->vi1 par(las=2) barplot(sort(vi1),main="Variable importance",space=0.5,cex.names=0.6) Прошу помощи у бывалых: 1) Каким должен быть код для вычисления кондиционной важности переменных при наличии пропусков данных? В приложении к [Hapfelmeier, 2013] есть код, но там другие задачи исследования и я не могу разобрать что к чему. 2) Как нарисовать индивидуальные деревья в cforest и правильно интерпретировать их? Может быть, кто - то соизволит поковыряться с моими данными, на этот случай прилагаю таблицу. Также прилагается Hapfelmeier, 2013 Сообщение отредактировал don - 25.09.2014 - 09:51
Прикрепленные файлы
Hapfelmeier_2013.pdf ( 368,37 килобайт )
Кол-во скачиваний: 387
gen2.txt ( 32,9 килобайт ) Кол-во скачиваний: 493 |
|
28.09.2014 - 19:41
Сообщение
#6
|
|
Группа: Пользователи Сообщений: 219 Регистрация: 4.06.2013 Из: Тверь Пользователь №: 24927 |
Попробуйте прочитать о важности переменных здесь http://gewissta.livejournal.com/ / Метод случайного леса.../
Код для вычисления обычной важности переменных там точно есть. Я немного работал с платформой Salford Predictive Modeler http://lib.giveawayoftheday.com/b/predictive_analytics/ Стандартные вычисления в среде этой платформы очень легко выполнять, просто используя интерфейс. Но и с кодами также можно работать. Сообщение отредактировал anserovtv - 28.09.2014 - 19:44 |
|
28.09.2014 - 19:44
Сообщение
#7
|
|
Группа: Пользователи Сообщений: 219 Регистрация: 4.06.2013 Из: Тверь Пользователь №: 24927 |
Коды, генерируемые программой, можно преобразовать в коды R или SAS и в другие форматы.
Сообщение отредактировал anserovtv - 28.09.2014 - 21:05 |
|
29.09.2014 - 14:00
Сообщение
#8
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
anserovtv
Большое Вам спасибо за первую ссылочку (до второй пока не добрался, рабочая сеть не позволяет)! Наглядно и доступно представлен алгоритм "классического" random forests. Однако, автор верно заметил, что классический вариант RF "не переваривает" анализ смешанных количественных и категориальных данных, в этом случае необходим cforest. Коллеги, если кто-то натыкался на пример использования cforest, будьте добры, поделитесь! Сообщение отредактировал don - 29.09.2014 - 14:02 |
|
29.09.2014 - 15:07
Сообщение
#9
|
|
Группа: Пользователи Сообщений: 116 Регистрация: 20.02.2011 Пользователь №: 23251 |
anserovtv Большое Вам спасибо за первую ссылочку (до второй пока не добрался, рабочая сеть не позволяет)! Наглядно и доступно представлен алгоритм "классического" random forests. Однако, автор верно заметил, что классический вариант RF "не переваривает" анализ смешанных количественных и категориальных данных, в этом случае необходим cforest. Коллеги, если кто-то натыкался на пример использования cforest, будьте добры, поделитесь! C кодом cforest ответил Вам в личку. Не совсем понятно, зачем Вам рассматривать каждое дерево по отдельности и что-то интерпретировать, ведь каждое дерево имеет свою индивидуальную структуру и свой прогноз, а конечный результат сводится к голосованию за результат, который дает большинство деревьев. Можно вывести структуру каждого из 100 (как в Вашем случае) деревьев, но потом все это интерпретировать будет заданием очень непростым и вряд ли полезным. После того, как cforest выдает Вам наиболее значимые предикторы на уровне СНИПов, Вам интересно ведь посмотреть какой генотип ассоциируется с риском развития заболевания или имеет наоборот протективный эффект. Для этого очень удобно использовать комбинацию методов: случайный лес, MDR, логистическая регрессия. Случайный лес - в качестве фильтра наиболее важных СНИПов, MDR - оценка взаимодействия СНИПов и оценка каждого из генотипов, а также их комбинаций, нахождения типа связи между СНИПами. Лог. регрессия - финальный этап, создание модели, которая включает лишь факторы, которые определились как важные из первых двух методов, валидация результатов. P.S. В Вашей базе - что прогнозируется? nat или cr? Если cr - получается, что исследуете очень редкое заболевание? Сообщение отредактировал TheThing - 29.09.2014 - 15:12 |
|
29.09.2014 - 21:09
Сообщение
#10
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
P.S. В Вашей базе - что прогнозируется? nat или cr? Если cr - получается, что исследуете очень редкое заболевание? Прогнозируется cr (некое заболевание). Это урезанная база, здесь cr 24 из 165, в оригинальной базе 55 из 439. Но в оригинальной базе 45 предикторов, включая 13 снипов, и громадное количество пропусков в некоторых из них (более 50%). Поэтому пришлось несколько сузить масштабы анализа, убрав часть предикторов. Возможно, это неверно. Поправьте меня, если я не прав. |
|
29.09.2014 - 21:41
Сообщение
#11
|
|
Группа: Пользователи Сообщений: 116 Регистрация: 20.02.2011 Пользователь №: 23251 |
Прогнозируется cr (некое заболевание). Это урезанная база, здесь cr 24 из 165, в оригинальной базе 55 из 439. Но в оригинальной базе 45 предикторов, включая 13 снипов, и громадное количество пропусков в некоторых из них (более 50%). Поэтому пришлось несколько сузить масштабы анализа, убрав часть предикторов. Возможно, это неверно. Поправьте меня, если я не прав. Мне сложно Вас поправить, поскольку я практически не знаю дизайна исследования, целей, задач и т.д. - если Вы решили урезать массив, значит на то были причины. Посмотрел немного Вашу базу, прикольная (в том смысле, что собрано столько материала!!). Однако, например, первые три СНИПа IL1b_rs1143634, IL2_rs2069762, IL4_rs2070874, а точнее распределение их генотипов в группах больных и здоровых (cr) не позволяет их нормально анализировать, поскольку наблюдается или perfect separation (то есть классы больных и здоровых идеально расходятся и большинство алгоритмов сходит с ума поэтому) или алгоритм просто не сходится (did not converge), поскольку присутствуют не все комбинации генотипов в базе данных (может быть потому что Вы ее сократили) - это беда всех параметрических методов, в том числе и логит-регрессии(curse of dimensionality). Что-то нужно думать с базой.. |
|
30.09.2014 - 06:55
Сообщение
#12
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
...если Вы решили урезать массив, значит на то были причины. Представьте следующую картину. У нас есть большая база, порядка полутора тысяч чел., в которой собраны наблюдения по разным параметрам, но с разной степенью "перекрывания" данных между ними. Нам известны: 1) генотипы по нескольким снипам для 400 чел., 2) параметр Х1 для 200 чел, из которых со снипами "перекрывается" 150, 3) параметр Х2 для 300 чел, из которых со снипами перекрывается 250, а с Х1 - только 15 4) параметр Х3 для 100 чел, из которых все перекрываются со снипами, но с Х2 и Х1 только 10-20. 5) и таких параметров еще десяток, с разной степенью "перекрывания" с другими. 6) для каждого из пациентов известен статус по болезни. Чтобы лучше понять, о чем идет речь, приложу файл с большой первоначальной базой (don_base111.txt). Причина урезания базы заключается в попытке отобрать тех пациентов, для которых известно как можно больше параметров. Делается это из предположения, что для адекватного анализа RF пропущенные значения конечно могут быть, но должны быть сведены к какому то разумному пределу (не более 30 %). То есть, до урезания базу можно представить так (база1.png), после - как-то так (база2.png). Сообщение отредактировал don - 30.09.2014 - 07:19
Прикрепленные файлы
|
|
30.09.2014 - 08:54
Сообщение
#13
|
|
Группа: Пользователи Сообщений: 116 Регистрация: 20.02.2011 Пользователь №: 23251 |
Представьте следующую картину. У нас есть большая база, порядка полутора тысяч чел., в которой собраны наблюдения по разным параметрам, но с разной степенью "перекрывания" данных между ними. Нам известны: 1) генотипы по нескольким снипам для 400 чел., 2) параметр Х1 для 200 чел, из которых со снипами "перекрывается" 150, 3) параметр Х2 для 300 чел, из которых со снипами перекрывается 250, а с Х1 - только 15 4) параметр Х3 для 100 чел, из которых все перекрываются со снипами, но с Х2 и Х1 только 10-20. 5) и таких параметров еще десяток, с разной степенью "перекрывания" с другими. 6) для каждого из пациентов известен статус по болезни. Чтобы лучше понять, о чем идет речь, приложу файл с большой первоначальной базой (don_base111.txt). Причина урезания базы заключается в попытке отобрать тех пациентов, для которых известно как можно больше параметров. Делается это из предположения, что для адекватного анализа RF пропущенные значения конечно могут быть, но должны быть сведены к какому то разумному пределу (не более 30 %). То есть, до урезания базу можно представить так (база1.png), после - как-то так (база2.png). Наверное собирали базу одну люди, а анализировать приходится Вам Может быть для начала изучить полностью вопрос с генетикой по максимальной базе без добавления доп. клин. параметров? Если что-то найдем по снипам, будем хотя-бы знать, что некоторые ассоциированы с заболеванием. Если снипы обнаружить не получиться, это значительно облегчит проблему перекрытия в базе, поскольку тогда стоит копаться лишь в клин. параметрах. Пока гляну сегодня по урезанной базе, что с остальными снипами.. |
|
30.09.2014 - 10:34
Сообщение
#14
|
|
Группа: Пользователи Сообщений: 24 Регистрация: 11.06.2014 Пользователь №: 26460 |
Наверное собирали базу одну люди, а анализировать приходится Вам Да, я пытаюсь объединить в одном анализе сразу несколько показателей из разных баз. Может быть для начала изучить полностью вопрос с генетикой по максимальной базе без добавления доп. клин. параметров? Если что-то найдем по снипам, будем хотя-бы знать, что некоторые ассоциированы с заболеванием. Я пробовал с помощью SNPStats (http://bioinfo.iconcologia.net/SNPstats) анализировать только снипы, по одиночке. Алгоритм заключается в анализе Хиквадратом сводной таблицы (если без кофакторов) или лог-регрессией (с кофакторами). Там получилась ассоциация со снипами в TNFa и IL4. Если с поправкой Бонферрони, то только TNFa. Если снипы обнаружить не получиться, это значительно облегчит проблему перекрытия в базе, поскольку тогда стоит копаться лишь в клин. параметрах. А если имеет место взаимодействие снипов с параметрами? То есть, сам по себе снип и среди других снипов не информативен, а среди тех, у кого повышен какой то параметр, он играет роль? Сообщение отредактировал don - 30.09.2014 - 10:38 |
|
30.09.2014 - 14:07
Сообщение
#15
|
||||
Группа: Пользователи Сообщений: 116 Регистрация: 20.02.2011 Пользователь №: 23251 |
А если имеет место взаимодействие снипов с параметрами? То есть, сам по себе снип и среди других снипов не информативен, а среди тех, у кого повышен какой то параметр, он играет роль? Ну может и в той части базы данных, которую Вы урезали содержится как раз информация про взаимодействие, но мы этого никогда не узнаем (пока ) . Поэтому мне кажется стоит с чего-то начать, например с "голой" генетики.. Вот что получилось по урезанной базе по всем снипам при прмощи метода MDR: Классификационная способность на выборке для обучения: 78% Классификационная способность на выборке для тестирования: 63,8% (ну уже что-то ) Наилучшая модель - 3-х локусная: IL6_G(-174)C,IL8_rs4073,IL12_rs3212227 кросс-валидация 10 из 10 Правила для комбинаций генотипов: IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 2 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 1 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 1 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 1. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 2 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 0 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 1 THEN CLASSIFY AS 0. IF IL6_G(-174)C = 0 AND IL8_rs4073 = 2 AND IL12_rs3212227 = 0 THEN CLASSIFY AS 1. 0 - AA 1 - Aa 2 - aa Classify as 0 - Control Classify as 1 - Case IL6_G(-174)C,IL8_rs4073,IL12_rs3212227 Cross-validation Statistics: Training Balanced accuracy: 0.789 Training Accuracy: 0.6788 Training Sensitivity: 0.9444 Training Specificity: 0.6336 Training Odds Ratio: 29.3935 (4.5023,191.8959) Training Χ²: 24.9117 (p < 0.0001) Training Precision: 0.3049 Training Kappa: 0.3091 Training F-Measure: 0.461 Testing Balanced accuracy: 0.6348 Testing Accuracy: 0.6121 Testing Sensitivity: 0.6667 Testing Specificity: 0.6028 Testing Odds Ratio: 3.0357 (0.169,54.5173) Testing Χ²: 0.6057 (p = 0.4364) Testing Precision: 0.2222 Testing Kappa: 0.1473 Testing F-Measure: 0.3333 Cross-validation Consistency: 10/10 |
|||
|