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

Для анализа файла с чтениями, картированными на геном митохондрии и хлоропласта резуховидки был создан файл в формате .bcf:

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

Для получения списка SNP и инделей использовалась программа bcftools:

bcftools view -vcg 1.bcf > 1.vcf

На выходе получился файл со списком. Число строк, содержащих 'DP='соответствует числу SNP, а содержащих 'INDEL;' - числу инделей.

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

Обнаружено 641 SNP и 302 инделя.

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

С помощью пакета Velvet были собраны геномы хлоропласта и митохондрии. N50 - это характеристика сборки генома, показывающая длину контига, при которой 50% гипотетической длины последовательности покрываются контигами длины, равной N50. Для этого использовались следующие команды:

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

В папке velveth_dir25 есть файл stats.txt с необходимой информацией о контигах. Найдем 10 наибольших значений длин контигов. Сами контиги находятся в файле contigs.fa. 10 самых длинных я выделил в отдельный файл.

Далее был сделан локальный blast по реферсному геному:

makeblastdb -in organellas.fasta -dbtype nucl
blastn -task megablast -query longest.fasta -db organellas.fasta -out blast.txt -outfmt 7

Информация из файла, выданного на выходе blastn, была занесена в таблицу 1.