Картирование подготовленных данных на геном
Задание 1: Картирование чтений на референсный геном
Для выполнения этого задания я использовал программу
$ bwa mem -t 12 ../bwa/chr3.fna ../dnaseqreads/SRR10720412_1_paired.fq.gz ../dnaseqreads/SRR10720412_2_paired.fq.gz > aln-se.sam 2> aln-se.log
-t 12 : этот параметр поясняет программе, что надо использовать 12 ядер сервера. Нужно это для ускорения процесса;../bwa/chr3.fna : путь к файлу с референсом, НА который картируем;../dnaseqreads/SRR10720412_1_paired.fq.gz ../dnaseqreads/SRR10720412_2_paired.fq.gz : пути к файлам с парными чтениями, КОТОРЫЕ мы картируем;> aln-se.sam : перенаправление картированных чтений в файл;2> aln-se.log : перенаправление логов в файл.
Задание 2: Анализ sam -файла с выравниванием
Внутри
Собственно внутри файл устроен таким образом: первые 2 строчки являются заголовком (отличить их от остальных можно по наличию "@" в начале строки) и вкратце сообщают нам информацию о референсе, чтениях и программе с командой, которые были использованы для картирования. Если говорить об остальной части файла, то там каждая строка сообщает нам информацию о чтении и о том, как оно картировалось на референс.
| Таблица 1. Информация о первых 11-ти столбцах |
||
|---|---|---|
| Номер столбца | Название столбца | Описание |
| 1 | QNAME | Уникальное имя рида. Если вместо уникального имени рида мы обнаруживаем символ "*", то это значит, что рид есть, но информация о нём недоступна |
| 2 | FLAG | Некая комбинация побитовых (если можно так сказать) меток/"флагов" (на самом деле, вообще не понятно, что это), выраженная в одном числе |
| 3 | RNAME | Имя референсного генома, на который мы все риды картировали. Если вместо имени референса мы здесь увидим символ "*", то это значит, что данный рид не закартировался |
| 4 | POS | Номер самого левого нуклеотида в выравнивании рида с референсом, где нуклеотид №1 - первый нуклеотид референса. Если на этой позиции стоит 0, то это значит, что рид не закартировался |
| 5 | MAPQ | Качество картирования чтения на референс |
| 6 | CIGAR | Отчёт о выравнивании рида с референсом в специальной кодировке. Пример: запись в формате CIGAR "2S5M2D2M" стоит расшифровывать так: 2S - два несовпадения, 5M - 5 совпадений, 2D - 2 делеции, 2M - 2 совпадения, и всё это расположено последовательно, друг за другом |
| 7 | RNEXT | Референсное название самого левого выравнивания последовательности следующего рида образца с референсом, причём если у нас рид последний, то следующим для него станет первый рид. Если оно у нас "*", то информация недоступна. Если оно у нас "=", то RNEXT совпадает с RNAME |
| 8 | PNEXT | Позиция самого левого выравнивания последовательности следующего рида с референсом. Если в этом поле стоит "0", то информация недоступна |
| 9 | TLEN | Длина образца. Если в этом поле стоит "0", то информация об длине недоступна |
| 10 | SEQ | Собственно последовательности рида. Если в этом поле стоит "*", то последовательность рида не сохранилась. А если в этом поле стоит "=", то это значит, что у нас рид полностью совпал с референсом |
| 11 | QUAL | Строка качества каритрования рида на геном (Phred+33) |
Размер
Задание 3: получение и индексирование bam -файла
Для начала я перевёл наш
$ samtools sort -@ 12 -o aln-se.bam aln-se.sam 2> samtoolsort.log
samtools sort - обращение к программеsort из набора утилит для работы с картированными чтениямиsamtools ;-@ 12 - позволение программе использовать 12 ядер для ускорения работы;-o aln-se.bam - название файла, куда скидывать новоиспечённый отсортированныйbam ;aln-se.sam -sam -файл, который мы должны конвертировать и отсортировать;2> samtoolsort.log - перенаправление логов в отдельный файл.
Вес получившегося
Далее
$ samtools index -@ 12 aln-se.bam
samtools index - обращение к программеindex из набора утилит для работы с картированными чтениямиsamtools ;-@ 12 - позволение программе использовать 12 ядер для ускорения работы;aln-se.bam -bam -файл, который мы индексируем.
Задание 4: Анализ bam -файла
Сейчас нам нужно провести статистический анализ отсортированных и проиндексированных картированных чтений. Сделаем мы это с помощью следующей команды:
$ samtools flagstat -@ 12 aln-se.bam > aln-se.flagstat
samtools flagstat - обращение к программеflagstat из набора утилит для работы с картированными чтениямиsamtools ;-@ 12 - позволение программе использовать 12 ядер для ускорения работы;aln-se.bam -bam -файл, который мы индексируем.> aln-se.flagstat - перенаправление полученных данных после анализа в отдельный файл.
Собственно, что мы получили в результате:
- Картировано чтений: 13.81%;
- Картировано чтений в корректно картированных парах чтений: 10.72%;
Выше приведённые проценты отличаются из-за того, что у нас чтения могут закартироваться на геном несколько раз, и тогда может выйти так, что некоторые чтения одной пары могут находиться далеко друг от друга и в самых разных направлениях, что и делает их "некорректно" каритрованными.
Задание 5: Получение картированных чтений
Для начала нам стоит выяснить название нашей хромосомы. Сделаю я это с помощью следующей команды:
$ samtools faidx chr3.fna
После выполнения этой команды я заглянул в получившийся файл
Следующим шагом надо получить чтения, картированные только на мою хромосому. Сначала я, собственно, выделю нужные мне чтения:
$ samtools view -h aln-se.bam NC_000003.12 > aln-se.chr3.sam
samtools view - обращение к программеview из набора утилит для работы с картированными чтениямиsamtools ;-h - указать, что в выходном файле я хочу видеть заголовок;aln-se.bam - наш файл со всеми картированными чтениями, откуда мы и будем выделять только нужные нам чтения;NC_000003.12 - этим параметром мы указываем, что в выходном файле хотим видеть чтения, картированные только на нашу хромосому;> aln-se.chr3.sam - таким образом мы перенаправляем нужные нам чтения в отдельный файл.
Как мы могли заметить, выходной файл получается в формате
$ samtools view -b aln-se.chr3.sam > aln-se.chr3.bam
samtools view - обращение к программеview из набора утилит для работы с картированными чтениямиsamtools ;-b - этим параметром мы указываем, что выходной файл хотим видеть в формате.bam ;aln-se.chr3.sam - полученный нами в прошлом пунктеsam -файл с чтениями, картированными только на нашу хромосому;> aln-se.chr3.bam - указание программе, чтобы она конвертированные данные положила в файл с указанным названием.
Вот мы и получили чтения, картированные только на нашу хромосому. Теперь посмотрим статистику:
$ samtools flagstat -@ 12 aln-se.chr3.bam > aln-se.chr3.flagstat
- Картировано чтений: 88.39%;
- Картировано чтений в корректно картированных парах чтений: 68.87%;
Как мы можем заметить, наша операция по экстрагированию нужных нам чтений сработала, так как проценты картированных и "корректно" картированных чтений увеличились: вместо 12.04% картированных чтений мы имеем 88.39%, а "корректно" картированных чтений - 68.87% вместо 10.72%.
Далее получим только "корректно" картированные пары чтений:
$ samtools view -@ 12 -f 0x2 -bS aln-se.chr3.bam > aln-se.chr3.corr.bam
И посмотрим статистику файла
$ samtools flagstat -@ 12 aln-se.chr3.corr.bam > aln-se.chr3.corr.flagstat
- Картировано чтений: 100%;
- Картировано чтений в корректно картированных парах чтений: 100%;
Собственно, у нас вместо 88.39% картированных чтений и 68.87% "корректно" картированных чтений по 100% и того, и того, что говорит нам о том, что мы получили файл, где абсолютно все чтения картировались на нашу хромосому, причём правильно, "корректно". Собственно, этого мы и добивались, это прекрасно :)
Осталось полученные прекрасные чтения проиндексировать для будущих программ:
$ samtools index -@ 12 aln-se.chr3.corr.bam
Задание 6: Подготовка выравнивания к поиску вариантов
Промаркируем дублированные чтения:
$ picard MarkDuplicates -M metrix.file.txt -I aln-se.chr3.corr.bam -O aln-se.mark.dupl.bam 2> markdupl.log
Посмотрим статистику полученного файла с промаркированными дубликатами:
$ samtools flagstat -@ 12 aln-se.mark.dupl.bam > aln-se.mark.dupl.flagstat
Как мы можем видеть с помощью