Практикум 15

Подготовка чтений и геномная сборка

Мне достались чтеня с AC проекта SRR4240358.

Я конкатенировала все последовательности адаптеров в один файл, отрезала адаптеры от чтений, триммировала чтения, удалив с конца нуклеотиды с качеством менее 20, и отфильтровала, оставив только чтения не менее 32 нуклеотидов. После я проиндексировала все k-меры длиной 31 и собрала на их основе геном. Код привожу ниже.

cat ../../adapters/* >> ./all_ad.fa
TrimmomaticSE -phred33 SRR4240358.fastq.gz without_ad.fq.gz ILLUMINACLIP:all_ad.fa:2:7:7
TrimmomaticSE -phred33 without_ad.fq.gz fin_trim.fq.gz TRAILING:20 MINLEN:32
velveth velv_dir 31 -fmtAuto -short fin_trim.fq.gz
velvetg velv_dir/
			

На первом этапе очистки было оставлено 10 368 884 из 10 543 839 ридов (174 955 (1.66%) удалено).
На втором этапе очистки было оставлено 8 016 437 из 10 368 884 ридов (2 352 447 (22.69%) удалено).
До обработки файл весил 493 Мб, после первой очистки — 485 Мб, после второй — 257 Мб.

Результат сборки

Файлы, полученные в результате сборки, можно скачать по ссылкам: файл с контигами и файл со статистикой.

N50, выданный velvetg, равен 8600 (в k-мерах, видимо, а в нуклеотидах 8630). Длины трех самых длинных контигов — 19821, 18714 и 16436 (в нуклеотидах 19851, 18744 и 16466). Их покрытия — 29,5, 29,9 и 30,8.

Среди контигов длиной более 1 k-мера были контиги с покрытием от 1 до 553. Не все контиги из файла со статистикой есть в файле с контигами (по информации из мануала, в этот файл входят только контиги длиной более 2k).

Плюс-минус длинные контиги имели покрытие около 30. Я выбрала три контига длиной более 1 k-мера с наибольшим покрытием, чья последовательность была в файле с контигами. Это были контиги с ID 18, 97 и 129 (их длины — 60, 53 и 129 k-меров, покрытия — 412, 405 и 333). Я выровняла алгоритмом blastn последовательности этих контигов на хромосому Buchnera aphidicola. Контиг с ID 18 выровнялся 1 раз, два остальных — более одного раза.

Три самых длинных контига имели ID 56, 34 и 40, их длины — 19821, 18714 и 16436 k-меров. Результаты выравнивания при помощи megablast привожу в таблице 1.

Табл. 1. Результаты megablast.
ID Длина контига в нуклеотидах Количество выравниваний megablast Длины выравниваний Число различий нуклеотидов Число гэпов Координаты участка хромосомы
56 19851 3 8614, 4396 и 4325 1756, 733 и 912 345, 83 и 156 500370 — 508806, 510438 — 514772 и 496111 — 500325
34 18744 6 2220, 3781, 3228, 2530, 1850 и 478 294, 702, 683, 488, 293 и 77 30, 144, 92, 60, 49 и 9 17962 — 20171, 23067 — 26764, 14727 — 17919, 8599 — 11103, 20358 — 22183 и 13994 — 14465
40 16466 2 6962 и 5019 1412 и 991 206 и 164 467412 — 474242 и 462496 — 467424

На рис. 1, 2 и 3 можно видеть dot plot'ы результатов выравнивания контигов с ID 56, 34 и 40 на геном Buchnera aphidicola.
Вообще довольно странно, что контиги настолько непохожи на геном: и ложатся они с промежутками, и идентичность внутри одного выравнивания около 75%, хотя принадлежат они тому же виду.

Рис. 1. Dot plot выравнивания контига с ID 56 на хромосому Buchnera aphidicola.
Рис. 2. Dot plot выравнивания контига с ID 34 на хромосому Buchnera aphidicola.
Рис. 3. Dot plot выравнивания контига с ID 40 на хромосому Buchnera aphidicola.