Учебная страница курса биоинформатики,
год поступления 2020
constructed
....
3.b Второе правило Чаргаффа
Для применения статистики нужно выдвинуть так называемую нулевую гипотезу о случайности появления буквы A или буквы T из комплементарной пары A-T на одной цепочке ДНК с вероятностями по 1/2 и независимо друг от друга. Т.е. принять распределение Бернулли и вычислить вероятность случайного появления наблюдаемого или большего отклонения от ожидаемого - число A = число T. См ниже как.
4. GC состав генома
В статье https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3053387/ в Table 1 перечислены GC составы некоторых видов бактерий из разных отделов бактерий. В обсуждении сравните GC состав вашей бактерии с таковым у других представителей того же отдела.
5. Анализ k-меров с маленьким k
Найдите k-меры (слова длины k) экстремальные по отклонению от ожидаемого, недопредставленные и перепредставленные. Пороги, принимаемые на основе опыта: cb < 0.8 - слово недопредставлено, cb > 1.2 - перепредставлено. Математического обоснования нет, так что эти пороги не более, чем орентировочные. Для экстремалов по cb надо искать причину. Иногда известна, чаще - нет.
6. Найдите и опишите повторяющиеся последовательности в геноме, появление которых нельзя объяснить случайностью
Нельзя объяснить случайностью совпадение длинных слов. Например разных слов длины 20 всего 420 = 240 примерно 10004 = 1012. В геноме число слов длины 20 примерно равно его длине N, т.к. слово может начинаться с 1 нукл, 2го, 3го и т.д. до (длина генома -20)го. Вероятность совпадение двух слов при N имеет порядок (1/10^12)*N, т.е очень мала.
Найдем длинные слова длины 30, встречающиеся в геноме не менее 20 раз! В моем примере числа 30 и 20 подобрал путем перебора, чтобы получить не слишком много находок, но более одной:
wordcount -wordsize 30 -mincount 20 -srevers1 neverov_genome.fasta stdout
stdout можно заменить на имя выходного файла
wordcount для слов такой длины работает заметное время, много минут.
Выберем одно из найденных слов и найдем координаты всех таких слов в геноме:
fuzznuc -pattern GTTTGTAGCTTACCTATAAGGGATTGAAAC -srevers1 neverov_genome.fasta stdout
Получили результат, довольно удивительный в этом примере!
Можно разрешить, например, два несовпадения в слове с геномом с тем которое ищется:
fuzznuc -pattern GTTTGTAGCTTACCTATAAGGGATTGAAAC -srevers1 -pmismatch 2 neverov_genome.fasta stdout | wc -l
wc поставил в pipeline чтобы посчитать число находок.
Число находок увеличилось - стало 194 а было 129. Случайное совпадение остается маловероятным.
Далее что интересно. Сравнить координаты находок анализируемого слова с координатами генов. Если внутри генов, то это дупликации генов? Каких?
Если нет, то на одной ли цепочки или на разных? Идут ли они близко расположенными парами? Инвертированные повторы (если в паре они на разных цепочках) или тандемные повторы (на одной цепочке). Всякие такие повторы крайне интересны, даже если объяснить их вам не удается.
Так открытие знаменитых CRISPR-Cas систем, использующиеся для генной инженерии и даже генной терапии людей, началось с того, что в 90х годах японский учёный (сейчас не помню фамилии, надо посмотреть) обнаружил в геноме повторяющиеся последовательности - странные потому что не мог их объяснить. Потом их назвали CRISPR (R от repeat), лет через 10 кое какое объяснение было найдено, и еще лет через пять стали использовать для генной инженерии. Дата по памяти, могу ошибиться +/- пять лет. Если спросите - посмотрю и отвечу
Старые подсказки
Постройте гистограмму числа "квазиоперонов" по числу генов. У бактерий и архей оперон - участок ДНК с одним или несколькими генами белков, транскрибируемый в одну матричную ДНК. Таким образом, гены в одном опероне закодированы на одной цепочке ДНК. Как правило, расстояние между ними небольшое.
Иногда используют простейший способ предсказания оперонов. "Квазиопероном" назовем максимальную последовательность генов, закодированных на одной цепочке с промежутками между генами не более порога, например, 100 п.н. Квазиоперон может состоять и из одного гена.
Постройте гистограмму числа квазиоперонов по числу генов в квазиопероне в вашем геноме.
Постройте гистограмму числа "квазиоперонов" в зависимости от порога на расстояние.Советую сделать ячейку с параметром порог длины с числом 100. Тогда изменение числа квазиоперонов при изменении порога получается изменением значения этого параметра.
Число генов в квазиопероне легко посчитать с помощью СЧЁТЕСЛИ. И гистограмму недолго построить.
Проверьте гипотезу о том, что гены распределены по цепочкам случайно с вероятностями 0,5 Пусть в вашем геноме 5000 генов белков и 2450 - на одной цепи, 2550 - на другой. Отклонение от ожидаемого числа 2500 равно 50. Надо ответить на вопрос возможно ли получить такое или большее отклонение при независимом случайном выборе цепочки для каждого гена.
- Подбросим монетку 5000 раз, посчитаем сколько раз выпал орел, каково отклонение от ожидаемого 2500 и сравним с наблюдаемым нами отклонением 50. Если больше или равно - ставим плюс, если меньше - ставим минус.
- Повторим эксперимент, допустим, 100 раз (а лучше - 1000).
- Пусть из 100 экспериментов только один раз отклонение оказалось больше или равно наблюдаемому нами 50. Значит, вероятность увидеть отклонение 50 или больше примерно равна 1/100 = 0.01. Вывод. Если мы готовы считать, что событие с вероятностью 0.01 маловероятное, то полученное нами отклонение 50 противоречит гипотезе о независимом случайном равновероятном выборе цепочки для гена. Значит, надо искать причины.
Первое испытание - в колонке 1. Используйте СЛУЧМЕЖДУ нулем ("решка") и единицей ("орел"). Функция выдает 0 или 1 с равной вероятностью.
Распространите формулу вниз столько раз, сколько генов в вашем геноме. В этом же столбце (например, в верхних ячейках) рассчитайте число орлов (СЧЁТЕСЛИ) и отклонение числа орлов от ожидаемого - без знака!.
Распространите все формулы в сто соседних столбцов. Посчитайте сколько раз отклонение больше или равно тому, которое вы обнаружили в своем геноме.
Сделайте вывод.
Замечание В Excel есть формула, которая сразу выдает нужный вам ответ - вероятность получить такое же или большее число отклонений от среднего, какое вы обнаружили. Если вы что-то знаете по теории вероятности, то можете ее найти и применить. При применении тоже есть некоторые фокусы, которые - не понимая что делаете - вы не сумеете учесть! Поэтому ...
Мой мальчик! Тебе эту песню дарю. Рассчитывай силы свои. И, если сказать не умеешь "хрю-хрю", - Визжи, не стесняясь: "И-и!" С.Маршак
Представьте статистические данные о пересечениях генов (если таковые обнаружатся). Постройте гистограмму длин пересечения генов белков в вашем геноме пересекающихся
Длинное пересечение генов - удивительная вещь, но в природе встречается. Также случаи пересечения генов могут быть результатом ошибки в предсказании кодирующих последовательностей. Для одного или трех типов пересечения голова к хвосту, голова к голове, хвост к хвосту. Отметить гены, пересекающиеся с предыдущим, можно в новой колонке с помощью ЕСЛИ. В следующей колонке можно вычислить сдвиг рамки и ориентацию пересекающихся генов друг относительно друга. Придумайте, как это сделать!
Если сделали, то посчитать число пар пересекающихся генов можно с помощью СЧЁТЕСЛИ.
(Самое важное и трудное) Представьте статистику белков по категориям достоверности их существования. Cведения о том, каким образом подтвеждено существование гена можно получить только из базы данных белков Uniprot
В таблице с генами есть колонка product_accession. Зайдите на сайт Uniprot и выберите Retrieve/ID mapping. Этот сервис служит для перекодировки из одной системы идентификации в другую. Отфильтруйте идентификаторы кодирующих последовательностей (а не РНК - у РНК нет "продуктов"), скопируйте колонку идентификаторов и вставьте в окно Provide your identifiers. Выберите FROM: EMBL/GenBank/DDBJ CDS, TO: Uniprot KB => Go.
Получите таблицу, которую можно скачать в формате Excel. Однако сначала надо отредактировать колонки таблицы => Columns. Оставьте колонки Entry name, добавьте Protein Existence (из колонки Miscellaneous), Length, Protein names.
Скачайте в формате Excel, скопируйте на страницу своего файла, и сделайте сводную таблицу по полю "Protein exsistence". Прочитайте где-нибудь что значит каждая из категорий.
Сведения о том, каким образом подтвеждено существование гена можно получить только из базы данных белков Uniprot. Как это сделать - см. в подсказках.
Тот же финт с прекодировкой можно применить к столбцу GeneID, выбрав соответственно БД GeneID (Entrez Gene) в окошке From. (Enrez - так называется совокупность баз данных и сервисов на странице NCBI)