Для загрузки данного мне файла с чтениями была использована команда:
wget "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/008/SRR4240358/SRR4240358.fastq.gz" -O SRR4240358.fastq.gz
Далее я решил создать файл в формате fasta, содержащий все возможные последовательности адаптеров, которые необходимо вырезать из предоставленных прочтений. Это было сделано следующей командой:
cat /mnt/scratch/NGS/adapters/* > adapters.fasta
Для сборки генома нужно предварительно подготовить чтения. Сначала удалим адаптеры.
Для удаления адаптеров была использована команда:
java -jar /usr/share/java/trimmomatic.jar SE -threads 8 -phred33 -trimlog LOG/trimmomatic.log SRR4240358.fastq.gz trimmomatic/SRR4240358.fastq.without_adapt.gz ILLUMINACLIP:adapters.fasta:2:7:7
Из 10543839 чтений, поданных на вход, осталось 10368884 (98.34%), число оснований при этом сократилось на 1.76%, как показывает суммирование значений в колонке выходного файла trimmomatic.log, отвечающей за число нуклеотидов, сохранённых для каждого чтения.
Далее были удалены нуклеотиды низкого качества с правых концов, а также последовательности, состоящие из менее 32 нуклеотидов.
java -jar /usr/share/java/trimmomatic.jar SE -threads 8 -phred33 -trimlog LOG/trimmomatic2.log trimmomatic/SRR4240358.fastq.without_adapt.gz trimmomatic/SRR4240358.fastq.correct.gz TRAILING:20 MINLEN:32
В итоге, было удалено 2352447, то есть 22.69% чтений. Размер исходного файла в формате .fastq.gz составлял 470M, после удаление адаптеров он уменьшился на 7М, а отбрасывание некачественных нуклеотидов и коротких чтений позволило получить файл размером 341M.
Первый этап сборки - подготовка k-меров (в нашем случае k = 31). Для этого запустим программу velveth:
velveth Assem_SRR4240358 31 -short -fastq.gz trimmomatic/SRR4240358.fastq.correct.gz
Затем командой velvetg была произведена сама сборка:
velvetg Assem_SRR4240358
Эта команда создала папку velveth и в ней следующие файлы:
velveth Assem_SRR4240358 31 -short -fastq.gz trimmomatic/SRR4240358.fastq.correct.gz
Последний файл содержит информацию о покрытии и длине контигов. В результате получили набор контигов с N50 = 8630. 3 самых длинных котига: 116466, 18744, 19851.Им соответствуют следующие покрытия: 30.79, 29.92, 29.47.
Контиг с длиной 19851 картировался на хромосому тремя участками.Координаты участков: 496111..500325; 500370..508806; 510438..514772. Идентичность для первых двух составляет 75-76%, а содержание гэпов 3-4% , для третьего участка показатели лучше: 81% совпадающих оснований и гэпов 1%.
Второй по длине контиг карировался с 5 разрывами. Координаты участков: 8599..11103; 13994..14465; 14727..17919; 17962..20171; 20358..22183; 23067..26764. Идентичность для всех участков находится в интервале 76-85%, а содержание гэпов составляет 1-3%.
Третий по длине контиг карировался двумя участками. Координаты участков: 462496..467421; 467412..474242. Идентичность обоих участков референсу составила 77%, а содержание гэпов — 3% и 2% соответственно. Данный конти инвертирован (картируется на геном в противоположном направлении).