Бенчмаркинг детекции геномных полиморфизмов

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

Специфичность обычно обозначает долю результатов, идентифицированных как отрицательные среди истинно отрицательных событий. Определение всевозможных "отрицательных" событий в случае с геномными вариантами становится расплывчатым. В связи с этим можно (1) охарактеризовать долю ложноположительных вариантов среди всех вариантов, т.е. подсчитать FDR, либо, (2) если рассматриваются одни и те же геномные регионы, сравнить количество ложноположительных результатов в абсолютном выражении.

Если для тестов используются результаты секвенирования некоторой панели генов или полноэкзомного секвенирования, то важно рассматривать варианты внутри только этих регионов, т.к. остальные районы генома не будут иметь покрытия ридами.

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

Команды, используемые для генерации сравнений, приведены ниже.

Тестирование герминального пайплайна

Тестовые образцы

Новейшим "золотым стандартом" для сложных вариантов и регионов стал образец HG002 (также известный как GM24385 и huAA53E0). Клеточная линия HG002 была секвенирована, собрана и улучшена консорциумом "Telomere-to-Telomere", совместно с консорциумом "Human Pangenome Reference" и консорциумом "Genome in a Bottle".

В 2016 году прошло соревнование PrecisionFDA Truth Challenge, в рамках которого участникам была предоставлена возможность протестировать свои пайплайны NGS на неохарактеризованном ранее образце HG002. Мы сравнили результаты нашего герминального пайплайна с опубликованной таблицей результатов соревнования.

При подготовке данного раздела была использована версия v1.1 сборки T2T-HG002.

Тестовые образцы

Известным источником "золотого стандарта" для поиска геномных вариантов является консорциум "Genome in a Bottle". Этот консорциум создан на базе Национального Интитута Стандартов и Технологий США и, в частности, занимается характеризацией человеческих генов для создания стандартов для валидации.

В рамках этого консорциума подробно охарактеризованы образец NA12878/HG001 из проекта HapMap и два семейных трио: трио евреев Ашкенази и трио китайцев Хань.

При подготовке данного раздела была использована версия v4.2.1 датасета "Genome In a Bottle". В качестве сырых данных использовались продукты полноэкзомного секвенирования панелью Nextera Rapid Capture.

Результаты сравнения

Для образца HG002 получила результаты бенчмаркинга, как это было реализовано в соревнования от GIAB для этого образца (https://precision.fda.gov/challenges/truth/results).

Свела данные соревнований с нашими в одну таблицу:

Метка	Отправитель	Организация	SNP-Fscore (%)	SNP-recall (%)	SNP-precision (%)	INDEL-Fscore (%)	INDEL-recall (%)	INDEL-precision (%)
-	-	NOVEL	98.8578	98.4426	99.2765	97.2653	95.3087	99.3039
anovak-vg	Adam Novak et al.	vgteam	98.4545	98.3357	98.5736	70.4960	69.7491	71.2591
astatham-gatk	Aaron Statham et al.	KCCG	99.5934	99.2091	99.9807 🏆	99.3424 ⭐	99.2404 ⭐	99.4446 ⭐
asubramanian-gatk	Ayshwarya Subramanian et al.	Broad Institute	98.9379	97.9985	99.8954	98.8418	98.5404	99.1451
bgallagher-sentieon	Brendan Gallagher et al.	Sentieon	99.9296 ⭐	99.9673 🏆	99.8919	99.2678	99.2143 ⭐	99.3213
cchapple-custom	Charles Chapple et al.	Saphetor	99.8448	99.8832	99.8063	99.1388	98.8448	99.4346 ⭐
ccogle-snppet*	Christopher Cogle et al.	CancerPOP
ciseli-custom	Christian Iseli et al.	SIB	97.7648	98.8356	96.7169	83.5453	82.5314	84.5844
ckim-dragen	Changhoon Kim	Macrogen	99.8268	99.9524 ⭐	99.7015	99.1359	99.1574 ⭐	99.1143
ckim-gatk	Changhoon Kim	Macrogen	99.6466	99.4788	99.8150	99.2271	99.1551 ⭐	99.2992
ckim-isaac	Changhoon Kim	Macrogen	98.5357	97.1616	99.9494 ⭐	95.8099	93.7006	98.0163
ckim-vqsr	Changhoon Kim	Macrogen	99.2866	98.6511	99.9303 ⭐	99.2541	99.0614	99.4476 ⭐
dgrover-gatk	Deepak Grover	Sanofi-Genzyme	99.9456 ⭐	99.9631 ⭐	99.9282 ⭐	99.4009 🏆	99.3458 🏆	99.4561 ⭐
egarrison-hhga	Erik Garrison et al.	-	99.8985	99.8365	99.9607 ⭐	97.4253	97.1646	97.6874
eyeh-varpipe	ErhChan Yeh et al.	Academia Sinica	99.4670	99.9638 ⭐	98.9751	92.5779	91.3854	93.8021
gduggal-bwafb	Geet Duggal et al.	DNAnexus Science	99.7820	99.8619	99.7021	96.9474	95.5004	98.4390
gduggal-bwaplat	Geet Duggal et al.	DNAnexus Science	98.8646	98.0471	99.6958	92.6621	87.0843	99.0034
gduggal-bwavard	Geet Duggal et al.	DNAnexus Science	99.3249	99.0431	99.6083	87.3464	87.1769	87.5166
gduggal-snapfb	Geet Duggal et al.	DNAnexus Science	99.2501	99.8026	98.7037	92.2602	90.5733	94.0112
gduggal-snapplat	Geet Duggal et al.	DNAnexus Science	99.0030	98.6815	99.3266	76.4210	69.0418	85.5664
gduggal-snapvard	Geet Duggal et al.	DNAnexus Science	99.0871	98.9341	99.2406	83.0264	83.4429	82.6139
ghariani-varprowl	Gunjan Hariani et al.	Quintiles	99.3496	99.8685	98.8361	87.2025	87.3272	87.0781
hfeng-pmm1	Hanying Feng et al.	Sentieon	99.9496 ⭐	99.9227 ⭐	99.9766 ⭐	99.3397 ⭐	99.0289	99.6526 ⭐
hfeng-pmm2	Hanying Feng et al.	Sentieon	99.9416 ⭐	99.9254 ⭐	99.9579 ⭐	99.3119 ⭐	99.0152	99.6103 ⭐
hfeng-pmm3	Hanying Feng et al.	Sentieon	99.9548 ⭐	99.9339 ⭐	99.9756 ⭐	99.3628 ⭐	99.0161	99.7120 🏆
jlack-gatk	Justin Lack	NIH	99.7200	99.9393 ⭐	99.5016	98.6899	98.8138	98.5664
jli-custom	Jian Li et al.	Roche	99.9382 ⭐	99.9603 ⭐	99.9160	99.3675 ⭐	99.0788	99.6580 ⭐
jmaeng-gatk	Ju Heon Maeng	Yonsei University	99.6144	99.4608	99.7686	99.1098	99.0216	99.1981
jpowers-varprowl	Jason Powers et al.	Q2 Solutions	99.5004	99.5447	99.4561	86.4885	85.2886	87.7226
ltrigg-rtg1	Len Trigg	RTG	99.8754	99.8921	99.8587	99.0160	98.3355	99.7061 ⭐
ltrigg-rtg2	Len Trigg	RTG	99.8749	99.8935	99.8562	99.2539	98.8759	99.6347 ⭐
mlin-fermikit	Mike Lin	DNAnexus Science	98.8629	98.2311	99.5029	95.5997	94.8918	96.3183
ndellapenna-hhga	Nicolas Della Penna	ANU	99.8818	99.8118	99.9519 ⭐	97.3838	97.0938	97.6756
qzeng-custom	Qian Zeng	LabCorp	99.4966	99.2413	99.7533	96.8316	96.8703	96.7929
raldana-dualsentieon	Rafael Aldana et al.	Sentieon	99.9260 ⭐	99.9131 ⭐	99.9389 ⭐	99.1095	98.7566	99.4648 ⭐
rpoplin-dv42	Ryan Poplin et al.	Verily Life Sciences	99.9587 🏆	99.9447 ⭐	99.9728 ⭐	98.9802	98.7882	99.1728

Результаты сравнения

Была проверена чувствительность методики детекции вариантов. Оказывается, что если рассмотреть все варианты, попадающие в заявленную панель Nextera и в регионы уверенности (т.е. те, в которые консорциум гарантирует верные коллы), то чувствительность оказывается 93.23%.

Рассмотрение незадетектированных вариантов позволяет заметить, что одной из основных причин является отсутствие ридов, покрывающих этот вариант. Пример такого варианта представлен на рис. 1:

Рисунок 1. Скриншот программы IGV, показывающий пример варианта, незадетектированного по причине отсутствия покрывающих его ридов в выравнивании (обведен красным кругом). Самые верхние трэки представлены вариантами "золотого стандарта" и задетектированными вариантами. Ниже расположены выравнивания, и в самом низу есть трэк регионов, которые заявлены как покрываемые используемой панелью Nextera.

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

Рисунок 2. Скриншот программы IGV, на котором показан пример варианта, незадетектированного по причине расположения в "сложном" регионе (обведен красным кругом). Самые верхние трэки представлены вариантами "золотого стандарта" и задетектированными вариантами. Ниже расположены выравнивания. В самом нижнем блоке есть трэк регионов, которые заявлены как покрываемые используемой панелью Nextera, а также трэк регионов, покрытых хотя бы двумя ридами.

В связи с этими эффектами при подсчете чувствительности можно сузить область рассмотрения до тех регионов, где есть, скажем, хотя бы 2 рида. В таком случае чувствительность возрастает до 97,65% (по сравнению с 93,23%). Или, например, если сузить область до регионов, не входящих в "сложные", чувствительность возрастает до 95,16%. Наконец, если рассмотреть только регионы, в которых имеется ненулевое покрытие (хотя бы 2 рида) и которые не относятся к "сложным" (т.е. пересечение этих двух областей), то чувствительность составляет 99.22%.

Рисунок 3. Столбчатая диаграмма чувствительности герминативного пайплайна. Показана чувствительность (в %) задетектированных вариантов "золотого стандарта" с разбивкой на варианты, (1) попадающие в заявленные регионы панели ("Все варианты"), (2) попадающие в регионы панели, но не в "сложные" районы, (3) попадающие в регионы панели и при этом покрытые хотя бы 2 ридами в выравнивании, (4) попадающие в пересечение всех этих областей. Во всех случаях рассматриваются только регионы высокой надежности, заявленные консорциумом GIAB.

Также важно знать, какое количество ложных вариантов продуцирует пайплайн. Для оценки можно, например, использовать FDR (False discover rate), т.е. долю ложноположительных вариантов среди всех обнаруженных вариантов. Значения для вышеупомянутых регионов можно посмотреть на рис. 4:

Рисунок 4. Столбчатая диаграмма FDR герминативного пайплайна. Показана FDR (в %) задетектированных вариантов "золотого стандарта" с разбивкой на варианты, (1) попадающие в заявленные регионы панели ("Все варианты"), (2) попадающие в регионы панели, но не в "сложные" районы, (3) попадающие в регионы панели и при этом покрытые хотя бы 2 ридами в выравнивании, (4) попадающие в пересечение всех этих областей. Во всех случаях рассматриваются только регионы высокой надежности, заявленные консорциумом GIAB.

В первую очередь бросается в глаза огромная разница в количестве ложноположительных вариантов между всеми участками, покрытыми панелью, (колонка "Все варианты" на рис. 4) и теми, которые не относятся к "сложным" регионам. Если во всех заявленных регионах панели FDR составляет 5,46% (т.е. более 5% всех обнаруживаемых в этих регионах вариантов - ложноположительные), то с отбрасываением "сложных" регионов FDR падает до 0.87%.

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

Влияние фильтрации на полученные значения

В нашей структуре пайлайна по умолчанию используются два типа фильтрации, предлагаемые утилитой GATK: это т.н. хардфильтрация (hard filtering) и фильтрация с помощью нейросетей (CNN подпрограммы).

На рис. 5 можно увидеть чувствительность и FDR для прогонов пайплайна с различными методиками фильтрации. Видно, что FDR удается снизить с 5,46% до 3,30% при использовании хардфильтрации. При этом чувствительность падает с 93,23% до 86,17%. Также видно, что совместное использование хардфильтрации и CNN фильтрации приводит к наименьшему значению FDR, хотя добавочный вклад CNN методов к отфильтровыванию ложноположительных вариантов в дополнение к хардфильтрации невелик (FDR падает с 3,30% до 3,14%).

Количество ложноположительных вариантов, отнесенное к исследуемой геномной дистанции в млн пар оснований, в целом повторяет картину FDR:

Рисунок 5. Столбчатая диаграмма чувствительности (слева), FDR (посередине) и количества ложноположительных вариантов на 1 млн п.о. (TN per Mbp, справа) герминативного пайплайна в зависимости от способа фильтрации вариантов. На каждой из диаграмм рассматриваются (1) варианты без фильтрации, (2) варианты, прошедшие фильтрацию с помощью CNNScoreVariants+FilterVariantTranches (CNN фильтрация), (3) хардфильтрацию и (4) обе фильтрации. Во всех случаях рассматриваются только регионы высокой надежности, заявленные консорциумом GIAB.

Тестирование соматического пайплайна

Для детекции соматических вариантов нами используется коллер Mutect2 из тулов GATK.

Тестовые образцы

В качестве источника стандартных образцов, содержащих известный набор соматических мутаций, мы использовали данные консорциума SEQC2 (²Fang et al., 2021). Его подход в некоторой степени аналогичен подходу, примененному для создания стандартов детекции герминальных вариантов, описанному выше.

Этот консорциум содержит как данные полногеномного секвенирования, так и полноэкзомные данные. При этом и те, и другие секвенировались в различных местах и институтах. Нами использовались данные архива SRR7890883 (WES_IL_T_1) в качестве образца опухоли и данные архива SRR7890874 (WES_IL_N_1) в качестве соответствующего нормального образца. Проводился как анализ опухолевого образца в т.н. tumor-only режиме работы пайплайна (без соответствующего неопухолевого или "нормального" образца), так и анализ в режиме tumor-normal (с применением неопухолевого образца).

Сравнение с "золотыми" вариантами

При сравнении вариантов внутри регионов надежности и покрытия панели чувствительность составила 97.44% для пайплайна в режиме tumor-only. Рассмотрение незадетектированных вариантов показывает, что основной причиной ложноотрицательных результатов является низкая частота мутаций. Пример такого незадетектированного варианта приведен на рис. 6:

Рисунок 6. Скриншот программы IGV, показывающий пример варианта, незадетектированного по причине низкой частоты соматического варианта. Самые верхние трэки представлены вариантами "золотого стандарта" и задетектированными вариантами. Ниже расположены выравнивания. В самом низу есть трэк регионов, которые заявлены как покрываемые используемой панелью Sureselect, и трэк регионов высокой достоверности.

Эта черта характерна именно для соматического пайплайна, но не для герминального. В случае герминальных вариантов мы ожидаем хотя бы половину молекул ДНК содержать вариант (для гетерозиготы), однако доля молекул, несущих соматический вариант, может быть очень мала.

Применение анализа регионов "сложности" от консорциума "Genome in a Bottle" показывает, что сужение области анализа до "несложных" регионов позволяет повысить чувствительность до 97.53%:

Рисунок 7. Столбчатая диаграмма чувствительности соматического пайплайна. Показана чувствительность (в %) задетектированных вариантов "золотого стандарта", (1) попадающих в регионы панели секвенирования ("Все варианты") и (2) попадающих в регионы панели и не попадающих в "сложные" регионы по классификациии консорциума GIAB. Во всех случаях рассматриваются только регионы высокой надежности, заявленные консорциумом SEQC2.

Специфичность, фильтрация и использование неопухолевого образца

Попытка посчитать метрики, характеризующие ложноположительные результаты (например, FDR), приводит к очень высоким относительным количествам ложных результатов в tumor-only анализе, что связано с тем, что герминальных вариантов намного больше, чем соматических. Фильтрация с помощью баз известных вариантов (dbSNP) и утилит вроде FilterMutectCalls позволяет частично отфильтровать эти варианты.

Конечно, лучше всего для этого пользоваться анализом в режиме tumor-normal, хотя в клинической практике далеко не всегда имеется такая возможность.

Результаты по чувствительности и FDR обоих режимов пайплайна с использованием фильтрации и без показаны на рис. 8:

Рисунок 8. Столбчатая диаграмма чувствительности (слева) и FDR (справа) соматического пайплайна в зависимости от фильтрации и наличия неопухолевого образца. На каждой из диаграмм рассматриваются варианты без фильтрации и прошедшие фильтрацию, полученные в tumor-only или tumor-normal анализе. Во всех случаях рассматриваются только регионы высокой надежности, заявленные консорциумом SEQC2.

По графикам видно, что применение фильтрации (используя все доступные фильтры FilterMutectCalls) понижает чувствительность как в tumor-only, так и в tumor-normal режиме анализа. Например, для tumor-only анализа чувствительность падает с 97,44% до 89,42%.

FDR имеет наибольшее значение для tumor-only анализа без фильтрации (98,48%). Это означает, что из всех обнаруживаемых с помощью Mutect2 вариантов лишь малая часть (менее 2%) является истинно соматическими вариантами. FDR падает до уровня 5,85% при использовании сиквенсов неопухолевой ткани. Применение же методов фильтрации in silico понижает FDR лишь незначительно до уровней выше 80%.

Команды, использованные при генерации файлов

Сравнение вариантов производилось утилитой bcftools isec. Предварительно файлы нормализовались с помощью bcftools norm, чтобы унифицировать представление вариантов. Кроме того, VCF "сужались" до необходимых регионов с помощью bedtools intersect.

Например, для сравнения герминальных "золотых вариантов", прошедших хардфильтрацию и CNN фильтрацию, нужны следующие команды (использовались при построении рис. 5):

cat pipeline_output.vcf | bcftools norm -a -f reference.fasta -m -both -O v | bcftools norm -d exact -f reference.fasta -O v | bcftools sort  > norm_pipeline_out.vcfbedtools intersect -u -header -a norm_pipeline_out.vcf -b nexterarapid_NA12878_hg38_lifted.regchrom.sorted.bed   | bedtools intersect -u -header -a - -b HG001_GRCh38_1_22_v4.2.1_benchmark.bed | bgzip -c > norm_pipeline_out.considered_regions.vcf.gzbcftools view  --exclude 'FILTER~"HDFLT_FS"||FILTER~"HDFLT_MQ"||FILTER~"HDFLT_MQRankSum"||FILTER~"HDFLT_QD"||FILTER~"HDFLT_QUAL"||FILTER~"HDFLT_ReadPosRankSum"||FILTER~"HDFLT_SOR"||FILTER~"CNN_2D_INDEL_Tranche_99.40_100.00"||FILTER~"CNN_2D_SNP_Tranche_99.95_100.00"' norm_pipeline_out.considered_regions.vcf.gz  | bgzip -c > norm_pipeline_out_HARDnCNN.considered_regions.vcf.gztabix norm_pipeline_out_HARDnCNN.considered_regions.vcf.gzgunzip -c HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz | bcftools norm -a -f reference.fasta -m -both -O v | bcftools norm -d exact -f reference.fasta -O v | bcftools sort  | bgzip -c > HG001_GRCh38_1_22_v4.2.1_benchmark.norm.vcf.gzbedtools intersect -u -header -a HG001_GRCh38_1_22_v4.2.1_benchmark.norm.vcf.gz -b nexterarapid_NA12878_hg38_lifted.regchrom.sorted.bed   | bedtools intersect -u -header -a - -b HG001_GRCh38_1_22_v4.2.1_benchmark.bed | bgzip -c > giab.considered_regions.vcf.gzbcftools isec norm_pipeline_out_HARDnCNN.considered_regions.vcf.gz giab.considered_regions.vcf.gz -p norm_pipeline_out_HARDnCNN_vs_giab

На выходе получается папка norm_pipeline_out_HARDnCNN_vs_giab, содержащая общие варианты и варианты, уникальные для полученного и скачанного файлов VCF.

Здесь использовались следующие файлы, полученные с FTP проекта GIAB:

• HG001_GRCh38_1_22_v4.2.1_benchmark.vcf.gz - файл VCF с вариантами;
• HG001_GRCh38_1_22_v4.2.1_benchmark.bed - регионы высокой достоверности;
• nexterarapid_NA12878_hg38_lifted.regchrom.sorted.bed - файл панели, примененной при получении использованных здесь файлов.

При исследовании соматических вариантов используется аналогичный подход. Например, чтобы сравнить фильтрованные варианты с вариантами "золотого" стандарта, нужны следующие команды (использовались при построении рис. 8):

cat pipeline_output.vcf | bcftools norm -a -f reference.fasta -m -both -O v | bcftools norm -d exact -f reference.fasta -O v | bcftools sort  > norm_pipeline_out_nofilt.vcfbedtools intersect -u -header -a norm_pipeline_out_nofilt.vcf -b S07604624_Covered.bed   | bedtools intersect -u -header -a - -b High-Confidence_Regions_v1.2.bed | bgzip -c > norm_pipeline_out_nofilt.considered_regions.vcf.gzbcftools view -f .,PASS norm_pipeline_out_nofilt.considered_regions.vcf.gz | bgzip -c > norm_pipeline_out_filt.considered_regions.vcf.gztabix norm_pipeline_out_filt.considered_regions.vcf.gzbcftools merge high-confidence_sINDEL_in_HC_regions_v1.2.vcf.gz high-confidence_sSNV_in_HC_regions_v1.2.vcf.gz | bcftools norm -f reference.fasta -m -both -O v | bcftools norm -d exact -f reference.fasta -O v | bcftools sort > somhiconf.vcfbedtools intersect -u -header -a somhiconf.vcf -b S07604624_Covered.bed   | bedtools intersect -u -header -a - -b High-Confidence_Regions_v1.2.bed | bgzip -c > somhiconf.considered_regions.vcf.gztabix somhiconf.considered_regions.vcf.gzbcftools isec norm_pipeline_out_filt.considered_regions.vcf.gz somhiconf.considered_regions.vcf.gz -p norm_pipeline_out_filt_vs_somhiconf

Аналогично, здесь использовались следующие файлы, полученные с FTP проекта GIAB:

• high-confidence_sINDEL_in_HC_regions_v1.2.vcf.gz и high-confidence_sSNV_in_HC_regions_v1.2.vcf.gz - варианты высокой достоверности;

• High-Confidence_Regions_v1.2.bed - регионы высокой достоверности;

• S07604624_Covered.bed - регионы панели Sureselect V7.

---

1. Чувствительность и специфичность↩
2. Fang L.T., Zhu B., Zhao Y. et al. Establishing community reference samples, data and call sets for benchmarking cancer mutation detection using whole-genome sequencing. Nat Biotechnol 39, 1151–1160 (2021)↩