На главную

Практикум 14

Часть 1 : Подготовка и очистка чтений

Команда	Эффект от выполнения
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 SRR4240359.fastq SRR4240359_deladap.fastq ILLUMINACLIP:adapters.fa:2:7:7	Перед началом работы из всех файлов с адаптерами в папке /P/y16/term3/block3/adapters создается единый файл с адаптерами adapters.fa. Затем была запущена программа Trimmomatic с шагом ILLUMINACLIP, которая обрезала адаптеры (записанные в файле adapters.fa) с концов чтений. Параметры 2:7:7 соответственно обозначающие максимальное число несовпадений для выравнивания адаптера и кусочка рида, минимальное качество выравнивания между ридом и двумя адаптерами для парноконцевых чтений (чтобы выравнивание было признано верным и кусочки рида были обрезаны как адаптеры), какое минимальное качество выравнивания должно быть между ридом и адаптером, если риды непарные (опять же чтобы обрезать выровнявшиеся кусочки). Файл содержащий последовательности ридов с обрезанными адаптерами - SRR4240359_deladap.fastq Изначально ридов: 13557938 Осталось: 13502066 (99,59%) Было удалено: 55872 (0,41%) (видимо состояли только из адаптеров) Размер файла до удаления адаптеров: 1442035098 Мб, размер файла после удаления адаптеров: 1435577866 Мб
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 SRR4240359_deladap.fastq SRR4240359_processed.fastq TRAILING:20 MINLEN:30	Программа принимает на вход образованный на предыдущем этапе файл с последовательностями ридов из которых удалены адаптеры, обрезает с конца чтений нуклеотиды с качеством меньше 20 и убирая риды, длиной меньше 30. Результат сохраняется в файл SRR4240359_processed.fastq Изначально ридов: 13502066 Осталось: 12549389 (92,94%) Было удалено: 952677 (7,06%) Размер файла до чистки по качеству: 1435577866 Мб, размер файлов после чистки по качеству: 1317276262 Мб
velveth ./ 29 -short -fastq SRR4240359_processed.fastq	Программа принимает на вход очищенные от адаптеров и фильтрованные чтения и создает в текущей директории несколько файлов:Roadmap, Log, Sequences. Суть работы программы заключается в создании хэш таблиц с k-мерами, в моем случае k-меры имеют максимально возможную длину 29=30-1.
velvetg ./	Программа принимает созданные в текущей директории на предыдущем шаге файлы Roadmap, Log, Sequences и собирает согласно созданном на предыдущем шаге хэш-таблицам риды в контиги. Помимо файла с контигами contigs.fa, создаются другие файлы Graph, PreGraph, LastGraph, описывающие процесс построения контигов при помощи графов а так же файл stats.txt. В дальнейшей работе с контигами учитывались те, которые попали в файл contigs.fa. (так как на них стоит порог по длине, в этот файл судя по help не должны попадать контиги короче удвоенной длины к-мера, однако по факту почему то наименьшие контиги, попавшие в это файл имеют длину неудвоенного к-мера.

Всего получилось 647 контигов, N50 - 47361bp (сумма длин всех контигов больше или равных этого больше половины суммы длин всех контигов сборки) Таблица с распределением контигов по длине показывает что большая часть контигов меньше 1000 пар нуклеотидов. Три наибольших контига имеют длины и покрытия 78160 bp, 53.878651; 83127 bp, 51.228207; 91528 bp, 57.461807 соответственно.

Распределение числа получившихся контигов по покрытию

На вышепредставленной диаграмме, отобращающей число контигов, попадающих в определенные интервалы по покрытию видно, что самым типичным является покрытие 2-6 ридов на нуклеотид контига, что очень немного.

Распределение длин контигов в зависимости от их покрытия

Из вышепредставленной диаграммы ясно, что более длинные и осмысленные контиги имеют бОльшее покрытие, чем более короткие и менее осмысленные, но так как маленькие контиги составляют большую часть всех собравшихся контигов, их маленькое покрытие является типичным. Если не рассматривать малые контиги, становится понятно, что большая часть осмысленных имеет очень похожие значения покрытия в районе 50-70. Контиг длиной 2657 и покрытый лучше всех ~ 119 ридами вовсе не картируется на геном. Еще один более-менее длинный контиг 3804 покрытый 69,562042 ридами на нуклеотид в среднем хорошо ложится на геном единым куском с 75% идентичных нуклеотидов на участок 104175-107985 и судя по геномному браузеру ncbi на этом месте находятся гены ферментов а не фаговые вставки. Еще один хорошо покрытый ридами (65,888748) и осмысленный относительно его длины (18588) контиг тоже неплохо ложится на геном, так как он ложится всего лишь тремя кусками с 83% 74% 77% идентичности на участки 35124-44693 ; 26737-28363; 30013-32745 соответственно (часть контига в целом участвующая в выравнивании среференсом - 74%, 83% нуклеотидов участвующих в выравнивании идентичны референсным). Хорошо заметный на диаграмме пик довольно длинного (7738) но плохо покрытого (32.657406) контига соответствует последовательности, которая не ложится на геном при использовании megablast.

Такие короткие контиги, на мой взгляд могли получиться по двум причинам:
1) Покрытие генома ридами оказалось неравномерным и те места генома, которые оказались покрыты ридами плохо собрались в короткие контиги. (при прохождении через граф де Брейна программа не находит других ридов, содержащих конкретное "слово" и обрывает контиг)
2) Короткие уникальные последовательности в геноме оказались разделены повторами и программа сборщик не смогла разрешить эти повторы

Характеристика выравниваний

Контиг #1

Самый длинный контиг, длиной 91528 (bp) и со средним покрытием на нулеотид 53.878651 ложится на геном не очень хорошо, так как в результате получается 18 отдельных выравниваний. При этом в выравнивании с референсом задействован 51% контига (query cover) и 74% нуклеотидов контига, задействованных в выравнивании с референсом, идентичны референсным нуклеотидам. Такое кусочное соотнесение контига с банковским геномом говорит о том, что в контиге присутствуют участки, отсутствующие в референсной хромосоме что говорит о наличии геномных перестроек, которые довольно часты у бактерий (в собранном нами геноме есть какие то вставки/делеции относительно рассматриваемого референса). Однако участки контига, выровнявшиеся с референсом, расположены в контиге том же порядке, что и в референсе.

Номер выравнивания	Координаты в референсе (соответствующие тем, что записаны для контига)	Координаты в контиге	Длина выравнивания	Процент идентичных нуклеотидов	Число мисматчей	Число гэпов	e-value
1	445895-442877	3379-6406	3054	80%	542	61	0.0
2	442817-441135	6536-8197	1694	79%	326	29	0.0
3	440944-440755	8333-8514	190	89%	13	8	1e-58
4	440732-440652	8604-8685	82	98%	1	1	5e-32
5	438139-435267	11310-14233	2954	77%	566	111	0.0
6	435241-431839	14295-17724	3495	75%	710	157	0.0
7	429483-425412	20326-24404	4151	76%	837	151	0.0
8	421327-420477	28514-29367	857	79%	170	9	2e-165
9	417677-413081	32299-36893	4665	80%	817	138	0.0
10	412321-410512	37980-39794	1828	81%	320	31	0.0
11	406218-404050	44042-46210	2177	82%	368	16	0.0
12	403823-398904	46453-51436	5031	80%	868	158	0.0
13	398726-389348	51681-61122	9592	74%	2123	363	0.0
14	386887-385425	63845-65312	1478	79%	288	25	0.0
15	385420-384182	65376-66659	1304	74%	252	85	3e-129
16	375867-374136	75484-77214	1743	82%	290	23	0.0
17	369544-366821	81975-84741	2785	77%	552	79	0.0
18	362954-361926	88647-89672	1038	77%	220	21	3e-159

Наш контиг соответствует участку генома с координатами 361926 - 445895. Средняя длина выровнявшегося участка (конец-начало) - 2623,11. Как мы видим растущие координаты в контиге соответствуют убывающим координатам генома, что отображено на dot-plot в виде убывающей прерывистой линии. Стало интересно, к какому участку генома принадлежит выравнивание с 1% (25 штук) гэпов. Оказалось, что почти все выравнивание (14) попадает на ген аспарагинил-тРНК синтетазы, что и обеспечивает высокую консервативность двух последовательностей (контига и референса) в этом выравнивании. Похожая ситуация соответствия контига конкретному гену наблюдается и для хорошо выровнявшегося контига №4, который соответствует GlntRNA

Изображение почти полного соответствия координат гена аспарагиновой АРСазы и контига, полученное при помощи геномного браузера ncbi

Изображение почти полного соответствия координат гена Gln tRNA и контига №4

Dot Matrix выравнивания контига и референса

Контиг #2

Второй по длине контиг 83127(bp), покрытый в среднем 51.228207 ридами на нуклеотид, ложится похожим образом (так же как контиг №1). Megablast выдает тоже 18 выравниваний с референсом, которые в контиге разделены некими последовательностями, не имеющими сходства с референсом. Обосновать это можно аналогично тому рассуждению, что я написала для первого контига. Однако доля второго контига, задействованного в выравнивании с референсом составляет 71%, что значительно больше, чем для первого контига. Это говорит о том, что в данном случае в сборке присутствует меньше участков, не имеющих аналогов в референсе, это может быть обосновано тем, что данный участок генома подвергается меньшему числу геномных перестроек. Доля идентичных нуклеотидов не сильно отличается от 1 контига и составляет 76%.

Номер выравнивания	Координаты в референсе (соответствующие тем, что записаны для контига)	Координаты в контиге	Длина выравнивания	Процент идентичных нуклеотидов	Число мисматчей	Число гэпов	e-value
1	273055-275551	1-2487	2574	75%	485	164	0.0
2	275566-283706	2546-10767	8396	76%	1596	429	0.0
3	283963-285070	11015-12124	1132	76%	223	46	3e-158
4	285200-286535	12821-14155	1349	76%	297	27	0.0
5	288181-291560	15945-19311	3423	78%	670	99	0.0
6	294227-295755	22733-24259	1534	81%	280	12	0.0
7	295935-303252	24568-31921	7428	77%	1549	184	0.0
8	307878-312179	38755-43066	4370	77%	886	127	0.0
9	312679-315982	43711-47019	3351	77%	678	89	0.0
10	318826-323043	50053-54266	4304	74%	944	176	0.0
11	324747-326950	55586-57792	2240	75%	488	69	0.0
12	327227-330003	58084-60853	2828	76%	571	109	0.0
13	330333-331006	61234-61907	675	83%	115	2	7e-170
14	333222-341508	63898-72248	8459	76%	1773	280	0.0
15	341781-343052	72607-73879	1297	78%	239	49	0.0
16	343228-346547	74025-77364	3391	76%	680	122	0.0
17	348233-349674	78487-79929	1450	78%	297	15	0.0
18	349918-352456	80178-82752	2606	75%	541	98	0.0

Контиг №2 соответствует участку генома 273055-330003. Средняя длина выровнявшегося участка больше, чем для первого контига и составляет 3312,778, что говорит о меньшем числе геномных перестроек для этого участка генома (ведь он меньше по размеру и вероятность геномной перестройки на нем тоже меньше).

Dot Matrix выравнивания контига и референса

Контиг #3

Контиг №3 имеет длину 78160 (bp), и среднее покрытие, большее, чем у 1 и 2 контигов - 57.461807. Этот контиг (вероятно из за большего покрытия и меньшей длины) ложится на геном 7 выравниваниями. Это может говорить о том, что данный контиг в меньшей степени (относительно числа событий) прерывается негомологичными референсу последовательностями, а следовательно содержит меньшее число геномных перестроек (из-за своей длины). Доля контига, участвующая в выравнивании составляет всего 56%, что немногим больше чем для №1 и меньше, чем для №2. К тому же соотнесение контига с геномом в данном случае отличается тем, что первые ~19000 нуклеотидов контига не содержат выравниваний с референсом (что явно не говорит в пользу качества такого контига, ведь отбросив информацию о первых 19000 нуклеотидах контига мы не потеряем никакой информации о его принадлежности к участву генома), в отличии от контигов №1 и №2, содержащих гомологичные референсу участки в своих начале и конце. Процент идентичных нуклеотидов в выравниваниях не отличается от №1 и №2 и составляет 76%.

Номер выравнивания	Координаты в референсе (соответствующие тем, что записаны для контига)	Координаты в контиге	Длина выравнивания	Процент идентичных нуклеотидов	Число мисматчей	Число гэпов	e-value
1	126623-127815	18868-20061	1199	84%	184	11	0.0
2	127825-140555	20126-32866	13010	75%	2711	548	0.0
3	144368-151796	36505-43904	7536	78%	1434	243	0.0
4	153752-161738	45885-53969	8168	78%	1549	264	0.0
5	161898-166752	54174-59034	4914	80%	891	112	0.0
6	166750-173180	59065-65508	6517	76%	1391	159	0.0
7	181712-185289	74226-77841	3652	76%	764	110	0.0

Этот контиг соответствует участку генома с координатами 126623-185289. Средняя длина выровнявшегося участка больше, чем для первых двух контигов и составляет 6333,43. Участки контига, выровнявшиеся с геномом, как и для контигов №1 и №2 идут внутри контига в том же порядке, что и в геноме.

Dot Matrix выравнивания контига и референса

Что происходит со сборкой при уменьшении длины хэш-слова с 29 до 25?

Уменьшается наибольшая длина контига, теперь она 17922, N50 тоже уменьшилось и стало равняться 2972. Распределение контигов по длине практически не изменяется (см таблица).
Три самых длинных контига теперь имеют довольно похожие длины (меньшие, чем для длины слова 29):17146, 17615, 17922, но большие чем для длины слова 29 покрытия (соответственно) 87.504494, 89.201416, 95.343765.

Распределение числа контигов с определенным покрытием

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

Что будет, если отбросить половину ридов?

При отбрасывании половины ридов количество контигов, попадающих в файл contigs.fa заметно снижается: c 647 до и 276 - после. Однако судя по таблице распределений длин контигов, заметное снижение наблюдается только среди контигов длиной меньше 1000 нуклеотидов. Длинных контигов как было порядка 30 штук, так и осталось, что логично, ведь если они изначально собирались длинными, значит были хорошо покрыты ридами и при удалении половины ридов существует вероятность, что оставшихся будет достаточно, чтобы собрать длинный контиг. Маленькие же контиги обычно покрыты плохо и поэтому убирая половину ридов, мы мешаем таким ридам собираться.
N50 при отбрасывании половины ридов остается почти таким же и даже немного увеличивается - 52117. Длина сборки (по контигам из файла contigs.fa) уменьшается незначительно c 713529 п.н. до 674523 п.н. и это уменьшение, как мы поняли происходит по большей части за счет отбрасывания коротких, плохо покрытых ридами контигов. Длины и покрытия трех самых крупных контигов составляют в данном случае 111442 п.н 32.386444; 69513п.н 31.142204; 68958п.н 28.292671. Логично что среднее покрытие падает примерно в 2 раза, а вот наибольшая длина контига даже выросла. Это можно связать с тем, что при отбрасывании части ридов, некоторые "слова" при сборке стали уникальными, повторы разрешились и собрались более длинные контиги. Однако судя по координатам в геноме, на которые они ложатся при использовании megablast, ни один из этих трех контигов не является частью или продолжением трех самых длинных контигов, получившихся при использовании всех ридов. 4й же контиг (61413 п.н 29.705778) является укороченной версией третьего по величине контига из обязательного задания. Более которкий вариант полностью ложится на более длинный, оставляя непокрытым один из концов. Это может означать, что непокрытый конец имел плохое покрытие ридами и не собрался при удалении половины из них. В целом я не заметила никакого значительного падения качества сборки.