Сборка de novo¶

Сборка генома бактерии Buchnera aphidicola de novo.

Риды SRR4240360 для последующей сборки были получены с помощью команды:

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/000/SRR4240360/SRR4240360.fastq.gz

Они находтся в src.

Далее нужно удалить из ридов все возможные адаптеры. Для начала был создан объединенный файл с адаптерами с помощью команды:

cat /mnt/scratch/NGS/adapters/*.fa > adapters.fa

Очистка ридов¶

С использованием trimmomatic необходимо удалить все адаптеры из ридов. Делается это командой:

java -jar /usr/share/java/trimmomatic.jar SE -threads 8 ./src/SRR4240360.fastq.gz SRR4240360_UnA.fastq.gz ILLUMINACLIP:./src/adapters.fa:2:7:7 &> ./log/un_adap.log

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

После удаления адаптеров можно фильтровать риды по качеству (не менее 20) и длине (не менее 32):

java -jar /usr/share/java/trimmomatic.jar SE -threads 8 SRR4240360_UnA.fastq.gz filtered.fastq.gz TRAILING:20 MINLEN:32 &> ./log/trimm.log

Команда удалила 297300 ридов, а размер файлов уменьшился с 202 843 525 байт у исходного до 192 395 357 байт у отфильтрофанного файла.

Таким образом были получены риды без адаптеров длинной не менее 32 и качеством не менее 20.

velveth¶

Почитав -help для velveth были получены k-меры длины 31 и сохранены в директории kmers с помощью команды:

velveth kmers 31 -fastq.gz filtered.fastq.gz -short &> ./log/velveth.log

velvetg¶

Почитав --help для velvetg в директории kmers была создана сборка на основне k-меров, полученных выше. Команда:

velvetg kmers &> ./log/velvetg.log

По результатам работы программы можно сказать, что N50 - 43070.

При помощи pandas был проанализирован файл stats.txt

In [1]:

import pandas as pd

df = pd.read_csv('stats.txt', sep = '\t')
df.head()

	ID	lgth	in	short1_cov	short1_Ocov
0	1	113474	0	33.525460	33.525460
1	2	41715	0	36.311950	36.311950
2	3	11607	0	31.551822	31.551822
3	4	64155	2	35.847323	35.847323
4	5	83603	0	33.646065	33.646065

В таблице ниже представлены 3 наиболее длинных контига, а также их ID и покрытие

In [2]:

df.sort_values('lgth', ascending = 0)[['ID', 'lgth', 'short1_cov']].head(3)

	ID	lgth	short1_cov
0	1	113474	33.525460
4	5	83603	33.646065
3	4	64155	35.847323

Далее представлены 4 контига с аномально большим покрытием.

In [3]:

df.sort_values('short1_cov', ascending = 0)[['ID', 'lgth', 'short1_cov']].head(4)

	ID	lgth	short1_cov
575	576	0	inf
172	173	1	134953.0
393	394	1	590.0
344	345	2	297.5

Интересным является контиг длиной 0. В принципе, логично, что если контиг имеет длину 0, то он подходит везде, однако странно само по себе, что такое бывает. Остальные же покрытия представленных контигов объясняются тем, что они очень короткие (1-2 нуклеотида), а значит - много куда подходят.

Сравнение трех самых длинных контигов с хромосомой Buchnera aphidicola¶

Сравнение проводилось в megablast между хромосомой Buchnera aphidicola (CP009253) и тремя самыми длинными контигами (ID: 1, 5, 4).

ID контига	Max Score	Total Score	Query Cover	Per. Ident	Acc. Len	Выдача
1	17265	51702	76%	81.43	628164	выдача
4	5749	28200	70%	78.38	628164	выдача
5	5465	26995	58%	74.95	628164	выдача

1-ый контиг выровнялся с промежутком от 528794 до 550219 с 545 гэпами (2%)

4-ый контиг выровнялся с промежутком от 2004 дo 11103 с 256 гэпами (2%)

5-ый контиг выровнялся с промежутком от 127825 дo 140555 с 548 гэпами (4%)