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

[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

[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 и чувствительности метода - я сказать не могу, пробую сама разобраться. Наверное, это можно сказать только глядя весь массив данных и зная особенности технических ошибок.

[tapka.livejournal.com]

Ну видно что в пвэлью они тоже подставляли это свое 0,001.
Вообще не знаю, как тут правильно.
Лучше всего конечно сделать порог по N1 какой-нибудь разумный.
Вообще я с секвенированием сама еще не сталкивалась, не знаю, какие там критерии, можно ли сравнивать гены между собой.
Если да, то может имеет смысл не прямо фолдчендж считать, а какую-нибудь ранговую статистику. Тут нули не должны мешать.
Типа если отсортировать гены по экспрессии в обоих тканях, то какие-то гены будут всегда вверху или внизу, а какие-то поменяют свое положение с сотой позиции на двухтысячную. Возможно при сравнении разных тканей это более перспективно чем фолдчендж, это же не одиночное воздействие, где можно рассчитывать что отработали гены какого-то определенного пасвея, а весь комплекс взаимоотношений.

[progenes.livejournal.com]

Вы говорите-говорите. Я все внимательно читаю и записываю. Если возникнут вопросы - спрошу.

[tapka.livejournal.com]

Ну я не настоящий пожарник, но поговорить могу :) Мне самой интересно как тут правильно поступать.
Кстати.
Бывшие коллеги-программисты вроде именно дип сиквенсингом занимаются сейчас, надо их попытать.

[dr-tambowsky.livejournal.com]

Да всё правильно Вы говорите, какие там такие особенные сложности?? Технологии меняются, голова должна оставаться. В первом приближении задачка ничем не отличается от старого доброго майкроэррэя образца 1914 года. Если на одном эррэе замеренная интенсивность флюоресценции 15, а на другом эррэе мы на той же пробе замерили 0.05 (или, не дай бог 0), то что мы делаем? Считаем log-ratio или проставляем в соответствующей колонке missing value? Ниже предела чувствительности. Или экспериментальный глюк. Или probe failure. Универсальный ответ - NA ;) У эррэя есть минимальная интенсивность, всё что ниже - шум. Почему в RNA-Seq всё должно быть по-другому? Ниже какого-то количества шум всё побеждает. Очень интересно, конечно, что в одном образце мы насчитали 15, а в другом 0, но квантифицировать это довольно сложно.

[bret.livejournal.com]

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

[dr-tambowsky.livejournal.com]

Ах, и Вы, конечно правы... Только определите понятие "уверенно". Структура шума - безусловно другая, если влезать в детали. Ошибки секвенирования, ошибки в alignment, всё это со страшной силой context-dependent - там где homopolymer runs, low complexity/repeat/orthologous sequences всё совсем плохо. И надёжной процедуры для определения mapping quality (вот это самое "уверенно") никто пока так и не предложил, это очень нетривиальная задача на самом деле. Не поймите меня превратно, технология отличная и гораздо "чище" эррэев - на уровне глобальном, когда вы смотрите на expression profile в целом. И dynamic range гораздо лучше. И возможность находить вещи de novo. Но когда Вы добираетесь до отдельных транскриптов и пытаетесь делать какие-то выводы на этом уровне - приходится быть более и более осторожным. Так же как с эррэями - в первом приближении проще отбросить всё что вызывает подозрения, потом уже лезть в детали. При использовании технологий, основанных на секвенировании, подозрений, возможно, возникает количетсвенно меньше, но они есть и принцип тот же.

[bret.livejournal.com]

тут я ни с чем не могу поспорить :)

Re: UPDATE II

[progenes.livejournal.com]

да, в знаменателе. Просто не знала, как вбить формулу.

[progenes.livejournal.com]

Спасибо, добрый человек.
"Если у нас разумный эксперимент, то N1 обычно равно N2 - (они определяются работой секвенатора, в первую очередь)" Судя по всему изначально N1 и N2 не равно, но сырые риды как бы "нормализованы" в показатель RPKM (read per kb per million reads), поэтому в конце концов N1 и N2 приравнены (если я правильно понимаю).
Отдельное спасибо, что объяснили про бином, а от я как баран на факториалы смотрела.
Ок, я могу ранжировать по P-value без проблем. Кроме P-value есть еще FDR для определения трешхолда для P-value. Причем у меня уже были отсортированные по FDR<0.001 и по log2>1 списки генов, которые лихо растыканы в красках по метаболитичесих путях. Но мне все-равно как-то надо интерпретировать, что за разница и на чем она базируется. То, что в одной ткани совсем нет ридов для каких-то генов, а в другой есть, причем это наблюдается на всех стадиях развития, для меня архиважная информация.

вроде бы все остальное правильно

[vasja-iz-aa.livejournal.com]

http://en.wikipedia.org/wiki/Bonferroni_correction

она в их формуле не фигурирует

[vasja-iz-aa.livejournal.com]

а должна

Вот в этом то все и дело. А без нее формула -- бред. Причем не просто смысла не имеющий бред, но заведомо приводящий к неверным выводам бред.

Re: она в их формуле не фигурирует

[progenes.livejournal.com]

Вы себе даже не представляете, КАК я Вам благодарна за Ваши комментарии. Хоть я все еще пока до конца и не поняла. Например, как определить threshold для самих значений FDR. Как и то, что, то есть как определяет FDR достоверность (?) P-value.

Если что, то у меня не просто риды, а "нормализованные" (? ) в reads per kb per million reads (RPKM). Опять таки, навскидку 80% унигенов (конечных кнтигов) с максимальной длиной около 200 баз. Референтного генома, чтобы определить, сколько реально генов и которые из них действительно экспрессия, а не контаминация - нет. Такой короткий контиг тупо считается "геном", от него и пляшут. Насколько к такому "ошметку" "прилипшие" риды отображают реальную картину экспрессии - это мне никакое P-value не скажет. А может это куски консервативных доменов с разных генов? Какие они там gaps допускали? Черти что!

Re: она в их формуле не фигурирует

[progenes.livejournal.com]

Я получила Ваш комментарий и прочитала статью. Ауч! У меня вопрос на засыпку. Если, допустим, сейчас взять и по новой повторить вашу работу, то, положа руку на сердце, насколько вы ожидаете повторяемость вообще? Вы сравниваете две библиотеки. У меня таких 48 постадийно и в разных тканях. И должна сказать, что есть пару фактов для отрезвления.

1. Функциональная аннотация по функциональным категориям и прочим ОГ во всех 48 библиотеках почти один в один. Это отражает только особенности и предел автоматической аннотации вообще. При более внимательном рассмотрении оказывается, что автоматическая аннотация покрывает процентов 60 (это я очень щедра сегодня).

2. Вы берете щедрый допуск E-value для BLASTa. Если я не ошибаюсь, то 10 в минус 6ти. После длительных и продолжительных боев я остановилась на 10 в минус 20 и все еще нахожу всякий мусор.
3. Возвращаясь к болезненной теме репродуктивности ризалтс. Смотрите внимательно на мои результаты:
стадия 1, ткань 1 vs. ткань 2 - 9412 differentially expressed genes (DFG)
стадия 2, ткань 1 vs. ткань 2 - 46347 DFGs
стадия 3, ткань 1 vs. ткань 2 - 7332 DFGs
стадия 4, ткань 1 vs. ткань 2 - 5601 DFGs
стадия 5, ткань 1 vs. ткань 2 - 5719 DFGs

Вам ничего не кажется подозрительным в таком раскладе? Ткани одни и те же, собранные с интервалом в несколько дней. Все библиотеки проходили идентичную нормализацию. Вам не кажется странной стадия 2? Если смотреть с биологической точки зрения, то там не проиходит ничего страшно драматичного в экспрессии, чтобы настолько отражалось в уровне экспрессии генов. Если сравнивать все стадии между собой, то если убрать стадию 2, то получаются логичные кластеры. А стадия 2 как бельмо на глазу, там DFGs показывают противоположную тенденцию. У меня для этого есть одно единственное объяснение - мы наблюдаем техническую ошибку, причем у меня сто идей, откуда она могла взяться. И вместо того, чтобы ее искать, я просто исключу эту стадию. Скажите, у вас есть какие-то критерии, чтобы исключить подобную ошибку в вашей статье при разовом сравнении двух библиотек? В остальном придраться к статье невозможно, чисто сделано и аккуратно.

Re: она в их формуле не фигурирует

[progenes.livejournal.com]

продолжаю изучать статью. От меня также ускользнуло вот что - вы чистые риды пересчитывали на длину потенциального гена, к которому рид относится? Допустим, потенциальный ген длиной в две тыщи баз. У нас два рида, которые на самом деле являются кусками одной и той же молекулы кДНК, которая физически разорвались в процессе приготовления? Мы их посчитаем как два независимых рида, а на самом деле это один.

Отрадно, что в картинке набмер 4 вы ограничились солидными контигами, длиной свыше 300 баз. Я тоже так делаю.

Таблица номер 5 особенно интересна. Я правильно понимаю, что контиги с двумя ридами определялись как кандидаты на потенциально DEG? А контиги с ридами больше восьми дают аж 320 DEG? Неплохо бы проверить гипотезу не статистическими методами, а техническими вроде qRT-PCR. Не удивительно, что мои 27 ридов в контиге показались кому-то странными. Коллеги проверили qRT-PCR и нозерном пару пилотных генов и пришли к неутешительному выводу, что трешхолд у них проходит по границе в... 100 ридов на контиг. Тогда результаты по генной экспрессии не то, чтобы коррелируют, а хотя бы напоминают отдаленно реальную картину с 454 транскриптомиксом. Конечно, это зависит от количества ридов, чистоты секвенирования, но такие факты полезно знать.