Cначала геномы хлоропласта и митохондрии резуховидки были скачаны и объеденены в общий файл arabidipsis_chl_mito.fasta (ссылка на файл не предоставлена, не могу скопировать в свою директорию "Disk quota exceeded"). Затем очищенные чтения были откартированы на геномы митохондрии и хлоропласта с помощью программы BWA.
Для работы программы необходимо проиндексировать геномы с помощью команды:
bwa index arabidipsis_chl_mito.fasta
Использовался алгоритм BWA-MEM, так как он предназначен для чтений длиной от 70 пар оснований и выше (см ссылку на результат работы FastQ. Почти все риды длиной 90-100 нуклеотидов ) и подходит для секвенирования с помощью Illumina. В результате картирование было произведено при помощи команды:
bwa mem arabidipsis_chl_mito.fasta Ath_tae_CTTGTA_L003_R1_002_trimmed.fastq> bwa_align.sam
Затем использовалась программа samtools для интерпретирования результатов, полученных на предыдущем шаге. Для того, чтобы выяснить, сколько чтений откартировано на каждую органеллу, были выполнены следующие команды
1) файл был переведён из формата .sam в .bam. (-b выход в формате BAM, -h включает заголовок на выходе, насколько я поняла отсутсвие заголовков вызывает трудности при работе с многими программами, -S на вход подан файл в формате SAM)
samtools view -b -S -h bwa_align.sam > bwa_align.bam
2)затем отсортирована (sort) информация об откартированных ридах
samtools sort bwa_align.bam bwa_align_sorted
(если ввести samtools sort bwa_align.bam bwa_align_sorted, то на выходе файлы bwa_align_sorted.0000.bam и bwa_align_sorted.bam.bam)
2) проиндексирована (index)
samtools index bwa_align_sorted.bam
2) получена нужная статистика (idxstats)
samtools idxstats bwa_align_sorted.bam
После выполнения последнего шага, на экран была выведена следующая статистика:
ENA|AP000423|AP000423.1 154478 665652 0 ENA|Y08501|Y08501.2 366924 72779 0 * 0 0 3099687
Каждая строка содержит информацию о референсной последовательности AP000423.1 - геном хлоропласта, Y08501.2 - митохондриальный геном), длина референсной последовательности и количество ридов, откартированных на неё, количество неоткартированных ридов (на практике это значение нужно только в последней строчке).
Полученные результаты показывают, что в 2.5 раза более короткий геном хлоропласта (154478 bp), откартировась значительно больше (примерно в 9 раз) чтений: 665652, чем на геном митохондрии(366924 bp): 665652 ридов . Неоткартированных ридов осталось 3099687. Из этого можно сделать вывод, что при проведении эксперимента секвенировали ДНК из клеток зелёных частей растений, так как в них больше хлоропластов, чем митохондрий.
Покрытие каждого нуклеотида было рассчитано с помощью команды:
samtools depth bwa_align_sorted.bam > nucl_cover
На выходе был получен файл nucl_cover, каждая строка которого содержит название референсной последовательности, на которую было произведено картирование, порядковый номер нуклеотида этой последовательности и его покрытие (количество закартированных ридов, покрывающих этот нуклеотид). С помощью LibreOffice Calc было посчитано среднее покрытие для нуклеотидов хлоропласта (427,7) и митохондрии (19,5).