ФББ 2013-2014

Однонуклеотидные полиморфизмы, индели и сборка de novo

Поиск однонуклеотидных полиморфизмов и инделей

Однонуклеотидные полиморфизмы - это единичные отличия в ДНК, возникающие в результате точечных мутаций. Именно они отличают геном двух представителей одного вида. Индель - это инсерция (т.е. вставка) или делеция ("выпадение" или удаление) участка в геноме. В данном практикуме индели и SNP будем искать в ридах, картированных на геномы хлоропласта и митохондрии Arabidopsis thaliana. Для этого используют команды из пакетов samtools и bcftools.

samtools mpileup -u -g -f sequence.fasta map.bam > map.bcf - создание .bcf файла. Здесь sequence.fasta - файл с геномами хлоропласта и митохондрии (референсные), а map.bam - бинарный файл с картированием ридов на геномы.

После создания файла .bcf можем работать с пакетом программ bcftools. Индели и SNP ищутся командой: bcftools view -v -c -g map.bcf > indel.vcf. Получим файл формата .vcf, где в каждой строчке будет стоять информация об инделе или SNP. Теперь посчитаем отдельно количество инделей и SNP, с помощью команд grep и wc -l. Результаты подсчёта: инделей - 287, SNP - 649.

Сборка хлоропласта и митохондрии

Для того, чтобы собрать геном, нужно сначала собрать риды в контиги. Для этого используем программы velveth и velvetg. Velveth - вспомогательная программа для velvetg, на вход ей подаются риды и длина желаемого k-мера, по этим данным она строит базу данных k-меров. Velvetg строит граф де Брёйна и соединяет риды в контиги. Выполним команду для нескольких длин k-меров. При попытке поставить значение выше 31 программа говорит, что ей слишком сложно, и ставит длину 31 (видимо, максимально возможная?). При попытке выполнить команду для k-мера чётной длины, берётся ближайшее нечётное значение. Начнём с максимальной длины k-мера: для него получается N50 = 280. Далее будем уменьшать значение (возьмём 25, 19, 11 и 29). Стоит отметить, что N50 во всех случаях оставалось одинаковым, что странно. Остановимся на длине k-мера в 29 нуклеотидов. Команды для этого значения запишутся так:

		velveth velveth_dir 31 -fastq  trim.fastq
		velveth_dir -cov_cutoff auto
	

Программа создаёт директорию velveth_dir, в которой лежат файлы с графами, контигами, референсной последовательностью, а также файл stats.txt, в котором хранится информация о полученных контигах. Теперь построим выравнивание контигов на последовательность хлоропласта и митохондрии. Сделаем это при помощь локального бласта. Сначала создадим базу данны с помощью команды makeblastdb -in sequence.fasta -dbtype nucl , затем проведём выравнивание командой blastn -task blastn -query contigs.fa -db sequence.fasta -outfmt 7 -num_alignments 1 -out alignment.fa .

Файл с выравниванием находится в директории /ngs/partyhard/pr11 и не приводится здесь из-за слишком большого объёма. Теперь отберём 10 самых длинных контигов с помощью файла со статистикой и Excel. Импортируем stats.txt в Excel, сортируем по столбику с длиной. Данные о самых длинных контигах приведены в таблице: