Определение вариантов (variant calling)

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

Данная статья посвящена собственно процессу получения набора вариантов в исследуемом геноме/экзоме. Для этой цели существует великое множество программ, среди которых стоит упомянуть наиболее широко используемые пакеты samtools, freebayes, а также используемый нами Genome Analysis ToolKit. Этап определения вариантов является вторым в протоколе GATK Best Practices, схема которого приведена в статье.

Определение вариантов в GATK: общие соображения

Качественное определение вариации в геноме из данных секвенирования является непростой алгоритмической и статистической задачей. Для того, чтобы увеличить чувствительность и точность, применяются разнообразные алгоритмические идеи: полногеномная сборка (реализовано в пакете fermikit за авторством Хенга Ли из Института Броада), локальная сбока (пакеты Platypus, Scalpel), а также используемая в GATK "модель уверенности в референсе" (reference condifence model) в сочетании с когортным анализом. Основная идея, лежащая за используемым в GATK алгоритмом, достаточно проста: из данных производится локальная сборка гаплотипов с использованием Де Брейн-подобного графа, каждый из гаплотипов выравнивается на геном при помощи алгоритма Смита-Уотермана, а после этого для каждого из выровненных гаплотипов вычисляется отношение правдоподобия (вероятность получить такие данные при таком гаплотипе) на основе ридов. Наиболее правдоподобный генотип вычисляется как наиболее правдоподобное сочетание гаплотипов. Для того, чтобы определить больше вариантов и увеличить их общее качество, в GATK HaplotypeCaller (GATK-HC) используется reference confidence model и когортный анализ. Для каждого из образцов определяются все сайты, в которых возможно наличие варианта исходя из данных (с соответствующим качеством и прочими параметрами). Эти сайты записываются в файл формата GVCF. После получения одиночных GVCF-файлов для всех образцов производится их совместное (когортное) генотипирование. При когортном генотипировании сильно увеличивается вероятность определить вариант при низком покрытии засчет информации о том, встречается ли такой вариант в популяции (когорте образцов). Для успешной работы алгоритма BROAD Institute рекомендует использовать когорты, состоящие из 30+ образцов. Если в исследуемой когорте недостаточно образцов, ее можно дополнить при помощи публично доступных данных (например, засчет образцов из проекта 1000 Genomes).

Определение вариантов по протоколу GATK

Для получения файлов GVCF-формата используется GATK-HC. Для работы этот walker требует BAM-файлы, полученные после стадии BQSR (смотри предыдущий этап). Запуск программы может осуществляться так:

java -Xmx4g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T HaplotypeCaller \ 
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -L /Molly/barbitoff/reference/illumina.intervals \ 
     -I sample_X.recal.bam \ # Your last-stage BAM file
     -o ../vcfs/sample_X.g.vcf 
     -ERC GVCF &> ./gatk_logs/sample_X.HaplotypeCaller.log &

HaplotypeCaller создает файлы с расширением .g.vcf, используемые в дальнейшем для генотипирования. Определение вариантов также стоит вести исключительно по интервалам экзома для увеличения среднего качества выходных вариантов. После запуска GATK-HC необходимо совместно генотипировать когорту walker-ом GenotypeGVCFs из GATK. Его можно запустить следующим образом:

java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T GenotypeGVCFs \ 
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -V gvcfs.list \ #Y our list of .g.vcf files
     -o experiment_name.raw.vcf

GenotypeGVCFs принимает на вход список директорий для анализа в файле (по одной директории в строке - оттуда будут взяты GVCF-файлы). Итоговый VCF-файл содержит неотфильтрованные результаты генотипирования, которые требуют дальнейшего процессинга при помощи как других элементов Genome Analysis ToolKit, так и при помощи других программ анализа вариантов. Процесс рекалибровки вариантов в GATK описан ниже.

Рекалибровка и фильтрация генотипов

После изначального определения вариантов VCF-файл содержит существенное количество вариантов, являющихся ложно положительными результатами (такова специфика любых данных). Для фильтрации вариантов и отбора высококачественных существует две основные стратегии: ручная фильтрация (производится для одиночных образцов или вариантов, полученных по данным РНК-сек) и рекалибровка, основанная на модели (variant quality score recalibration - VQSR). Второй подход наиболее широко используется в GATK и является протокольным стандартом. Для построения модели фильтрации используются базы данных вариантов: dbSNP, HapMap, Omni, Mills. После построения позитивной модели (модель хорошего варианта) и негативной модели (модель плохого варианта) эти модели применяются к набору вариантов в эксперименте и для каждого варианта посчитывается VQSLOD-score - логарифм отношения вероятностей отнесения варианта к позитивной и негативной модели. На основании VQSLOD-значений и производится конечная фильтрация.

Рекалибровка вариантов производится в два последовательных запуска - сначала для SNP, а затем для indel-ов. Для построения модели используется walker VariantRecalibrator, а для применения модели - ApplyRecalibration. Запустить VariantRecalibrator можно, например, таким способом:

java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T VariantRecalibrator \
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -input experiment_name.raw.vcf \ # Your raw VCF file
     -resource:hapmap,known=false,training=true,truth=true,prior=15.0 \
              /Molly/barbitoff/gatk-bundle-b37/hapmap_3.3.b37.vcf \
     -resource:omni,known=false,training=true,truth=true,prior=12.0 \
              /Molly/barbitoff/gatk-bundle-b37/1000G_omni2.5.b37.vcf \
     -resource:mills,known=false,training=true,truth=true,prior=12.0 \
              /Molly/barbitoff/gatk-bundle-b37/Mills_and_1000G_gold_standard.indels.b37.vcf \
     -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 \
              /Molly/barbitoff/gatk-bundle-b37/dbsnp_138.b37.vcf \
     -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an InbreedingCoeff \
     -mode SNP \ # SNP or INDEL for respective types of variation
     -recalFile experiment_name.recal \ # File with model
     -tranchesFile experiment_name.tranches \ # File with tranch characteristics
     -rscriptFile experiment_name.plots.R # File with visualisation script

Для каждого из режимов на вход подается информация об используемых базах для построения модели (включены в GATK bundle). На выходе записывается файл с параметрами рекалибровки и файл с "траншами", содержащий информацию о различных "срезах" исходных данных, а также скрипт для визуализации этой информации. Основная идея состоит в том, что каждый срез или транш представляет собой кусок исходного набора вариантов, в той или иной степени подходящий под построенную модель. Для каждого среза характерен определенный процент истинно и ложно положительных вариантов. Файл с характеристиками траншей имеет примерно такую структуру:

# Variant quality score tranches file
# Version number 5
targetTruthSensitivity,numKnown,numNovel,knownTiTv,novelTiTv,minVQSLod,filterName,model,accessibleTruthSites,
callsAtTruthSites,truthSensitivity
90.00,82241,5446,2.9206,1.6414,2.9449,VQSRTrancheSNP0.00to90.00,SNP,60228,54205,0.9000
99.00,94865,6139,2.8446,1.6498,-0.5785,VQSRTrancheSNP90.00to99.00,SNP,60228,59625,0.9900
99.90,111293,10947,2.7034,1.5998,-11.9256,VQSRTrancheSNP99.00to99.90,SNP,60228,60167,0.9990
100.00,113855,13451,2.6707,1.4327,-39931.7923,VQSRTrancheSNP99.90to100.00,SNP,60228,60228,1.0000

Первая колонка задает уровень чувствительности, колонки со второй по четвертую - характеристики набора данных, отфильтрованного по этому уровню, пятая - уровень упомянутого выше VQSLOD-порога. Для лучшего понимания происходящего можно воспользоваться скриптом для визуализации. На выходе будет примерно такой результат:

Tranches.png

Каждый столбец соответствует определенному срезу данных, удовлетворяющему порогу чувствительности из транш-файла выше. Окрашенные синим области соответствуют реальным вариантам, оранжевые - ложно положительным. Соответственно, необходимый порог фильтрации выбирается исходя из требуемого отношения TP:FP экспериментатором. Срез, соответствующий этому порогу, и будет являться итоговым набором вариантов. !!! В нынешней версии программы визуализация не соответствует реальной структуре данных, то есть фильтрацию стоит производить по транш-файлу, выбирая подходящий минимальный VQSLOD и количество вариантов Для применения модели можно запустить GATK следующим образом:

java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T ApplyRecalibration \
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -input experiment_name.raw.vcf \
     -mode SNP \ # SNP or INDEL
     -recalFile experiment_name.recal \
     -tranchesFile experiment_name.tranches \
     -o experiment_name.recal.vcf \
     -ts_filter_level 90.0 # Your chosen TScore for filtration

Важно, что для построения модели и применения оной при втором прогоне нужно использовать файлы, полученные на выходе после запуска в первом режиме!

После рекалибровки вариантов отфильтрованные можно исключить из VCF-файла примерно такой командой:

 
java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T SelectVariants \
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -V experiment_name.recal.vcf \
     -select 'vc.isNotFiltered()' \
     -o experiment_name.out.vcf

После окончания рекалибровки INDEL-ов скрипт автоматически отфильтрует варианты, не прошедшие рекалибровку, при помощи программы SelectVariants из GATK, а также отсортирует файлы по директориям.

Перекалибровка генотипов и фильтрация по GQ

По окончанию рекалибровки вариантов рекомендуется также произвести уточнение качества генотипирования образцов при помощи "априорных вероятностей", полученных с использованием родословной в формате PED или по базам данных вариации. Эта стадия позволяет более точно оценивать индивидуальное (не когортное, как в поле QUAL!) качество генотипа и производить фильтрацию низкокачественных генотипов присвоением им флага 'lowGQ'. Описанные выше процедуры можно провести так:

java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T CalculateGenotypePosteriors \
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \
     -V experiment_name.out.vcf \
     --supporting /Molly/barbitoff/gatk-bundle-b37/1000G_phase3_v4_20130502.sites.vcf \
     -o experiment_name.GR.vcf

java -Xmx64g -jar /Molly/barbitoff/software/gatk-protected/target/GenomeAnalysisTK.jar \
     -T VariantFiltration \
     -R /Molly/barbitoff/reference/GATK_b37/human_g1k_v37.fasta \ 
     -V experiment_name.GR.vcf \
     -G_filter "GQ < 20.0" -G_filterName lowGQ \
     -o experiment_name.GF.vcf

По окончанию этого этапа на выходе остается файл с высококачественными вариантами, которые дальше могут подвергаться дополнительной аннотации, фильтрации и анализу. Об этом будет рассказано в следующих статьях.