В этом задании нужно было картировать риды, очищенные в предыдущем задании, на геномы хлоропласта и митохондрии резуховидки Таля. Геномы были объединены в файл organelles.fasta. Для самого картирования использовалась программа BWA, использующая преобразование Барроуза — Уилера. Сначала файл с референсным геномом должен быть проиндексирован командой

bwa index organelles.fasta

Картирование производится командой

bwa mem organelles.fasta Trimmomatic_cleaned.fastq > mapped.sam

(.sam - стандартный формат для откартированных на геном чтений).

Для анализа произведенного картирования использовался пакет программ samtools. Вначале файл переводился в бинарный формат .bam, с которым работают samtools:

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

После этого файл чтения были сортированы командой

 samtools sort mapped.bam mapped.sorted 

и индексированы

 samtools index mapped.sorted.bam 

(последняя команда создает файл mapped.sorted.bam.bai, в котором лежит индекс).

Проиндексированный и сортированный файл был подан на вход команде samtools idxstats, которая показывает, сколько ридов откартировалось на каждый геном. На выходе получается таблица:

 ENA|AP000423|AP000423.1 154478  666916  0
 ENA|Y08501|Y08501.2     366924  73105   0
 *       0       0       3101710 

Первый столбец показывает геномы (первый - хлоропласт, второй - митохондрия), второй - их длину, третий - количество картированных на них ридов. Видно, что хотя геном хлоропласта короче генома митохондрии, количество картированных на него ридов больше. Значит, в образце было больше ДНК из хлоропластов - видимо, просто хлоропластов в растительной клетке было больше, чем митохондрий.

Командой samtools depth было вычислено покрытие на позицию для каждой из органелл:

 samtools depth mapped.sorted.bam > coverage.txt 

Полученный файл: [x]. Среднее покрытие было расчитано в Excel. Среднее покрытие на позицию для хлоропласта составило 428,606, а для митохондрии гораздо меньше - 19,627.