SNP, индели и сборка

Наша задача - сравнить полученные ранее риды и референсный геном и найти различия (однонуклеотидные полиморфизмы (SNP) и индели). Для сравнения нам понадобится программа mpileup внутри уже известной нам samtools. (флажки для урегулирования форматов выводного и входных файлов). Строка запроса:

samtools mpileup -u -g -f organellas.fasta bwa_sort.bam > comp.bcf

Затем нужно получить список SNP и инделей. Для этого нам понадобится другая программа bcftools:

bcftools view -vcg comp.bcf > comp.vcf

В итоге получаем файл, в котором строки, начинающиеся с "DP=", соответствуют SNP, а строки, начинающиеся с "INDEL;", соответствуют инделям. Посчитаем их:

grep 'DP=' comp.vcf | wc -l
grep 'INDEL;' comp.vcf | wc -l

Итог: 653 SNP и 288 инделей.

Следующий шаг - сборка хлоропласта и митохондрии. Для сборки генома из ридов используют пакет velvet, который состоит из velveth и velvetg. У собранного в результате генома есть важная характеристика "N50", которая показывает длину контига, при которой 50% гипотетической длины последовательности покрываются контигами длины, равной N50. Мне было предложено попробовать изменять параметр "hash_length" в пределах от 15 до 97, чтобы найти ту длину k-мера, при которой N50 будет наибольшей. Итоговая строки запроса (hash_length=25 ; N50=246):

velveth velveth_dir25 25 -fastq result.fastq
velvetg velveth_dir25 -cov_cutoff auto

Из полученного (папки "velveth_dir25") нам нужен файл "stats.txt", в котором содержится информация о контигах, и файл "contigs.fa" с самими контигами, из которых мы отберем 10 по длине в файл "10contigs.fasta". И последнее, что нужно сделать - это blast контигов по нашему файлу с геномом хлоропласты и митохондрии:

makeblastdb -in organellas.fasta -dbtype nucl
blastn -task megablast -query 10contigs.fasta -db organellas.fasta -out blast_result -outfmt 7

Получаем файл "blast_result". Подробнее в таблице 1.