Учебный сайт
Владимира Ноздрина

В манифесте анаконды не укажут имена
Тех, кто видел свет жиронды из окна.
Кобыла и трупоглазые жабы, "Будущее вечно"

Поиск полиморфизмов у человека.

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

Картировать будем триммированные чтения на проиндексированный референсный геном десятой хромосомы человека. Всё это было получено в предыдущем практикуме.
Итак, картирование производилось следующей командой: $ bwa mem -t 10 bwa/chr10.fna dna_reads/trim_p1.fastq.gz dna_reads/trim_p2.fastq.gz > map/out.sam 2> map/log.txt Что всё это значит:
bwa mem – вызов программы-картировщика.
-t 10 – использовать 10 ядер процессора, чтобы считалось быстрее.
bwa/chr10.fna – последовательность референса.
dna_reads/trim_p1.fastq.gz dna_reads/trim_p2.fastq.gz – входные файлы.
> map/out.sam – записать вывод в файл map/out.sam.
2> map/log.txt – записать логи в файл map/log.txt.

Анализ sam файла с выравниванием

Полученный файл весит 17 гигабайт. Он состоит из заголовка, выравниваний и информации о выравниванияо. На одно выравнивание выделяется ровно одна строка. Строки заголовка начинаются с @. В нём содержится общая информация о том, как было получен этот файл. Ниже можно посмотреть на его первые строки. 
$ head out.sam 
@SQ SN:NC_000010.11 LN:133797422
@PG ID:bwa  PN:bwa  VN:0.7.17-r1188 CL:bwa mem -t 10 bwa/chr10.fna dna_reads/trim_p1.fastq.gz dna_reads/trim_p2.fastq.gz
SRR10720412.2 77  * 0 0 * * 0 0 GCTGAGTGATGGAGGGGGATAAGTCTTCTCGATAATAAACCTTCCCCAGGCCATCGGTCGTGTCCAAGTCCGCCA GGGD@F>EEBDDFFDDEEDEGGEEGGBGGFGEGBBGDGDDGGGGGGEDGGGEF?EGIGGEGIIGIEGBEGEEEEEGFFFFICC#CCB?=AB AS:i:0  XS:i:0
SRR10720412.3 141 * 0 0 * * 0 0 TGTGTTCCTTTTGCTTTGTNGTGTCTCCTGAGCAACCGCCNNTNNNNACAGTAAACATATTTTCTTTAATTATTA GEGEGIGIHHDIHIGFFFB#C=CAAIIIGIBHIIHEFCEF##?####498A=?AAIIIIIE774#####8>##7889#4##<<7><#######0<>::@<=GEGIFIGHFBDGGDDIHIHI AS:i:0  XS:i:0
SRR10720412.6 77  * 0 0 * * 0 0 CATTTTTATTGCATTTTTGCAGTTNNTNCCATTTATTGAAAATNACAATTGTAACAATTTTCTTA HHHHHHEHHHHHDHHHHHHH>@???BBGHHEHEFB#EDBBBBHEHHHHFHHHHHHHH AS:i:0  XS:i:0
SRR10720412.6 141 * 0 0 * * 0 0 ATGGTTTCTTTACAGAANNNNNNTCNNCACCNTNNNTTCGNNNNNNNNNCTGTTCGAGGTGGAGTGACAGGACGT GGGGFIIIIIIIIHI88######59##?<:<#<###<8;<#########657559DDGG@GFED@DDEBE@;FA< AS:i:0  XS:i:0
Вот, что обозначают поля:
  1. QNAME – Название чтения, которое картируется.
  2. FLAG – Краткая информация о картировании чтения. (В двоичном формате это число имеет смысл...)
  3. RNAME – Название референсной последовательности. (Последовательности референса, к которой чтение картировалось)
  4. POS – Позиция нуклеотида референса, с которого начинается выравнивание.
  5. MAPQ – Качество картирования. Вычисляется по формуле -10lgP, где P – вероятность того, что чтение картировалось неверно.
  6. CIGAR – Сжатая информация о выравнивании в формате CIGAR.
  7. RNEXT – Название референсной последовательности для следующего чтения.
  8. PNEXT – POS для первого выравнивания следующего чтения.
  9. TLEN – Длина покрытия.
  10. SEQ – Последовательность чтения.
  11. QUAL – Понуклеотидное качество чтения.
  12. AS – Вес выравнивания. Опциональное поле.

Получение и индексирование bam файла

bam файл получим с помощью следующей команды: $ samtools sort -o out.bam out.sam samtools sort – вызов программы, которая сортирует выравнивания в sam файле. По умолчанию по увеличению левой координаты.
-o out.bam – вывод программы записать в файл out.bam.
out.sam – исходный файл.
Полученный файл весит 4.94 гигабайт. Теперь его нужно проиндексировать. Индексировать будем следующей командой: $ samtools index out.bam samtools index – вызов программы-индексировщика.
out.bam – исходный файл.

Анализ bam файла

Посмотрим, сколько правильных картирований получилось. Для этого используем следующую команду: $ samtools flagstat out.bam > flagstat.txt Вывод записался в файл flagstat.txt, но я могу вставить всё его содержимое прямо сюда: 
$ cat flagstat.txt 
79756299 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
49945 + 0 supplementary
0 + 0 duplicates
10002438 + 0 mapped (12.54% : N/A)
79706354 + 0 paired in sequencing
39853177 + 0 read1
39853177 + 0 read2
7625834 + 0 properly paired (9.57% : N/A)
8526868 + 0 with itself and mate mapped
1425625 + 0 singletons (1.79% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
Видно, что хоть как-то картировалось 10 002 438 чтений (12.54% от общего числа), а правильно парно картировалось только 7 625 834 (9.57%). Эти числа отличаются, потому что правильным считается картирование, при котором чтения направлены друг к другу, и располагаются недалеко. Но может случиться так, что скартируется только одно из двух чтений. Или, например, у какого-то из чтений будет альтернативное место, куда можно картироваться. В общем, можно придумать много ситуаций, когда эти числа будут отличаться.

Получение картированных чтений

Чтобы получить чтения, картированные на десятую хромосому человека(NC_000010.11), воспользуемся следующей командой: $ samtools view -h out.bam NC_000010.11 | samtools sort -o view_chr10.bam Получился сразу bam файл с картированными чтениями. Снова применим flagstat и посмотрим на выдачу: 
$ samtools flagstat view_chr10.bam > flagstat2.txt
$ cat flagstat2.txt 
11428063 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
49945 + 0 supplementary
0 + 0 duplicates
10002438 + 0 mapped (87.53% : N/A)
11378118 + 0 paired in sequencing
5689059 + 0 read1
5689059 + 0 read2
7625834 + 0 properly paired (67.02% : N/A)
8526868 + 0 with itself and mate mapped
1425625 + 0 singletons (12.53% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
Итого, картировалось 10002438 чтения (87.53% от всех оставшихся в этом файле чтений), из них правильно парно картировалось 7625834 (67.02%). Абсолютные числа не изменились по сравнению с предыдущими, хотя доли картированных чтений от общего числа чтений в файле увеличилась. По какой-то причине существуют некартированные чтения, которым почему-то причислена хромосома. У некоторых из таких стоит флаг 117, в одном из пунктов которого написано, что последовательность, комплементарная паре этого чтения, картировалась, что может объяснять, почему чтению дали хромосому. Но есть и флаг 69, который просто означает, что ничего никуда не картировалось.
Теперь получим только правильно парно картированные чтения. Для этого применим следующую команду: $ samtools view -f 0x2 -bS view_chr10.bam | samtools sort -o paired_chr10.bam Мы получили новый файл. Ещё раз посмотрим на flagstat: 
$ samtools flagstat paired_chr10.bam > flagstat3.txt
$ cat flagstat3.txt 
7627753 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
1919 + 0 supplementary
0 + 0 duplicates
7627753 + 0 mapped (100.00% : N/A)
7625834 + 0 paired in sequencing
3812917 + 0 read1
3812917 + 0 read2
7625834 + 0 properly paired (100.00% : N/A)
7625834 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
Всего в этом файле 7627753 картированных чтения (100%), а правильно парно картированных 7625834 (100%). Я понятия не имею, почему эти числа различаются, но оба 100%. Разница между ними 1919, что совпадает с числом под пунктом supplementary, но я не знаю, что это такое. В любом случае, число правильно картированных из этого файла совпадает с этим же числом из предыдущего файла.
Проиндексируем самый последний полученный файл: $ samtools index paired_chr10.bam

Подготовка выравнивания к поиску вариантов

Для того, чтобы маркировать дублированные чтения применим следующую команду: $ picard MarkDuplicates -M metrix.file.txt -I paired_chr10.bam -O mark_chr10.bam Использовав flagstat ещё раз, получим, что дублированных чтений 474484, но ещё раз выдачу я приводить не буду.