среда, 7 октября 2015 г.

Идентификация сайтов SNP по данным NGS: Использование R функции exactSNP.

Функция exactSNP является составной частью библиотеки Rsubread, разработанной сотрудниками Отдела биоинформатики Университета Мельбурна (Мельбурн, Австралия) Wei Shi и Yang Liao. Библиотека находится в свободном доступе на официальном сайте проекта Bioconductor (http://bioconductor.org/packages/release/bioc/html/Rsubread.html). Там же можно посмотреть руководство пользователя (http://bioconductor.org/packages/release/bioc/manuals/Rsubread/man/Rsubread.pdf) и виньетку “Rsubread package: high-performance read alignment, quantification and mutation discovery” (http://bioconductor.org/packages/release/bioc/vignettes/Rsubread/inst/doc/Rsubread.pdf). Кроме того, для более углубленного изучения широких возможностей библиотеки Rsubread по анализу данных геномного секвенирования читателям можно порекомендовать дополнительное руководство “Subread/Rsubread users guide” (http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf), а также авторскую публикацию “The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote” в журнале Nucleic Acids Research за 2013 год (PubMed ID 23558742).
Модель.
Для идентификации сайтов SNP в функции exactSNP реализован простой подход, основанный на точном тесте Фишера. На первом этапе оценивается локальный фоновый шум – количество альтернативных и референсных нуклеотидов по каждой геномной позиции, покрытой ридами. После этого с помощью точного теста Фишера идентифицируются те геномные позиции, которые содержат альтернативные нуклеотиды в ридах и при этом статистически значимо выделяются на фоне шума. Поскольку при анализе целого генома количество позиций, содержащих альтернативные нуклеотиды, может быть очень большим, и поскольку некоторые геномные позиции могут быть покрыты большим количеством ридов, то для ускорения вычислений основной код функции exactSNP реализован на С.
Аргументы функции.
Функция exactSNP имеет тринадцать аргументов.
q readFile – файл, содержащий картированные риды. По умолчанию принимаются файлы формата SAM, однако может быть задан и BAM-файл, но в этом случае аргументу isBAM (см. ниже) необходимо присвоить значение TRUE.
q isBAM – параметр, позволяющий использовать в качестве исходных BAM-файлы. По умолчанию имеет значение FALSE.
q refGenomeFile – файл, содержащий референсную последовательность (целый геном или его фрагмент) в формате FASTA. Относительно предоставленной в этом файле последовательности будет анализироваться экспериментальный геном.
q SNPAnnotationFile – файл в формате VCF, содержащий аннотации известных сайтов SNP. Аргумент не является обязательным и по умолчанию имеет значение NULL. При необходимости такого рода файл может быть создан на основе данных, содержащихся в релевантных базах (например, NCBI dbSNP или 1000 Genomes). Поскольку данные об известных сайтах SNP нужны не только на этапе идентификации, но и на этапе аннотирования структурных вариантов экспериментального генома, то проблемам создания таких фалов и их использования будет посвящен самостоятельный раздел “Файлы формата VCF: организация, создание и использование”.
q outputFile – имя файла в формате VCF, в котором будут сохранены данные обо всех идентифицированных сайтах SNP, а также инсерциях и делециях.
q qvalueCutoff – пороговый уровень значения q, выше которого геномная позиция классифицируется как структурный вариант. Значение q рассчитывается как отрицательный десятичный логарифм уровня значимости –log10(p), получаемого из точного теста Фишера. По умолчанию этот аргумент имеет   значение   12, при   этом   предполагается, что   геномная   позиция   имеет   50-кратное покрытие ридами. Однако для каждой текущей геномной позиции функция автоматически меняет пороговый уровень q в зависимости от реального покрытия этой позиции ридами.
q minAllelicFraction – минимальная доля ридов, содержащих альтернативный нуклеотид по текущей геномной позиции, при превышении которой позиция будет классифицирована как структурный вариант. Аргумент может принимать значения от 0 (по умолчанию) до 1.
q minAllelicBases – минимальное количество ридов, содержащих альтернативный нуклеотид по текущей геномной позиции, при достижении или превышении которого позиция будет классифицирована как структурный вариант. Аргумент может принимать целочисленные значения (1 по умолчанию).
q minReads – минимально необходимое количество ридов, покрывающих геномную позицию со структурным вариантом, что бы эта позиция была включена в окончательный список вариантов. По умолчанию этот аргумент имеет значение 1.
q maxReads – максимально допустимое покрытие геномной позиции, содержащей структурный вариант. При превышении этого порога позиция будет исключена из дальнейшего анализа. Аргумент позволяет контролировать частоту ложно-положительных результатов, обусловленных ПЦР-артефактами. По умолчанию имеет значение 3000.
q minBaseQuality – минимальное качество прочтения нуклеотида, ниже которого нуклеотид не будет использован при идентификации сайта SNP. По умолчанию этот аргумент имеет значение 13.
q nTrimmedBases – количество нуклеотидов, удаляемых по концам ридов перед тем, как риды будут использованы в анализе. По умолчанию удаляется по 3 нуклеотида с каждого конца рида. Очевидно, что этому аргументу можно присвоить значение 0 если риды прошли соответствующий предпроцессинг.
q nthreads – количество потоков, задействованных при проведении параллельных вычислений. По умолчанию этот аргумент имеет значение 1.
Получение первичного списка сайтов SNP.
Ниже представлен обобщенный R код для идентификации сайтов SNP с помощью функции exactSNP.

rm(list = ls())
setwd("/media/gvv/NGS_Illumina_MiSeq/")
library(Rsubread)
library(parallel)
cl = makeCluster(detectCores(logical = FALSE) – 1)
bamFile = "Files_BAM/ExampleBAM.bam"
refGenome = "Reference_FASTA_files/hg19.fa"
SNPAnnoFile = "GSVs_References/GRCh37p13_SNPs144_LeukemiaGenes.vcf.gz"
vcfFile = "GSVs_Detected/Example_primarySNPs.vcf"
SNP = exactSNP(readFile = bamFile,
   isBAM = TRUE,
   refGenomeFile = refGenome,
   SNPAnnotationFile = SNPAnnoFile,
   outputFile = vcfFile,
   qvalueCutoff = 12,
   minAllelicFraction = 0,
   minAllelicBases = 1,
   minReads = 1,
   maxReads = 3000,
   minBaseQuality = 13,
   nTrimmedBases = 0,
   nthreads = length(cl))

Мы подали на функцию небольшой BAM-файл, содержащий только 54855 ридов. Кроме того, функции был предоставлен заранее подготовленный файл с аннотациями известных SNP, локализованных в интересующих нас генах (262355 сайтов SNP, принадлежащих 59 генам). Что бы ускорить вычисления мы воспользовались функцией makeCluster из библиотеки parallel и провели вычисления в три потока: функция detectCores из этой же библиотеки позволила установить, что процессор на нашем рабочем компьютере имеет четыре физических ядра и три из них (detectCores(logical = FALSE) – 1) было передано в распоряжение функции makeCluster.
После запуска кода на экране будет отображаться результат его выполнения, как показано на нижеследующем рисунке. Для упрощения на рисунке приведена только часть вывода. Обращаем внимание читателя на то, что в этом выводе в разделе настроек (exactSNP setting) появилась строка Flanking windows size, где указан размер геномных интервалов, фланкирующих анализируемый сайт нуклеотидной вариации. По этим интервалам проводится оценка локального фонового шума, которая необходима для точного теста Фишера.

Окно вывода результатов работы функции exactSNP
Конечным результатом работы функции является файл формата VCF, а обобщенная информация представлена в разделе Summary окна вывода. Название и место записи итогового файла VCF мы определили с помощью аргумента outputFile. Этот файл содержит краткую аннотацию идентифицированных сайтов SNP в формате VCFv4.0 и может быть использован для пост-фильтрации. Фрагмент данных, содержащихся в полученном файле, показан ниже:

##fileformat=VCFv4.0
##INFO=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##INFO=<ID=BGMM,Number=1,Type=Integer,Description="Number of mismatched bases in the background">
##INFO=<ID=BGTOTAL,Number=1,Type=Integer,Description="Total number of bases in the background">
##INFO=<ID=MM,Number=1,Type=String,Description="Number of supporting reads for each alternative allele">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=SR,Number=1,Type=String,Description="Number of supporting reads for variants">
#CHROM
POS
ID
REF
ALT
QUAL
FILTER
INFO
chr12
49427283
.
CTGCA
CA
1
.
INDEL;DP=516;SR=34
chr12
49427678
.
GCTGT
GT
1
.
INDEL;DP=555;SR=20
chr12
49433413
.
AAG
AG
1
.
INDEL;DP=361;SR=17
chr12
49436812
.
C
A,G,T
40
.
DP=572;MM=3,2,275;BGTOTAL=5706;BGMM=29
chr12
49442179
.
G
A
40
.
DP=172;MM=84;BGTOTAL=1656;BGMM=9

Преимущества и недостатки функции exactSNP.
В завершение мы попробуем, исходя из собственного опыта работы, выделить преимущества и недостатки функции exactSNP по сравнению с функцией callVariants. К плюсам функции exactSNP относятся:
1) функция управляется с помощью набора простых аргументов, при этом значения большинства из этих аргументов могут быть оставлены по умолчанию;
2) исходным файлом, содержащим картированные риды, может быть неиндексированный BAM-файл;
3) в качестве референсной последовательности используется обычный файл в формате FASTA, а не индексированная последовательность в виде объекта класса GmapGenome;
4) функция позволяет идентифицировать не только сайты однонуклеотидного полиморфизма, но и такие варианты простого нуклеотидного полиморфизма как небольшие инсерции и делеции.
В тоже время эта функция:
1) не позволяет проводить пост-фильтрацию полученных результатов без предварительной записи итогового файла VCF – пост-фильтрация может быть проведена только после повторной загрузки этого файла в рабочее пространство R среды;
2) менее чувствительна, чем функция callVariants, и некоторые реальные структурные варианты с ее помощью не детектируются;
3) требует в 5.8 раза больше времени для анализа одного и того же набора ридов относительно одной и той же референсной последовательности, чем функция callVariants.

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

Отправить комментарий