Введение

в анализ данных NGS

Блок 3 - Занятие 12

Анастасия Жарикова

В прошлых сериях

  • чтения (fastqc)

  • проверяем качество

  • исправляем при необходимости

  • убеждаемся, что качество исправлено


  • готовим референс

  • картируем чтения на референсный геном

  • получаем sam (bam)

Что делать дальше?

Дано:

  • «очищенные» чтения хорошего качества (fastq)

  • Последовательность референсного генома (fasta)

Задача:

Каждому чтению найти свое место на геноме - картирование

Картирование

Программы:

  • bowtie

  • bwa

  • hisat2

Есть много других!

Шаг 0. Подготовка референса: индексирование Для каждой программы свой индекс!

Шаг следующий – картирование чтений на референс

Получаем .sam или .bam

SAM/BAM

Некоторые файлы содержат

  • “шапку”

  • основное содержание

SAM

“Шапка”

@HD     VN:1.5  GO:none SO:coordinate
@SQ     SN:chr1 LN:248956422
@SQ     SN:chr2 LN:242193529
@SQ     SN:chr3 LN:198295559
@SQ     SN:chr4 LN:190214555
@SQ     SN:chr5 LN:181538259
@SQ     SN:chr6 LN:170805979
@SQ     SN:chr7 LN:159345973
@SQ     SN:chr8 LN:145138636
@SQ     SN:chr9 LN:138394717
@SQ     SN:chr10        LN:133797422
@SQ     SN:chr11        LN:135086622
@SQ     SN:chr12        LN:133275309
@SQ     SN:chr13        LN:114364328
@SQ     SN:chr14        LN:107043718
@SQ     SN:chr15        LN:101991189
@SQ     SN:chr16        LN:90338345
@SQ     SN:chr17        LN:83257441

SAM

Содержание

NB551509:3:HGW7GBGX9:3:11609:7745:2302  99      chr1    10012   0       35M16S  =       10176   209     CTAACCCTAACCCTAACCCTAACCCTAACCCTAACGCTATCGGTACCGCTA     AAAAAE/EEEEEEEE6EEEEEEEA//E/A/A<////////</////////<     MC:Z:25S45M     MD:Z:35 PG:Z:MarkDuplicates     RG:Z:NextSeq5500        NM:i:0  MQ:i:0  AS:i:35 XS:i:39
NB551509:3:HGW7GBGX9:3:11609:7745:2302  147     chr1    10176   0       25S45M  =       10012   -209    AACCCTAACACCCACCCACACCCTCAACCTAACGCTCACCCTAACCCTAACCCTAACCCTAACCCTAACC  A///</////E////////E//E/////A//AE/E//EEEEEEE/EEEE/EEEAEEEEEEEEEEEAAAAA  MC:Z:35M16S     MD:Z:8C2A33     PG:Z:MarkDuplicates     RG:Z:NextSeq5500        NM:i:2  MQ:i:0  AS:i:35 XS:i:38
NB551509:3:HGW7GBGX9:1:12102:24585:13487        99      chr1    10817   0       75M     =       11013   220
     GGGGTGGAGGCGTGGCGCAGGCGCAGAGAGGCGCGCCGCGCCGGCGCAGGCGCAGAGACCCCTGCTACCGCGGCC     AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/E//EAE<A//EE/E//</6/E///<E/A/AAEEE//<     XA:Z:chr15,-101980076,75M,3;chr1,+10893,75M,5;chr15,-101980000,75M,5;chr12,+11022,75M,5;        MC:Z:51S24M     MD:Z:59A1A10T2  PG:Z:MarkDuplicates     RG:Z:NextSeq5500        NM:i:3  MQ:i:0  AS:i:62 XS:i:62
NB551509:3:HGW7GBGX9:1:12102:24585:13487        147     chr1    11013   0       51S24M  =       10817   -220
    GGGCCTTGGGGGGGTGGGGGGGGGCTGGCGCTGACCCGAACGCCTCGGGGGGGGGGTGGGGGGGGCGTGTGTTGC     <//////EEEA/6//E</6/EE////E///////////////A////EEE/EEEAA//EEAEEEEEE/EEAAAAA     XA:Z:chr1,-193764955,42S25M8S,0;chr16,+12525686,6S23M46S,0;chr15,-79311000,49S21M5S,0;chr15,+101979931,24M51S,1;        MC:Z:75M        MD:Z:6T17       PG:Z:MarkDuplicates     RG:Z:NextSeq5500        NM:i:1  MQ:i:0  AS:i:19 XS:i:25
NB551509:3:HGW7GBGX9:2:11111:18968:7557 163     chr1    12927   0       74M     =       13001   149     AGCTACAAGCAGCAAACAGTCTGCATGGGTCATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCAGGGGCCCCCA      AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEE      MC:Z:75M        MD:Z:74 PG:Z:MarkDuplicates     RG:Z:NextSeq5500

SAM

  • Read Name: NB551509:3:HGW7GBGX9:3:11504:22127:11632

  • SAM flag: 99

  • contig name or * for unmapped: chr1

  • mapped position of base 1 of a read on the reference sequence: 15678

  • MAPQ mapping quality: 0

  • CIGAR string describing insertions and deletions: 75M

  • Name of mate: =

  • Position of mate: 15726

  • Template length: 123

  • Read Sequence: CAGGGGCCAGCTGGCAAGAGCAGGGGGTGGGCAGAAAGCACCCGGTGGACTCAGGGCTGGAGGGGAGGAGGCGAT

  • Read Quality: AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE

  • Additional information in TAG:TYPE:VALUE format: MC:Z:75M MD:Z:75 PG:Z:MarkDuplicates RG:Z:NextSeq5500 NM:i:0 MQ:i:0 AS:i:75 XS:i:75

SAM: flag

flags

One of the reads is unmapped: 73, 133, 89, 121, 165, 181, 101, 117, 153, 185, 69, 137

Both reads are unmapped: 77, 141

Mapped within the insert size and in correct orientation: 99, 147, 83, 163

Mapped within the insert size but in wrong orientation: 67, 131, 115, 179

Mapped uniquely, but with wrong insert size: 81, 161, 97, 145, 65, 129, 113, 177

SAM: flag

explain

SAM: CIGAR

  • M: alignment match (can be a sequence match or mismatch) yes yes

  • I: insertion to the reference yes no

  • D: deletion from the reference no yes

  • N: skipped region from the reference no yes

SAM: tags

explain раз

explain два

explain три

VCF

Спецификация

VCF: структура

Три основные части:

  • Шапка (строки начинаются с ##)

  • Строка (одна) с заголовками столбцов (начинается с #)

  • Информация о вариантах

VCF: столбцы

manual раз

manual два

8 фиксированных колонок

  • CHROM – имя хромосомы

  • POS – позиция варианта

  • ID – может быть любая информация о варианте, но обычно пустой (.)

  • REF – референсная аллель

  • ALT – альтернативная аллель

  • QUAL – качество варианта (Phred-scaled)

  • FILTER – PASS (если ранее была осуществлена маркировка по каким-либо показателям: покрытие, качество и пр)

  • INFO – различные характеристики варианта

  • FORMAT – список параметров варианта для конкретного образца

  • HG00096 – значения параметров, указанных в столбце FORMAT для конкретного образца (в заголовке – ID образца)

VCF: метрики образца

  • Колонка FORMAT: GT:AD:DP:GQ:PL

  • Колонка ID образца: 0/1:21,4:25:99:1220,108,0

VCF: GT

Кодирует генотип варианта

Для диплоидных организмов:

  • 0 – референсный аллель

  • 1 – альтернативный аллель

Образец по варианту:

  • 0/0 – референсная гомозигота

  • 0/1 – гетерозигота

  • 1/1 – альтернативная гомозигота

VCF: AD, DP

Отражает покрытие

  • AD – количество чтений, которые поддерживают каждую из возможных аллелей; участвуют все чтения, использованные при поиске вариантов

  • DP – общее количество чтений, прошедших фильтрацию и поддерживающих каждую из представленных аллелей

VCF: PL, GQ

Отражает качество генотипа

  • PL – нормализованные «вероятности» возможных генотипов (по шкале Phred). Поле содержит 3 числа, что соответствует генотипам 0/0, 0/1, 1/1. PL наиболее вероятного генотипа = 0

  • GQ – вычисляется на основании PL, представляет собой разницу «вероятностей» двух наиболее вероятных генотипов (но не более 99). Низкие значения (т.е. << 99) – в генотипе нет уверенности

Multiple VCF

В одном VCF файле может быть представлена информация сразу о нескольких образцах

В конце будут добавлены столбцы на каждый образец

QUAL – максимальный из возможных

BED

спецификация

BedGraph

смотрим самостоятельно

BEDtools

manual раз

manual два

Очень хороший инструмент для работы с геномными интервалами

Более 35 опций + параметры

Система координат

Обозначить координаты интервала между 3 и 7 нуклеотидов

  • 0-based: [2,7) - BAM, BCFv2, BED

  • 1-based: [3,7] - SAM, VCF, GFF, Wiggle formats

Отладка кода

Некоторые программы отрабатывают довольно долго

Не всегда получается сразу написать корректно работающий код

При отладке кода лучше использовать небольшой тестовый набор данных

Например, взять 1000 пар чтений

tmux

Терминальный мультиплексор

Позволяет работать с несколькими сессиями в одном окне командной строки

Позволяет подключаться (и отключаться) к текущему состоянию сессии

Запущенные программы и процессы продолжают работать

tmux

  • tmux - без параметров создает сессию 0

  • tmux new -s myname - создает сессию с именем myname

  • tmux ls - проверить наличие сессий в tmux

  • tmux a -t myname - подключиться к сессии myname

  • tmux kill-session -t myname - закрыть сессию myname

tmux

Команды начинаются с префикса Ctrl+B

Сначала нажимаем Ctrl+B, а затем нужную команду:

  • с - новое окно

  • n - следующее окно

  • p - предыдущее окно

  • ” - деление окна горизонтально

  • % - деление окна вертикально

  • x - закрыть окно

  • d - отключиться от сессии

tmux

мануал раз

мануал два

Конец!