Выравнивание ридов на референсный геном

Материал из Институт биоинформатики
Перейти к: навигация, поиск
GATK Best Practices

Подготовка ридов к анализу

Как уже было сказано ранее, любой анализ данных начинается с сырых ридов, хранящихся в файле с расширением .fastq. Каждый файл имеет название, содержащее номер образца, последовательность индекса (баркода), номер дорожки в проточной камере и идентификатор, обозначающий номер рида в паре (1 - прямой рид, 2 - обратный). Обычно это выглядит так: 1_TAAGGCGA-GTAAGGAG_L001_R1_001.fastq.gz . Для более удобной работы весь пайплайн построен на унифицированном наименовании файлов, содержащих риды и прочую информацию (bam-файлы, логи, индексы etc.). Все эти файлы имеют названия вида run_1.sample_1.lane_1.(R1/2)* и содержат информацию о номере запуска прибора (Illumina HiSeq 2500), номере анализируемого образца, номере дорожки и рида в паре (только для .fastq файлов). Для переименования файлов прилагается скрипт rename_fastqs.sh

#!/bin/bash

RUN=$1

for FQ in *.fastq.gz
do
    echo $FQ > name
    mv $FQ run_${RUN}.sample_${FQ%%_*}.lane_$( grep -oP 'L00\K\d' name ).R$( grep -oP '_R\K\d' name ).fastq.gz
    rm name
done
Пример графика FastQC

Контроль качества секвенирования

Любая работа с данными секвенирования не может обойтись без контроля качества, который должен производиться на всех этапах. Для контроля качества мы используем общепринятый стандарт - программу FastQC (подробнее здесь). Программа существует как в консольной, так и в GUI-версии, и может предоставлять отчет в различных форматах, включая pdf и html. Отчет FastQC содержит несколько блоков (подробно о содержании их можно почитать на официальном сайте и в документации). В обработке экзомных данных мы, в основном, обращаем внимание на распределение качества оснований по длине рида, кластеру, а также на графики GC контента и адаптерного контента (особенно важен в случае работы с китом Illumina Rapid Capture). Некоторые примеры графиков из отчета приведены ниже.

Выравнивание ридов на геном

Выравнивание на референсный геном в протоколе GATK Best Practices производится программой bwa mem. По данным статей, выравниванием при помощи этого алгоритма дает наилучшие результаты при дальнейшем определении вариантов, в том числе, и засчет наиболее аккуратной оценки качества выравнивания (подробнее о сравнении софта для выравнивания и коллинга вариантов можно почитать в этой статье. При условии форматирования названия fastq-файла по образцу выше содержимое флага -R может выглядеть, например, таким образом:

FQ=sample_X.lane_1.R1.fastq.gz
TAG=${FQ%%.R1.fastq.gz}
bwa mem -M -t 16 -R "@RG\tID:$TAG\tSM:${TAG%%.lane*}\tLIB:1\tPL:illumina" /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta ${TAG}.R1.fastq.gz ${TAG}.R2.fastq.gz | samtools view -b - > bams/${TAG}.bam

samtools merge sample_X.bam sample_X.lane_1.bam sample_X.lane_2.bam # Merging multiple lanes if exist

/Molly/barbitoff/software/samtools-1.3/samtools sort -T tmp/sample_X -o sample_X.sorted.bam sample_X.bam ; rm sample_X.bam &

Визуализация выравнивания

Для визуализации выравнивания используется геномный браузер IGV (Integrated Genomics Viewer). IGV способен визуализировать различные форматы данных, как исходные BAM файлы, так и подготовленные из них трэки, содержащие информацию о покрытии по геному (основной формат таких трэков - tdf).

Изготовить такие трэки можно при помощи IGVTools. Команда может выглядеть, например, так: /Molly/predeus/IGVTools/igvtools count -z 5 -w 50 -e 0 $TAG.dedup.bam $TAG.tdf b37 &> ../igv_logs/$TAG.maketdf.log. Оценка покрытия по геному может быть весьма важна для контроля качества, и особенно в экзомном секвенировании.