Воркшоп по биоинформатике

[progenes]

Я предполагаю, что меня читают биологи и информатики. Поскольку я сейчас столкнулась с феерическими расчетами, от которых у меня волосы на загривке вздыбились, считаю, что будет неплохо, если я расскажу причастным где впредь быть предельно внимательным. Я, к сожалению, не могу дотянуться до того парня, который это посчитал, чтобы надавать ленейкой по рукам. Но должна сказать, что это сервис, претендующий на серьезный.

Больше чем месяц назад я получила результаты и, помнится, даже всхипнула от ужаса. Дело было вот как.

Перед биологом стоит задача - сравнить экспрессию генов в нескольких тканях. Для этого надо выделить РНК и каким-то из методов определить, с каких генов эта РНК считалась и в каком количестве. Методов есть несколько. От нозерна до микрочипов. Самый модный - это секвенирование 454. О нем и пойдет речь.

Я не буду вдаваться в подробности, что за ткани я анализирую, скажу только, что этой ткани столько, что невооруженным глазом не видно, нарубили лазером под микроскопом. Выделили РНК в количествах на пределе воображения и отправили на прочитку (секвенирование). Это приблизительно так, как сейчас читают геномы, только не ДНК, а РНК. Причем прочитка - это полноценный сервис, который включает все. На выходе, как я уже упоминала 300 гигов информации: сырые сиквенсы, сбивка в контиги, скаффолды и унигены, бласты, функциональная аннотация по геномной антологии, метаболитических путях, визуализация, статистическая обработка и дифференциальный анализ. Все растыкано по 2000 тыщам файлов. Задача биолога теперь все ОСМЫСЛИТЬ и интерпретировать и сделать выводы, как же отличается работа генов в разных тканях и почему.

Я сузила сначала задачу и из вороха файлов нашла исходник в экселе: сравнение генной экспрессии в двух (из 48ми) тканях. Теперь внимание и пристегнитесь. Результаты сравнения представлены в виде log2 значения соотношения экспрессии генов в ткани 1 vs. 2. Просто настолько, что можно понять и идиоту, верно? Значения колеблются от -15 до плюс 15 (это уже log2). Огого, сечете разницу в экспрессии? Всех генов несколько десятков тыщ.

В этом месте меня подвело банальное любопытство. Что ж это за ген такой, у которого разница в работе в двух близлежащих тканях 2¹⁴? Роюсь в ворохе файлов и нахожу сырые результаты. И тут, друзья, у меня глаза на переносице и сбежались. Потому что я знаю, что это за сырые значения. А это всего навсего количественный подчет ШТУК КУСКОВ РНК, которые принадлежат одному гену. Держитесь теперь крепче, пример из жизни.

Ген Х. Логарифмированное значение соотношения экспрессии в тканях 1 vs. 2 равно 14,72. Сырые данные в студию. В ткани Nr.1 насчитали 0 (ноль) кусков, в ткани Nr.2 насчитали 27 кусков. Формула расчета log2(27:0). Что, съели касатики?!!! Говорите на ноль нельзя делить? Ну нельзя так нельзя (хотя в результатах стоит ноль). Я прикинула, как могли бы рассуждать те, кто уныло смотрит на ноль. Ноль надо заменить на число, отличное от нуля. Я начала тупо подставлять и проверять логарифмом, как у них 14,72 получилось. Оказалось, что 0 приравняли до 0,001. log2(27:0,001)=14,72

Если вы еще не ржете, посчитайте, какое значение log2 будет, если в ткани Nr.1 насчитали 0 (ноль) кусков, в ткани Nr.2 насчитали 2 (два) куска. То есть вы догадываетесь, куда можно засунуть эти 300 гигов и ограничиться одной таблицей в экселе, да?

Проблема в том, что редкий биолог интересуется сырыми данными, если сервис поставил уже готовое соотношение в красочных схемах и диаграммах. И редкий информатик интересуется особенностями того, что ему поручено посчитать. Для него это голые абстрактные числа, а для меня это штуки кусков РНК. Особенно печально, если биолог мало знает о проблемах и недостатках того метода, которым хочет что-то проверить. После таких ляпов у меня возникли подозрения к алгоритмам сбивки в контиги (которые я встречала в других случаях), к проблемам аннотации (с которой долбилась годами). Эта же проблема также касается и анализов всяких там аффиметриксов и прочих биочипов.

Вырасту большой и научусь программированию сама.

Comments

[psilogic.livejournal.com]

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

[progenes.livejournal.com]

Через мои руки прошло уже столько информатиков, что у меня есть ощущение, что проще самой научиться программированию, чем в тысячный раз объяснять программисту, что я от него хочу. Язык болит, чесное слово. Все одни и те же грабли. В этот раз была надежда, что сервис прошел стадию диалогов с генетиками. Ан нет, надо все еще перепроверять.

[stoshagownozad.livejournal.com]

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

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

[progenes.livejournal.com]

Математика простая. В одной клетке сидит 27 кроликов. Другая клетка пустая. Вопрос - во сколько раз в одной клетке больше кроликов, чем в другой. Ответ: в 27 тыщ раз.

[vigna.livejournal.com]

А как анализировали, это ваши программисты какие-то свои собственные программы писали или что-то готовое, типа той же CLC genomics?

[progenes.livejournal.com]

Нет, это не наши программисты. Это готовое.

[tapka.livejournal.com]

А вот расскажите мне, как почти информатику, какое бы вы хотели видеть число в этой ситуации в той эксельной таблице? inf? Вообще наверное помимо logFC нужно там добавлять какое-нибудь p-value. Типа 20 кусков и 10 это хорошая разница, а 2 и 0 ничо не значит. Ну и порог какой-нибудь сколько должно быть кусков вообще, чтоб мы сказали, что ген активно экспрессируется.

[progenes.livejournal.com]

p-value есть. в этом примере он равен 5,84E-07 (что мало мне говорит) Его считали по формуле , где N1 - total clean tag number of sample 1 и N2 total clean tag number of sample 2.

А вот где проходит трешхолд абсолютных значений, p-value и чувствительности метода - я сказать не могу, пробую сама разобраться. Наверное, это можно сказать только глядя весь массив данных и зная особенности технических ошибок.

[mopexod.livejournal.com]

Ну, это общая проблема, не только у биологов. Очень мало кто умеет правильно разобраться в физической сущности измерений. Да что там разобраться, даже правильно померять мало кто может!

[malck.livejournal.com]

Нуачо, с точки зрения вычислений это очень даже нормально. Скажите "спасибо", что 0,001 (= 1E-3), а не 1E-10, например. :) А то результаты были б еще фееричнее. :)

[hoegni.livejournal.com]

Лучше бы как раз взяли 10^(-10). Тогда бы результат был очеыидно неправдоподобен, и Русланка бы его выкинула сразу без сомнений.

[oude-rus.livejournal.com]

такой фокус - это довольно типичный прием, когда не хотят делать отдельное поле "результат недействителен". Причем число-с-потолка берется такое, чтобы по результату было сразу видно: ерунда.

Посмотрите в описании пакета, это наверняка где-то оговаривается. Другое дело, что реальные пацаны описаний не читают.

[tapka.livejournal.com]

Что значит ерунда? С нулем может быть как интересный, так и не интересный вариант. Если в одной ткани 1000, а в другой 0, то это например интересно, а если 5 и 0, нет. Отправлять в отстой все результаты, где в одном случае 0, это значит лишиться возможно чего-то стоящего. Но и выискивать из кучи ерунды интересное тоже то еще занятие.

Вообще это ж вроде начала статистики. Все эти бросания монеток, разноцветные шары... Должна быть какая-то разумная метода оценки достоверности различий для случаев целочисленных измерений, где ноль значащее число, без приписок. Просто сервисцентр использует какой-то не тот критерий.

[hoegni.livejournal.com]

Ээ... Но ведь ты поняла, что результат - неправдоподобен? Таким образом, ситуация обработана более-менее правильно. Там, где у тебя вменяемое количество РНК на входе, там и вменяемый результат на выходе.

А вообще да, очень малые значения в одной из групп - это всегда источник головной боли при расчете относительных величин и нормализации, что поделать.

[progenes.livejournal.com]

Ну я то поняла, я ж не только им, я и себе не доверяю и перепроверяю. Но я видела буквально недавно. Товарищ полез в Genvestigtor и натаскал оттуда экспрессии по интересующему его гену. А там красиво, визуализация. На голубом глазу втирает - вот в корнях ген в сто раз сильнее экспрессируется чем в листках, это что-то ж значит! А с боку на картинках бар в абсолютных значениях. Говорю, - родной, ты ж видел, что там 0,1 и 10? А аффиметриксовый трешхолд на абсолютную экспрессию я ставлю 100, чтобы не нагрести шумов. А ты про 0,1 и 10. Отаке.

[ixaydan-xti.livejournal.com]

Можно попросить информатика пообъяснть.. что он делал...
И попросить переделать как вам надо :)))
Вот и всё :)))

Просто не успеешь всё выучить...

[progenes.livejournal.com]

Увы, это как раз не тот случай, что я могу попросить объяснить и попросить переделать. Я могу попробовать разобраться, как он делал и, если возможно, переделать самой.

[potar1.livejournal.com]

ржала :):)
и как выглядела рецензия? цензурных слов хватило?

[progenes.livejournal.com]

Ну я в процессе. Уже постепенно прозреваю, как все это по новой пересчитать.

ехидненько так

[malck.livejournal.com]

А вот интересно... а какова в этих результатах вероятность ошибок первого и второго рода? :D

Re: ехидненько так

[progenes.livejournal.com]

очень высокая. то есть в том виде, в каком у меня эти результаты, серьезной статистики проводить невозможно. скажу честно - эта проблема с нулями самая маленькая в этом эксперименте.

[krydelta.livejournal.com]

Хм, може просто змістити кільскість генів на одиницю?
0 генів трактується 1 ген, 1 ген - як два гени. І таким чином: log2(10/0) ~ log2(11/1). Це можна робити як для всіх пар, а можна тільки для тих, де є 0.

Це трошки спотворить результат в загальному, але значно полегшить роботу з нулями.

[progenes.livejournal.com]

Як варіант.

[paparazzzzi.livejournal.com]

Не будучи биологом, ничего в генах не поняла, но от методики офигела.

[tiy.livejournal.com]

На мой взгляд, это совершенно стандартная проблема с нулевыми ячейками. Если считать p-value в любом исследовании банальным методом chi-square test, то при наличии пустых ячеек к ним добавляют какую-нибудь цифирь поменьше (и называется это, если я не ошибаюсь, Haldane correction).

Разница прочиток (reads, как это будет по-русски?!) между 0 и 1 Вас, конечно, не интересует, а вот 0 и 27 уже, возможно, указывает на тканевую специфичность (зависит от покрытия (coverage), трудно сказать 27 это много или мало, если не знать как обстоят дела в среднем по госпиталю). Так что совсем выкидывать эти результаты я бы не стала. Можно установить границу на минимальную экспрессию в одной из тканей, и далее выкидывать все, что не соответсвует критериям.

В качестве способа иллюстрации уровня экспрессии наши товарищи используют в картинках дополнительную вертикальную ось, где отмечают уровень экспрессии определенной пробы по сравнению со всеми другими тестируемыми пробами - если открыть вот эту страничку - http://www.nature.com/ng/journal/v42/n4/extref/ng.543-S1.pdf и поискать в тексте supplementary figure 2, то под шапкой к рисунку будет много картинок, где справа на оси Y есть "риска" для каждой пробы. Чем она выше, тем выше экспрессия этой пробы по отношению ко всем остальным. В данном случае использовалось не секвенирование а genome-wide illumina gene expression chip, но можно придумать что-то подобное и в Вашем случае секвенирования.

Год назад в Nature была пара статей про RNA-секвенирование в качестве оценки экспрессии, возможно, там найдется полезная информация на тему того, какой они использовали cut-off для анализа экспрессии:
http://www.nature.com/nature/journal/v464/n7289/pdf/nature08872.pdf
http://www.nature.com/nature/journal/v464/n7289/pdf/nature08903.pdf

(статья Притчарда несколько более разумная, я бы с нее начинала). Впрочем, сейчас довольно много печатают на эту тему, из того что я помню - неплохое ревью на тему RNA секвенирования было в январском номере Nat.Rev.Genetics.

К вопросу об анализе данных - в Гронингене на кафедре генетики есть несколько открытых позиций (PhD, postdoc) в проектах связанных с экспрессией генов. Так вот на них с гораздо большей радостью возьмут математиков, чем биологов. Все-таки биологии легче научиться по ходу дела, чем математике и программированию (хехе, это мой хронический комплекс - почему я не программист...).

Успехов! :)

[progenes.livejournal.com]

tiy (Sasha Zhernakova) не понимаю, куда подевался Ваш комментарий, но на почту он мне пришел. Спасибо!

[tiy.livejournal.com]

вижу его прямо над Вашим комментарием :)

[shao-s.livejournal.com]

Эх, Руслана, читаю и плачу... Примерно после таких закидонов я изучил программу Wolfram Mathematica. И теперь, что характерно, периодически делаю матмоделирование для других. А они не проверяют, а верят так...

[progenes.livejournal.com]

Издеваешься, да?

[ad3002.livejournal.com]

Это еще что, я пытаюсь сравнить несколько ChIP-Seq - вот это да.

Даже RNA-Seq математика весьма сложная. И если использовать терминологию кроликов и клеток, то вы ловили зверей в одной очень большой клетке набитой зверями, и в другой. В первой поймали много зверей и среди них 27 с ушами как у кроликов. Во второй тоже много зверей, но среди них с ушами как у кроликов не было.

При таком раскладе ответить на вопрос "Насколько кроликов больше во второй клетке, чем в первой?" трудновато. Да даже корректно оценить вероятность того, что в первой их больше не очень просто.

[snysmymrik.livejournal.com]

Ну, качественные результаты же получены? В ткани А гена Х сильно больше чем в ткани Б. Чего ещё надо биологу, чтобы встретить старость?

[toi-samyi.livejournal.com]

нельзя никому верить
самому надо проверять

[sha-lun.livejournal.com]

Кстати да, попробуйте Wolfram Mathematica.
Из бесплатного есть языки программирования Python (python.org) и J (jsoftware.com)
Не с паскаля же сегодня начинать обучение.

[no1u1w1w6c]

а ещё есть Scheme, у которой в numeric tower изначально чёткое разделение на обычные «float», дроби и bignum'ы. и прозрачная работа с этим всем. и учится за два вечера с пивом. но советуют зачем-то гвидобейсик. а производные APL — они, конечно, забавны для матриц. однако ж учиться программировать с их помощью… я бы сказал, что это несколько экстравагантная идея.

хотел бы запрограммировать что-то такое

[freedom_of_sea.livejournal.com]

какую биопрограмму вы бы хотели иметь? Лучше в виде сайта.

[worldmanual.livejournal.com]

круто... Хорошо, что заметили))

а я-то все мечтаю о временах, когда можно будет загнать все результаты в супер-компьютерный-мозг, чтобы он там все пообмыслил и выдал закономерности, незамеченные пока людьми =)

при этом не знаю, как оценивать достоверность результатов из вполне серьезных статей.

ну, предположим, процентов на 30 достоверны. Еще 30% - злонамеренный обман (от легкой подтасовки "надо защититься/выступить на конфере с красивой кривой" до высосанных из пальца историй ГМ-алармистов). А еще 30% - ошибки типа вот этой... Еще 10% надо куда-то пристроить, ну да неважно)).

вопрос: а по-вашему, как эти проценты распределяются?..

еще раз, из любопытства, а без ссылок виден?

[tiy.livejournal.com]

На мой взгляд, это совершенно стандартная проблема с нулевыми ячейками. Если считать p-value в любом исследовании банальным методом chi-square test, то при наличии пустых ячеек к ним добавляют какую-нибудь цифирь поменьше (и называется это, если я не ошибаюсь, Haldane correction).

Разница прочиток (reads, как это будет по-русски?!) между 0 и 1 Вас, конечно, не интересует, а вот 0 и 27 уже, возможно, указывает на тканевую специфичность (зависит от покрытия (coverage), трудно сказать 27 это много или мало, если не знать как обстоят дела в среднем по госпиталю). Так что совсем выкидывать эти результаты я бы не стала. Можно установить границу на минимальную экспрессию в одной из тканей, и далее выкидывать все, что не соответсвует критериям.

В качестве способа иллюстрации уровня экспрессии наши товарищи используют в картинках дополнительную вертикальную ось, где отмечают уровень экспрессии определенной пробы по сравнению со всеми другими тестируемыми пробами - если открыть вот эту страничку - (ссылку удалила, Dubois et al, Nat Genet. 2010 Apr;42(4):295-302; supplementary figures) и поискать в тексте supplementary figure 2, то под шапкой к рисунку будет много картинок, где справа на оси Y есть "риска" для каждой пробы. Чем она выше, тем выше экспрессия этой пробы по отношению ко всем остальным. В данном случае использовалось не секвенирование а genome-wide illumina gene expression chip, но можно придумать что-то подобное и в Вашем случае секвенирования.

Год назад в Nature была пара статей про RNA-секвенирование в качестве оценки экспрессии, возможно, там найдется полезная информация на тему того, какой они использовали cut-off для анализа экспрессии:
Pickrell et al, Nature 464, 768-772
Montgomery et al, там же

(статья Пикрелла/Притчарда несколько более разумная, я бы с нее начинала). Впрочем, сейчас довольно много печатают на эту тему, из того что я помню - неплохое ревью на тему RNA секвенирования было в январском номере Nat.Rev.Genetics.

К вопросу об анализе данных - в Гронингене на кафедре генетики есть несколько открытых позиций (PhD, postdoc) в проектах связанных с экспрессией генов. Так вот на них с гораздо большей радостью возьмут математиков, чем биологов. Все-таки биологии легче научиться по ходу дела, чем математике и программированию (хехе, это мой хронический комплекс - почему я не программист...).

Успехов! :)

Re: еще раз, из любопытства, а без ссылок виден?

[agasfer.livejournal.com]

"Разница прочиток (reads, как это будет по-русски?!) между 0 и 1 Вас, конечно, не интересует, а вот 0 и 27 уже, возможно, указывает на тканевую специфичность (зависит от покрытия (coverage), трудно сказать 27 это много или мало, если не знать как обстоят дела в среднем по госпиталю)."

Это указывает прежде всего на большой разброс. Пока не сделаны множественные копии анализа уровней экспрессии с каждого образца данных, никаких выводов делать нельзя, тем более что 27 копий это близко к порогу чувствительности.

"Так что совсем выкидывать эти результаты я бы не стала."

А я бы стал. Это примерно как колориметрия вблизи стандартного буфера: сравнивают оптическую плотность 0.001 и 0.003, и говорят, "Видите? Во 2м образце в 3 раза больше!"

[blog.mykola.org (from livejournal.com)]

А чому результати порівняння надано саме у такому вигляді? Чому саме логарифм відношення?

[0lenka.livejournal.com]

С днем рождения :)))

[progenes.livejournal.com]

дякую!