Функция 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.
Комментариев нет:
Отправить комментарий