Учебный сайт Лидии Гаркуль

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

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

Для картирования чтений на хромосому была использована команда:

bwa mem -t 12 ../../bwa/chr8.fna ../../dna_seq/out_paired_SRR10720419_1.fastq.gz ../../dna_seq/out_paired_SRR10720419_2.fastq.gz > map_chr8.sam 2> log_mem

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

Файлы формата SAM состоят из заголовков и секции выравнивания. Строки заголовков начинаются символа "@". Раздел выравнивания имеет 11 обязательных столбцов. Из названия и описания находятся в таблице 1 ниже.

Таблица. 1. Первые 11 столбцов тела файла.
Номер столбца Название столбца Тип данных Описание
1 QNAME String ID чтения
2 FLAG int Характеристика выравнивания - сумма побитовых флагов
3 RNAME String Название референсной последовательности
4 POS Int Координата пары в выравнивании (1 - номер самой левой пары)
5 MAPQ Int Качество картирования
6 CIGAR String Информация о выравнивании в формате CIGAR
7 RNEXT String Название референса следующего рида в образце
8 PNEXT Int Позиция, с которой выравнивается следующий рид
9 TLEN Int Длина выравнивания
10 SEQ String Последовательность самого чтения
11 QUAL String Phred Score для нуклеотидов в чтении

Получившийся sam-файл весит около 17G.

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

С помощью команды samtools sort -o map_chr8.bam map_chr8.sam 2> log_bam переконвертикуем наш sam-файл в более легкий bam-формат. Опция -o в данной команде задает название выходного файла. Поток ошибок был перенаправлен в файл log_bam

Выходой bam-файл весит 4,8 Гб.

Далее с помощью команды samtools index -@ 12 map_chr8.bam bam-файл был проиндексирован.

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

Так как bam-файл бинарный, то для его анализа воспользуемся следующей командой: samtools flagstat -@ 12 map_chr8.bam > flagstat.

Из полученного файла flagstat известно слудующее: всего на референс картировано 11.40%; корректро картированных чтений 8.52%. Это возможно возникает из-за того что чтения могут быть картированы несколько раз.

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

a. Так как мы картируем чтения полного генома только на восьмую хромосому, у нас получились небольшие значения в предыдущем пункте. Для того чтобы посмотреть на чтения только восьмой хромосомы, найдем название нужной хромосомы следующей командной: samtools faidx ../bwa/chr8.fna -o chr8.fai. Из выходного файла узнаем название - NC_000008.11.

Зная название хромосомы, получим чтения картированные только на нее: samtools view -h map_chr8.bam NC_000008.11 > reads_chr8.sam. Далее переконвертируем bam-файл в sam-файл: samtools view -b reads_chr8.sam > reads_chr8.bam

b. Далее к полученному bam-файлу применим команду csamtools flagstat reads_chr8.bam > chr8_flagstat. Результат лежит в файле chr8_flagstat. Теперь картированно 86.98% чтений и корректно картировано 65.30%.

c. Абсолютное количество картированных чтений почти не изменилось, а относительное сильно выросло.

e. Теперь получим только правильно картированные пары чтений. Для этого воспользуемся командой: samtools view -f 0x2 -bS reads_chr8.bam > right_reads_chr8.bam

f. Далее к полученному файлу применим команду: сsamtools flagstat right_reads_chr8.bam > right_flagstat. Выходной файл лежит тут.

g. Теперь картировано 100% чтений (картированных и корректно картированных).

h. Проиндексируем наш bam-файл с корректрыми чтениями: samtools index right_reads_chr8.bam

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

Маркировка дублированных чтений с помощью команды picard MarkDuplicates -M metrix.file.txt -I right_reads_chr8.bam -O marked.bam

Далее опять применим команду flagstat к полученному файлу: samtools flagstat marked.bam > marked_flagstat.

Из файла marked_flagstat узнаем, что дублированных чтений оказалось 581752.

Текстовый файл со всеми командами тут.