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

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

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

Чтобы сравнить риды с референсным геномом, используем подпрограмму samtools mpileup. Опции: -g, -u для файла в формате bcf в несжатом виде, -f для входного файла референсного генома в формате fasta.

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

Чтобы получить список SNP и инделей, используем bcftools.

bcftools view -vcg 1.bcf > 1.vcf

Результат - файл со списком 1.vcf. Число строк, содержащих 'DP=' соответствует числу SNP, а содержащих 'INDEL;' - числу инделей.

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

Обнаружено 649 SNP и 287 инделей.

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

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

velveth velveth_dirXX XX -fastq out.fastq
velvetg velveth_dirXX -cov_cutoff auto

XX - длина k-mer. После работы программы velvetg она выдает на экране значение N50. Необходимо найти максимальное значение между 15 и 97. Запустим команды несколько раз. Было предложено стартовать от длины k-mer, равной 35. На это программа отреагировала так:

Velvet can't handle k-mers as long as 35! We'll stick to 31 if you don't mind.

Возражать было нечем, пришлось согласиться.
Максимальное N50 равно 245 при длине k-mer, равной 25.

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

Теперь нужно сделать локальный бласт контигов по файлу organellas.fasta, в котором находится референсный геном хлоропласта и митохондрии. Данные о 10 контигах и их расположении в референсном геноме приведены в таблице 1.

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

Результат: blast.txt.

Таблица 1. 10 самых длинных контигов, найденных Velvet