Картирование на референсный геном

Для начала из файлов с геномом митохондрий и хлоропласт Arabidopsis в fasta-формате был составлен один. Затем с помощью программы BWA нужно было откартировать на эти геномы очищенные чтения из предыдущего практикума.

BWA реализует алгоритм преобразования Барроуза — Уилера. Первым делом надо проиндексировать референсный геном:

bwa index organellas.fasta

Далее собственно картирование:

bwa mem organellas.fasta result.fastq > bwa.sam

Для анализа того, что получилось, надо использовать программу samtools с разными параметрами. Так как она работает только с данными в формате .bam, а у нас .sam, нам необходимо перевести bwa.sam в bwa.bam с помощью программы view:

samtools view -b -S -h bwa.sam > bwa.bam

Теперь надо выяснить, сколько чтений откартировалось на каждую органеллу. Для этого продолжаем использовать samtools. Внутри этой программы есть другая программа idxstats, которую для этого и нужно использовать, предворительно отсортировав (sort) и проиндексировав (index) наше картирование:

samtools sort bwa.bam bwa_sort
samtools index bwa_sort.bam
samtools idxstats bwa_sort.bam > out

Проанализируем полученный результат. На митохондрию (ENA|Y08501|Y08501.2) откартировалось 73186 ридов, а на хлоропласт (ENA|AP000423|AP000423.1) - 669767. Осталось только посчитать среднее покрытие для каждой из органелл. Для этого внутри samtools есть еще одна программа - depth:

samtools depth bwa_sort.bam > cover

Если посчитать, среднее покрытие митохондрии 19,7, а для хлоропласта 431,2.