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

Учебный сайт Софроновой Алины

Для этого практикума нам дан файл c геномом хлоропласта и митохондрии резуховидки Arabidopsis thaliana - genome.fasta, который будет использоваться в качестве референсного генома. Картировать будем очищенные риды из прошлого практикума.

Картирование осуществлялось программой BWA (Burrows-Wheeler Alignment Tool). Для всех алгоритмов работы этой программы необходимо проиндексировать референсный геном - команда bwa index genome.fasta
Для осуществления картирования был использован алгоритм BWA-MEM, выполняющий локальное выравнивание. Он подходит для очищенный ридов с прошлого практикума (их длина 100-101, секвенирование с помощью Illumina). Картирование было произведено командой:

bwa mem genome.fasta Ath_tae_CTTGTA_L003_R2_005_2.fastq > bwa.sam.

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

Команда: samtools view -b -S -h bwa.sam > bwa.bam

После сортируем полученный файл командой sort (сортировка по левым координатам выравниваний) и индексировать с помощью index.
Команда:

samtools sort bwa.bam bwa_sort samtools index bwa_sort.bam

Далее мы хотим узнать сколько чтений откартировалось на каждую органеллу. Для этого воспользуемся командой idxstats, которая выдает статистику в индексированном файле.
Команда: samtools idxstats bwa_sort.bam

На экран была выведена такая информация о статистике:


ENA|AP000423|AP000423.1 154478  667242  0
ENA|Y08501|Y08501.2     366924  72397   0
*       0       0       3111307

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

Команда: samtools depth bwa_sort.bam > cover.txt.

Был получен файл cover.txt. Каждая строчка этого файла содержит название последовательности, на которую было произведено картирование, номер нуклеотида и количество чтений,которое покрывает этот нуклеотид. При помощи Excel я нашла среднее покрытие: для митохондрии - 19,4, для хлоропласта - 427,5. Покрытие довольно большое (особенно для хлоропласта).

Вернуться к 3 семестру