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

Есть обработанные риды генома Arabidopsis thaliana. Задача - определить, какие из них картируются на геномы хлоропласта и митохондрии.
Файл с геномами хлоропласта и митохондрии: organellas.fasta.

Для картирования используется программа BWA (Burrows-Wheeler Alignment Tool). Она основана на преобразовании Барроуза — Уилера.
Вначале проиндексируем референсный геном:

bwa index organellas.fasta

Алгоритм BWA-MEM рекомендован для чтений Illumina длиной более 70 п.о., поэтому он подходит.

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

Чтобы проанализировать полученные выравнивания, используем программу samtools. Она работает с данными в формате bam. Переведем полученные данные в этот формат с помощью подпрограммы view. Опции: -S для входного файла в формате sam, -b для выходного файла в формате bam, -h для заголовка в выходном файле.

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

Чтобы выяснить, сколько чтений откартировалось на каждую органеллу, воспользуемся подпрограммой idxstats. Вначале нужно отсортировать данные (sort) и индексировать их (index).

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

Результат:

ENA|AP000423|AP000423.1 154478  670023  0
ENA|Y08501|Y08501.2     366924  73370   0
*       0       0       3115103
        

Это значит, что на хлоропласт (ENA|AP000423|AP000423.1) откартировалось 670023 рида, а на митохондрию (ENA|Y08501|Y08501.2) - 73370.

Чтобы выяснить среднее покрытие для каждой из органелл, используем подпрограмму samtools depth.

samtools depth bwa_sort.bam > cover.txt

Получаем, что среднее покрытие для хлоропласта составляет 431, для митохондрии 20.