Посчет количества SNP и инделей

С помощью samtools и bcftools я подсчитала количество SNP (однонуклеотидных замен) и инделей (вставок и делеций) в картированных в предыдущем практикуме ридах. Для этого выполнили команды:

 samtools mpileup -ugf organelles.fasta mapped.sorted.bam > mapped.bcf
 bcftools view -vcg mapped.bcf > mapped.vcf	

Первая команда создает bcf-файл со списком инделей и SNP, вторая переводит бинарный bcf в формат vcf. В файле vcf в каждой строчке записан либо SNP (обозначаемый 'DP'), либо индель ('INDEL'). Далее командой grep и wc -l можно подсчитать количество SNP и инделей в файле.

Количество SNP составило 675, а количество инделей - 292.

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

Используются программы пакета velvet. Сборка производится из очищенных ридов позапрошлого практикума.

velveth - программа, подготавливающая базу k-меров для программы velvetg. На вход подается название выходной директории, длина k-мера (hash length) и формат и название входного файла с ридами. Длина k-мера должна быть нечетной для избежания палиндромов, не выше 31 (т.к. записывается на 64 битах) и желательно выше 10 для достоверности. Я попробовала 3 длины: 31, 25 и 19. Например, для длины 25 команда выглядит так:

 velveth velveth_dir2 25 -fastq Trimmomatic_cleaned.fastq 

Далее с для каждой производилась команда:

 velvetg velveth_dir -cov_cutoff auto 

• Для k=31: Final graph has 335850 nodes and n50 of 216, max 5493, total 40556793, using 0/3841341 reads
• Для k=25: Final graph has 393027 nodes and n50 of 242, max 5487, total 46640785, using 0/3841341 reads
• Для k=19: Final graph has 842802 nodes and n50 of 157, max 2162, total 47889211, using 0/3841341 reads

Для дальнейшей работы был выбран граф, построенный на k-мерах длиной 25, т.к. для них получилась наибольшее N50 - 242 (т.е. 50% длины последовательности покрываются контигами длиной не менее 242).

В Excel был импортирован файл stats.txt из созданной velvet директории, и сортирован по длине контига (нужно добавить, что настоящая длина контига будет на k-1 больше, чем приведенная в stats.txt длина). 10 самых длинных были выделены из файла contigs.fa в отдельный файл: [x] и по базе, сделанной из геномов митохондрии и хлоропласта был произведен локальный бласт этих длинных контигов. Как ни удивительно, все они оказались из митохондрии.