Для загрузки данного мне файла с чтениями была использована команда:
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/009/SRR4240379/SRR4240379.fastq.gz
Далее я решил создать файл в формате fasta, содержащий все возможные последовательности адаптеров, которые необходимо вырезать из предоставленных прочтений. Это было сделано следующей командой:
cat /mnt/scratch/NGS/adapters/* > adapters.fasta
Для сборки генома нужно предварительно подготовить чтения. Сначала удалим адаптеры.
Для удаления адаптеров была использована команда:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240359.fastq.gz SRR4240359_1.fastq.gz
ILLUMINACLIP:adapters.fasta:2:7:7
В результате я получил следующие результаты. Исходно было 13557938 количество чтений, стало 13502066 (99.59%). То есть адаптерами оказались 55872 (0.41%) чтений.
Далее были удалены нуклеотиды низкого качества с правых концов, а также последовательности, состоящие из менее 32 нуклеотидов.
java -jar /usr/share/java/trimmomatic.jar SE -phred33 SRR4240359_1.fastq.gz SRR4240359_c1.fastq.gz
TRAILING:20 MINLEN:32
Количество чтений заметно сократилось с 13502066 до 12184080 (90.24%). То есть было удалено 1317986 (9.76%).
В итоге исходный файл весил 445 мб, а очищенный - 385 мб.
Первый этап сборки - подготовка k-меров (в нашем случае k = 31). Для этого запустим программу velveth:
velveth Assem 31 -short -fastq.gz SRR4240379_c1.fastq.gz
Затем командой velvetg была произведена сама сборка:
velvetg Assem
Эта команда создала папку velveth и в ней следующие файлы:
contigs.fa Graph LastGraph Log PreGraph Roadmaps
Sequences stats.txt
Последний файл содержит информацию о покрытии и длине контигов. В результате получили набор контигов с N50 = 70607. 3 самых длинных котига: 125674(44.55), 108447 (42.009), 71403 (39.41), в скобках указано покрытие. Также нашел последовательность с аносально высоким покрытием - 411220.00
Контиг с длиной 108447 картировался на хромосому только в одном месте:
127825 - 140555, Identities = 75%, Gaps = 4%
Контиг с длиной 125674 картировался на хромосому сразу на 25 участков:
14727 - 613671, Identities = 83%, Gaps = 1%
Контиг с длиной 71403 картировался на хромосому только в одном месте:
266073 - 273028, Identities = 80%, Gaps = 2%