Практикум 12. Анализ данных NGS: картирование чтений, подготовка к поиску вариантов
В ходе это практикума прочтения были картированы на референсную последовательность четвёртой хромосомы человека и подготовлены к поиску вариантов.
Картирование чтений на референсный геном
Полученные в предыдущем практикуме парные триммированные чтения были картированы на референс при помощи программы bwa mem. Использованная команда приведена ниже, логи её выполнения доступны здесь.
bwa mem -t 10 chr4.fna SRR10720416_1_paired.fastq.gz SRR10720416_2_paired.fastq.gz > mapping/reads_mapped.sam 2> mapping/bwa_mem_run.log
Здесь bwa mem - команда запуска программы, -t 10 - опция, указывающая использовать под вычисления 10 потоков, chr4.fna - название файла с референсным геномом (для работы также требуются файлы, полученные при индексации этого генома), SRR10720416_1_paired.fastq.gz и SRR10720416_2_paired.fastq.gz - прочтения, которые надо картировать, > mapping/reads_mapped.sam - вывод, перенаправленный в sam-файл.
Анализ sam-файла с выравниванием
Далее приведены первые несколько строк получившегося sam-файла.
@SQ SN:NC_000004.12 LN:190214555
@PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem -t 10 chr4.fna SRR10720416_1_paired.fastq.gz SRR10720416_2_paired.fastq.gz
SRR10720416.144 77 * 0 0 * * 0 0 TGCGGAGAGTTTGAGATGTATAGGAAAGCGACGAGGATACCCCAGAGGCCTCACACATCCGCAGTCTGCTTC >((&3*%++3D@D@5AA8:?>=?A=86=;?8AE8CABACCDD5BD-'-036=:=:0&%1-=;666CAAA?@B AS:i:0 XS:i:0
SRR10720416.144 141 * 0 0 * * 0 0 NACCNNAANNNNCCTAATAAANNNTGGANTTCGGCAATCATGTATCAATANNGGTTCACAAGTTGTGAC #(*(##,0####.//66644)###),(*#(77777@C@@@C@CC@86668##1115544/58888@@@@ AS:i:0 XS:i:0
SRR10720416.146 77 * 0 0 * * 0 0 AGTCCTTCCTATACTTTGCTCCTCCAGTTCAGCCGGCAGGCTTATTGTACCATTGATAGGTTCATGAACTGCAT >'((>11;-0EBB?8A?&?ABB7B?7=7=B*CCCAGGDGGEEGGGDGD8DBB@D;D5@;D9BB@BEBEBDBB9B AS:i:0 XS:i:0
SRR10720416.146 141 * 0 0 * * 0 0 NACAAGAANNNNAAAAGAAAGNAAAGGANGTCTCCCAAAAGCTCTCTTTTGNACTTTGTTACTATCGAAG #(,'''&)####-*'11,..0#666+--#-66666C@@CC@C@C@C@CC4-#555777514@C@C@@@@@ AS:i:0 XS:i:0
SRR10720416.150 77 * 0 0 * * 0 0 GAAAATGCAAACGTCACGGAGAAGCTCACCCGTACAACATCCTGACGGCTCTGTTAACAGCCAAACAGAG ;%-&.-(-(.,01-(..*3)66>;4;5:>:)1*9->>;>>-<<=<3'%,19<>;=BB/BB@D@@@@CD@B AS:i:0 XS:i:0
SRR10720416.150 141 * 0 0 * * 0 0 NGCTCGAGANTTTCCTCAAAATAAGTTTTCTGACTTGTAACCAGAAGGAGGTTATCAGTTCAGTTTTTAAACA #(((%&(''#//--/98899@@444@@@C@44@C@@@@@CCC@C@@@@@@@@@@@CC@CCC@@C@@@@@@@C@ AS:i:0 XS:i:0
SRR10720416.151 77 * 0 0 * * 0 0 CTACCCAAAACAAGCAAACAAACAAAAAACAACACTAGATTTAGGAAGAAAACAAAGCAATGGTTCAAAATTCTG >&1(23-0-5DDDD>GG:DGDCDEDHHHAHDHH@HHHHHBFH>HDDDD>GDDDDDHFH@HDEB=DFEFFEBHBHF AS:i:0 XS:i:0
SRR10720416.151 141 * 0 0 * * 0 0 NACAAAGGGNCTACAGGCCCCAAGCAAGTGTGAAATCCTTTTTTCATAAAGGATAAGTACATGGTAAAGTTTCTG #/2--++/-#66666@@CCCCC@C@@@@@@77775<<<<<@@C44C@@@@@C@7C@@C@@@@@?@:<:<:@8@@@ AS:i:0 XS:i:0
Файл формата sam (Sequence Alignment Map) состоит из строк заголовка, начинающихся с «@», и тела файла, в котором есть 11 обязательных и произвольное количество дополнительных полей; записи содержат информацию о проанализированных чтениях. Ниже перечислены первые 12 полей, присутствующих в полученном файле.
• QNAME - имя прочтения
• FLAG - побитовая сумма флагов (характеризуют определённые свойства выравнивания)
• RNAME - имя референсной последовательности (*, если нет выравнивания)
• POS - координата в референсе, с которой начинается выравнивание (0, если выравнивания нет)
• MAPQ - качество картирования (0, если выравнивания нет)
• CIGAR - представление выравнивания как последовательности некоторого количества совпадений, несовпадений и гэпов (формат CIGAR) (*, если нет выравнивания)
• RNEXT - имя референса для следующего прочтения (*, если выравнивания нет)
• PNEXT - координата выравнивания следующего прочтения (0, если выравнивания нет)
• TLEN - длина выравнивания (0, если выравнивания нет)
• SEQ - последовательность прочтения
• QUAL - phred-score нуклеотидов в прочтении
• AS - вес выравнивания (0, если выравнивания нет)
Отмечу, что, кроме AS, в файле встречаются и другие дополнительные поля - NM, MD, MC, XC.
Данный sam-файл достаточно большой и весит 16 Гб.
Получение и индексирование bam-файла
Файл в формате sam был переведён в сортированный двоичный bam-файл с помощью команды:
samtools sort -o reads_mapped.bam reads_mapped.sam &> sort_run.log
Здесь samtools sort вызывает используемую программу, -o reads_mapped.bam задаёт имя выходного файла, reads_mapped.sam - файл на входе. Можно посмотреть лог выполнения (там почти ничего нет). Файл в формате bam весит уже 4.6 Гб (почти в 3.5 раза меньше, чем sam). Далее полученный файл был проиндексирован:
samtools index reads_mapped.bam &> index_run.log
Здесь samtools index - вызывемая программа, reads_mapped.bam - входной файл. В результате работы программы был создан файл reads_mapped.bam.bai (тоже бинарный). Лог-файл на этот раз оказался абсолютно пустым.
Анализ bam-файла
Поскольку просто прочитать бинарный файл нельзя, снова была применена программа samtools:
samtools flagstat reads_mapped.bam > flagsbam.txt 2> flagstat_run.log
Лог опять пуст. В полученном txt-файле содержалась некоторая информация из исходного bam-файла:
78023950 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 69192 + 0 supplementary 0 + 0 duplicates 12054258 + 0 mapped (15.45% : N/A) 77954758 + 0 paired in sequencing 38977379 + 0 read1 38977379 + 0 read2 9037372 + 0 properly paired (11.59% : N/A) 10392818 + 0 with itself and mate mapped 1592248 + 0 singletons (2.04% : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)
Можно видеть, что картированными на референс оказались 15.45% прочтений, при этом правильно картированы в парах были 11.59% прочтений. Вероятно, это связано с тем, что часто прочтение может картироваться на несколько участков референса, а правильно картированным чтение считается, только если оно было картировано рядом с парным ему чтением и они направлены навстречу друг другу.
Получение картированных чтений
Поскольку обрабатываемые прочтения были получены в результате полноэкзомного секвенирования, а картировались они на единственную четвёртую хромосому, далее требуется отделить правильно картированные прочтения. Сначала с целью узнать название хромосомы была использована команда:
samtools faidx chr4.fna
В полученном файле chr4.fna.fai значился, судя по всему, RefSeq AC референсной хромосомы: NC_000004.12. Для получения прочтений, картированных на хромосому, использовалась команда:
samtools view -h -bS reads_mapped.bam NC_000004.12 > reads_mapped_chr4.bam
Чтобы получить читаемую информацию из bam-файла, была выполнена команда:
samtools flagstat reads_mapped_chr4.bam > flags_mapped_chr4.txt
В выходном файле можно видеть следующее:
13646506 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 69192 + 0 supplementary 0 + 0 duplicates 12054258 + 0 mapped (88.33% : N/A) 13577314 + 0 paired in sequencing 6788657 + 0 read1 6788657 + 0 read2 9037372 + 0 properly paired (66.56% : N/A) 10392818 + 0 with itself and mate mapped 1592248 + 0 singletons (11.73% : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)
Теперь картированы 88.33% прочтений, корректно картированы в парах - 66.56%. Оба значения примерно в 5.7 раза выше, чем для картированных прочтений всего экзома (15.45% и 11.59% соответственно). Далее необходимо выделить корректно картированные прочтения; для этого была использована следующая команда:
samtools view -f 0x2 -bS reads_mapped_chr4.bam > correct_mapped_chr4.bam
Затем, как обычно, применён samtools flagstat:
samtools flagstat correct_mapped_chr4.bam > flags_correct_chr4.txt
Содержимое полученного файла:
9039503 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 2131 + 0 supplementary 0 + 0 duplicates 9039503 + 0 mapped (100.00% : N/A) 9037372 + 0 paired in sequencing 4518686 + 0 read1 4518686 + 0 read2 9037372 + 0 properly paired (100.00% : N/A) 9037372 + 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)
Можно видеть, что в файле correct_mapped_chr4.bam остались только правильно картированные прочтения (100% от всех прочтений), в отличие от reads_mapped_chr4.bam, где ещё присутствовали как некорректно картированные прочтения, так и не картированные вовсе.
Наконец, полученный bam-файл был проиндексирован:
samtools index correct_mapped_chr4.bam
Подготовка выравнивания к поиску вариантов
Для этого требовалось маркировать дублированные прочтения. Была использована команда:
picard MarkDuplicates -M metrix.file.txt -I correct_mapped_chr4.bam -O correct_markd.bam &> picard_run.log
Лог-файл доступен здесь.
Далее - flagstat:
samtools flagstat correct_markd.bam > flags_markd.txt
Содержимое flags_markd.txt:
9039503 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 2131 + 0 supplementary 495290 + 0 duplicates 9039503 + 0 mapped (100.00% : N/A) 9037372 + 0 paired in sequencing 4518686 + 0 read1 4518686 + 0 read2 9037372 + 0 properly paired (100.00% : N/A) 9037372 + 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)
По-видимому, было найдено 495 290 дублированных прочтений.