воскресенье, 4 октября 2015 г.

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

В этом сообщении мы обсудим возможности R функции callVariants при идентификации сайтов простого нуклеотидного полиморфизма (SNP) по данным, полученным с помощью NGS. При этом предполагается, что у исследователя уже имеются готовые BAM- и BAI-файлы с индексами, и что у него есть навыки работы в среде такой операционной системы как Linux.
Функция callVariants входит в сосав библиотеки VariantTools, разработанной сотрудниками компании Genentech (Сан Франциско, США) Michael Lawrence, Jeremiah Degenhardt и Robert Gentleman. Последний из соавторов этой библиотеки, Robert Gentleman, к тому же является одним из отцов-основателей языка программирования R. Библиотека находится в свободном доступе на официальном сайте проекта Bioconductor (http://bioconductor.org/packages/release/bioc/html/VariantTools.html). Там же можно посмотреть руководство пользователя (http://bioconductor.org/packages/release/bioc/manuals/VariantTools/man/ VariantTools.pdf) и виньетку “An introduction to VariantTools” (http://bioconductor.org/packages/release/bioc/vignettes/ VariantTools/inst/doc/VariantTools.pdf). Любознательным читателям можно порекомендовать презентацию Michael Lawrence (https://www.bioconductor.org/help/course-materials/2014/CSAMA2014/3_Wednesday/lectures/VariantCallingLecture.pdf), посвященную возможностям библиотеки VariantTools, а также небольшой учебный видеокурс по проблемам идентификации структурных вариаций генома, разработанный сотрудниками Bioconductor (https://www.youtube.com/watch?v=An-4mX77hl8).
Модель.
В основу функции callVariants положено представление о том, что появление ридов с альтернативными нуклеотидами может быть описано с помощью биномиального распределения:
где n это число ридов, покрывающих данную геномную позицию, k – число ридов, имеющих альтернативный нуклеотид по данной геномной позиции, при этом k может принимать целочисленные значения от 0 до n, p – вероятность появления рида с альтернативным нуклеотидом и q – вероятность появления рида с референсным нуклеотидом.
Альтернативный нуклеотид в риде может появиться как закономерно (данный нуклеотид действительно является структурным вариантом экспериментального генома), так и в силу случайных событий (например, ошибки секвенирования). Что бы установить реальный источник альтернативного нуклеотида можно воспользоваться тестом на отношение правдоподобий двух биномиальных моделей – модели, предполагающей, что альтернативный нуклеотид является закономерностью, и модели, говорящей, что мы имеем дело со случайностью:
где p1 это нижняя граница частоты встречаемости альтернативного нуклеотида, выше которой он может быть классифицирован как структурный вариант и p2 – вероятность появления альтернативного нуклеотида в силу случайных причин.
В том случае, когда соотношение LRT получается больше 1, предпочтение отдается первой модели – в рассматриваемой позиции действительно есть структурная вариация и эта позиция может быть классифицирована как сайт SNP. Если же соотношение получается меньше 1, то наблюдаемые несоответствия между референсной последовательностью и ридами носят случайный характер.
Аргументы функции.
Функция callVariants имеет четыре аргумента, причем три из них являются комплексными.
q x – индексированный BAM-файл, объект класса BamFile или объект класса VRanges, возвращаемый функцией tallyVariants.
q tally.param – параметры,   контролирующие   этап   сбора   количественной   информации   о   вариациях   нуклеотидной последовательности в ридах.
q calling.filters – фильтры, используемые на этапе идентификации потенциальных кандидатур на структурные варианты.
q post.filters – фильтры, используемые на этапе отбора реальных структурных вариантов из списка потенциальных кандидатур.
Исходные данные в виде BAM-файла могут быть поданы на функцию callVariants напрямую или же через создание объекта класса BamFile. Объект класса BamFile является не самим BAM-файлом, загруженным в рабочее пространство R среды, а лишь ссылкой на этот файл и файл с индексами. Использование такого объекта вместо прямой загрузки самого BAM-файла позволяет обрабатывать данные из BAM-файла с помощью разных методов без повторных загрузок в R среду файла с индексами. Ниже представлен код для создания такого объекта:

rm(list = ls())
library(Rsamtools)
dataBAM = BamFile(file = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/ExampleBAM.bam",
        index = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/ExampleBAM.bam.bai")
print(dataBAM)
class: BamFile
path: D:/Transcriptome/ExampleBAM.bam
index: D:/Transcriptome/ExampleBAM.bam.bai
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

Функция BamFile, использованная в нашем примере, входит в состав библиотеки Rsamtools, именно поэтому мы начали с загрузки данной библиотеки в рабочее пространство R среды. Однако если в рабочее пространство уже была загружена библиотека VariantTools, то загружать отдельно библиотеку Rsamtools нет необходимости, так как она импортируется вместе с библиотекой VariantTools.
Аргумент tally.param является комплексным и задается с помощью функции низкого уровня TallyVariantsParam. Эта функция имеет несколько аргументов:
q genome – геном, сжатый в объект класса GmapGenome.
q read_pos_breaks – параметр,  определяющий  размер  геномных  интервалов,  по  которым  будет проводиться  сбор количественной информации о вариациях нуклеотидной последовательности в ридах. По умолчанию имеет значение NULL. Если задается пользователем, то может принимать только целочисленные значения.
q high_base_quality – порог качества прочтения нуклеотида, выше которого нуклеотид считается надежно прочтенным. По умолчанию все нуклеотиды принимаются как высококачественные, однако конечный результат анализа можно улучшить, выставив реальный пороговый уровень качества. В случае с ридами, полученными на таких секвенаторах как Illumina v.1.8, авторы библиотеки VariantTools рекомендуют использовать порог качества в 23 балла.
q minimum_mapq – порог качества рида, выше которого рид будет включен, а ниже которого он будет исключен из анализа. По умолчанию этот аргумент имеет значение 13, однако может быть изменен по усмотрению пользователя.
q variant_strand – количество цепей, по которым должен быть обнаружен альтернативный нуклеотид, что бы позиция была классифицирована как структурный вариант. Может принимать значения 0, 1 или 2. Практический опыт работы показывает, что приемлемый по размеру список первичных кандидатур можно получить, присвоив этому аргументу значение 1.
q ignore_query_Ns – могут ли игнорироваться не идентифицированные в ридах нуклеотиды. Аргумент может принимать два логических значения: TRUE (по умолчанию) или FALSE. Если задается значение FALSE, то при анализе низкокачественных данных это может привести к получению слишком обширного списка идентифицированных структурных вариантов.
q ignore_duplicates – могут ли игнорироваться риды, помеченные как дубликаты. Аргумент может принимать два логических значения: TRUE (по умолчанию) или FALSE.
q mask – список “замаскированных” геномных интервалов. Список задается с помощью стандартного объекта класса GRanges. Все структурные варианты, попадающие в “замаскированные” геномные интервалы, исключаются из конечного перечня.
q keep_extra_stats – дополнительная статистика, получаемая на этапе сбора количественной информации о структурных вариантах. Аргумент может принимать два логических значения: TRUE (по умолчанию) или FALSE. Отключение этого аргумента позволяет уменьшить объем используемой оперативной памяти и сократить время, затрачиваемое на расчеты. Однако в некоторых случаях (например, при разведочном анализе) получаемая с помощью данного аргумента информация может быть полезной.
q read_length – ожидаемая длина ридов. Эта информация используется для расчета медианы расстояния от текущего рида до конца ближайшего рида. В том случае, если аргумент не задан, то его значение будет основано на случайной выборке ридов из загруженного BAM-файла. Если длина ридов варьирует, то данный аргумент не используется и не проводится расчет всех зависимых от него статистик.
Ключевой проблемой при работе с функцией TallyVariantsParam является аргумент genome, который может быть задан только как объект класса GmapGenome. Объекты такого класса являются специализированными базами данных, где хранится целый геном (или иная нуклеотидная последовательность), индексированный для последующего использования в картировании, выравнивании и/или обобщении последовательностей мРНК и EST, а также коротких ридов с помощью программ GMAP и/или GSNAP (http://research-pub.gene.com/gmap). Функция TallyVariantsParam использует программу GSNAP, предназначенную для работы только с короткими ридами, но не в ее исходном варианте, а в виде R версии, скомпилированной в библиотеку gmapR, доступную на официальном сайте проекта Bioconductor (https://www.bioconductor.org/packages/release/bioc/html/gmapR.html).
К сожалению, объектов класса GmapGenome, готовых к использованию в R среде, нет (за исключением некоторых упрощенных примеров, импортируемых вместе с библиотеками gmapR и VariantTools), хотя для работы с оригинальными версиями программ GMAP и GSNAP они и существуют (например, 5,5-Гб архив пре-компилированной референс-сборки hg19 генома человека доступен для скачивания по адресу http://research-pub.gene.com/gmap/genomes/hg19.tar.gz). Однако пользователь может самостоятельно создать нужный ему объект прямо из R среды, воспользовавшись для этого либо объектом класса BSgenome, либо файлом формата FASTA с необходимой нуклеотидной последовательностью.
Объекты класса BSgenome являются специализированными базами данных, хранящими нуклеотидные последовательности целых геномов разных видов живых организмов. Такого рода объекты широко используются при работе с функциями, зависимыми от библиотеки Biostrings, причем уже имеющиеся объекты класса BSgenome преобразованы в стандартные библиотеки, доступные через официальный сайт проекта Bioconductor. Что бы воспользоваться предоставляемыми возможностями мы первоначально установим библиотеку BSgenome, которая создаст нам необходимую инфраструктуру. После этого мы сформируем запрос, что бы выяснить, какие геномы преобразованы в объекты класса BSgenome.

source("https://bioconductor.org/biocLite.R")
biocLite("BSgenome")
library(BSgenome)
GENOMES = available.genomes(splitNameParts = TRUE, type = getOption("pkgType"))
str(GENOMES)
'data.frame': 71 obs. of 5 variables:
$ pkgname                   : chr "BSgenome.Alyrata.JGI.v1" "BSgenome.Amellifera.BeeBase.assembly4" ...
$ organism                   : Factor w/ 22 levels "Alyrata", "Amellifera",..: 1 2 2 2 3 3 4 4 4 4 ...
$ provider                     : Factor w/ 7 levels "BeeBase", "JGI",..: 2 1 7 7 5 5 7 7 7 7 ...
$ provider_version       : chr "v1" "assembly4" "apiMel2" "apiMel2" ...
$ masked                     : logi FALSE FALSE FALSE TRUE FALSE FALSE ...
print(GENOMES[1:5, ])
pkgname
organism
provider
provider_version
masked
1
BSgenome.Alyrata.JGI.v1
Alyrata
JGI
v1
FALSE
2
BSgenome.Amellifera.BeeBase.assembly4
Amellifera
BeeBase
assembly4
FALSE
3
BSgenome.Amellifera.UCSC.apiMel2
Amellifera
UCSC
apiMel2
FALSE
4
BSgenome.Amellifera.UCSC.apiMel2.masked
Amellifera
UCSC
apiMel2
TRUE
5
BSgenome.Athaliana.TAIR.04232008
Athaliana
TAIR
04232008
FALSE

Как мы видим, всего доступен 71 геном (см. поле pkgname созданного нами объекта GENOMES), в том числе доступны 4 сборки генома человека: от hg17 до hg38. Установим одну из библиотек с геномом человека и загрузим ее в рабочее пространство R среды, что позволит нам создать ссылку на соответствующий объект класса BSgenome с его кратким описанием (для простоты вывод описательной информации сокращен).

source("https://bioconductor.org/biocLite.R")
biocLite("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)
print(Hsapiens)
Human genome
|
| organism: Homo sapiens (Human)
| provider: UCSC
| provider version: hg19
| release date: Feb. 2009
| release name: Genome Reference Consortium GRCh37
| 93 sequences:
|  chr1  chr2     chr3     chr4     chr5     chr6     chr7     chr8
|  ...      ...         ...         ...         ...         ...         ...         ...
|  chrUn_gl000245   chrUn_gl000246   chrUn_gl000247   chrUn_gl000248   chrUn_gl000249
| (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)
str(Hsapiens)
Formal class 'BSgenome' [package "BSgenome"] with 17 slots
..@ pkgname
: chr "BSgenome.Hsapiens.UCSC.hg19"
..@ single_sequences
: Formal class 'TwobitNamedSequences' [package "BSgenome"] with 1 slot
.. .. ..@ twobitfile
: Formal class 'TwoBitFile' [package "rtracklayer"] with 1 slot
.. .. .. .. ..@ resource
: chr "C:/Users/admin/Documents/R/win-library/3.1/BSgenome.Hsapiens. UCSC.hg19/
extdata/single_sequences.2bit"
..@ multiple_sequences
: Formal class 'RdaCollection' [package "XVector"] with 2 slots
.. .. ..@ dirpath
: chr "C:/Users/admin/Documents/R/win-library/3.1/BSgenome.Hsapiens.UCSC.hg19/extdata"
.. .. ..@ objnames.
: chr(0)
..@ source_url
: chr http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/
..@ user_seqnames
: Named chr [1:93] "chr1" "chr2" "chr3" "chr4" ...
.. ..- attr(*, "names")= chr [1:93] "chr1" "chr2" "chr3" "chr4" ...
..@ injectSNPs_handler
: Formal class 'InjectSNPsHandler' [package "BSgenome"] with 4 slots
.. .. ..@ SNPlocs_pkgname
: chr(0)
.. .. ..@ getSNPcount
: function ()
.. .. ..@ getSNPlocs
: function ()
.. .. ..@ seqname_translation_table
: chr(0)
..@ .seqs_cache
: <environment: 0x0d6ac240>
..@ .link_counts
: <environment: 0x0d6ab434>
..@ nmask_per_seq
: int 0
..@ masks
: Formal class 'RdaCollection' [package "XVector"] with 2 slots
.. .. ..@ dirpath
: chr NA
.. .. ..@ objnames
: chr(0)
..@ organism
: chr "Homo sapiens"
..@ species
: chr "Human"
..@ provider
: chr "UCSC"
..@ provider_version
: chr "hg19"
..@ release_date
: chr "Feb. 2009"
..@ release_name
: chr "Genome Reference Consortium GRCh37"
..@ seqinfo
: Formal class 'Seqinfo' [package "GenomeInfoDb"] with 4 slots
.. .. ..@ seqnames
: chr [1:93] "chr1" "chr2" "chr3" "chr4" ...
.. .. ..@ seqlengths
: int [1:93] 249250621 243199373 198022430 191154276 180915260 171115067 ...
.. .. ..@ is_circular
: logi [1:93] FALSE FALSE FALSE FALSE FALSE FALSE ...
.. .. ..@ genome
: chr [1:93] "hg19" "hg19" "hg19" "hg19" ...

Далее мы сожмем имеющийся объект класса BSgenome в объект класса GmapGenome, для чего мы воспользуемся одноименной функцией GmapGenome из библиотеки gmapR. Если у нас уже установлена библиотека VariantTools, то нам не нужно отдельно устанавливать библиотеку gmapR, так как она будет импортирована как зависимость библиотеки VariantTools.

library(gmapR)
genomePath = "/media/gvv/NGS_Illumina_MiSeq/Reference_Genomes/"
humanBSgenome = GmapGenome(genome = Hsapiens,
       directory = GmapGenomeDirectory(path = genomePath, create = TRUE),
       name = "BSgenomeGRCh37",
       create = TRUE)

В качестве аргументов для функции GmapGenome нужно указать название генома (видовая принадлежность объекта класса BSgenome, приведенная в названии соответствующей библиотеки), директорию, куда будет записана новая база данных, и название этой базы. Процесс создания новой базы данных будет отображаться на экране, а по его окончанию можно будет получить краткую характеристику базы, вызвав объект humanBSgenome. В дальнейшем нам достаточно будет указать название базы и путь к ней, что бы воспользоваться этой базой повторно. Следует отметить, что преобразование целого генома в объект класса GmapGenome требует значительных вычислительных возможностей компьютерной техники и временных затрат.
Альтернативным источником нуклеотидной последовательности целого генома может быть файл в формате FASTA, полученный различными путями. Например, референс-сборка hg19 генома человека может быть скачана с ftp-сервера геномного браузера UCSC в виде единого файла в формате 2bit (ftp://hgdownload.soe.ucsc.edu/goldenPath/currentGenomes/Homo_sapiens/ bigZips/hg19.2bit). Формат файлов 2bit является одним из стандартов хранения и передачи больших нуклеотидных последовательностей в компактном виде. После скачивания полученный 2bit-файл легко преобразуется в стандартный файл формата FASTA с помощью утилиты twoBitToFa. Эта утилита доступна для скачивания с сайта геномного браузера UCSC по адресу http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa и работает в среде операционной системы Linux (там же можно получить версии утилиты, написанные под UNIX или MacOSX). Ниже приведен пример набора команд для запуска утилиты и конвертации форматов файлов:

pwd
cd /media/gvv/NGS_Illumina_MiSeq/Reference_FASTA_files
$ chmod +x twoBitToFa
$ ./twoBitToFa hg19.2bit hg19.fa

После того, как нужный нам файл формата FASTA получен, преобразуем его в объект класса GmapGenome с помощью соответствующей функции. Будем считать, что библиотека gmapR уже загружена в рабочее пространство R среды. Более того, при загрузке библиотеки gmapR в рабочее пространство загружается также библиотека rtracklayer, возможности которой нам понадобятся.

dataFASTA = "/media/gvv/NGS_Illumina_MiSeq/Reference_FASTA_files/hg19.fa"
fastaFile = rtracklayer::FastaFile(dataFASTA)
genomePath = "/media/gvv/NGS_Illumina_MiSeq/Reference_Genomes/"
humanFASTAgenome = GmapGenome(genome = fastaFile,
 directory = GmapGenomeDirectory(path = genomePath, create = TRUE),
 name = "FASTAGRCh37",
 create = TRUE)

Аналогичным способом можно преобразовать в объект класса GmapGenome не весь геном, а только его часть: несколько хромосом, отдельно взятую хромосому или же набор геномных интервалов. Сиквенсы индивидуальных хромосом в формате FASTA могут быть получены с ftp-сервера геномного браузера UCSC или из другого релевантного источника. При необходимости несколько файлов формата FASTA с индивидуальными хромосомами могут быть объединены в один и только после этого поданы на функцию GmapGenome. В качестве примера мы преобразуем в объект класса GmapGenome хромосому 12 человека, сиквенс которой в формате FASTA мы скачали с ftp-сервера геномного браузера UCSC по адресу ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr12.fa.gz

FASTAchr12 = "/media/gvv/NGS_Illumina_MiSeq/Reference_FASTA_files/chr12.fa"
fastaFile = rtracklayer::FastaFile(FASTAchr12)
genomePath = "/media/gvv/NGS_Illumina_MiSeq/Reference_Genomes/"
humanFASTAchr12 = GmapGenome(genome = fastaFile,
          directory = GmapGenomeDirectory(path = genomePath, create = TRUE),
          name = "hg19Chr12",
          create = TRUE)

В тех же случаях, когда не нужен весь геном и не нужны целые хромосомы, а лишь набор геномных интервалов (например, участков генома, соответствующих изучаемым генам), первоначально необходимо создать файл, содержащий нуклеотидные последовательности таких участков в формате FASTA. Из множества подходов, позволяющих это сделать, мы рассмотрим лишь один – мы извлечем нужные нам последовательности из объекта класса BSgenome и сохраним их в файле формата FASTA. Код представлен ниже:

rm(list = ls())
setwd("/media/gvv/NGS_Illumina_MiSeq/")
library(BSgenome.Hsapiens.UCSC.hg19)
LeukemiaGenes = read.table(file = "Annotation of genes of the myeloid leukemias.txt",
         header = TRUE,
         sep = "\t",
         as.is = TRUE,
         quote = "\"")
str(LeukemiaGenes)
'data.frame': 59 obs. of 11 variables:
 $ Approved_symbol
: chr  "ABL1" "ASXL1" ...
 $ Approved_name
: chr  "ABL proto-oncogene 1" "additional sex combs like transcriptional regulator 1" ...
 $ Chromosome
: chr  "chr9" "chr20" ...
 $ Band
: chr  "q34.12" "q11.21" ...
 $ Locus_type
: chr  "gene with protein product" "gene with protein product" ...
 $ RefSeq_accession
: chr  "NM_007313" "NM_015338" ...
 $ Ensembl_gene_ID
: chr  "ENSG00000097007" "ENSG00000171456" ...
 $ UCSC_gene_ID
: chr  "uc004bzw.3" "uc002wxs.3" ...
 $ Strand
: chr  "+" "+" ...
 $ Gene_start
: int  133589333 30946147 108093559 76760356 ...
Regions = GRanges(seqnames = Rle(LeukemiaGenes$Chromosome),
        ranges = IRanges(start = LeukemiaGenes$Gene_start,
end = LeukemiaGenes$Gene_end,
width = LeukemiaGenes$Gene_end - LeukemiaGenes$Gene_start + 1,
names = LeukemiaGenes$Approved_symbol))
print(Regions)
GRanges object with 59 ranges and 0 metadata columns:
seqnames
ranges
strand
<Rle>
<IRanges>
<Rle>
ABL1
chr9
[133589333,
133763062]
*
ASXL1
chr20
[  30946147,
31027122]
*
ATM
chr11
[108093559,
108239826]
*
ATRX
chrX
[  76760356,
77041755]
*
BCOR
chrX
[  39909068,
40036582]
*
...
...
...
...
TET2
chr4
[106067032,
106200973]
*
TP53
chr17
[    7565097,
7590868]
*
U2AF1
chr21
[  44513066,
44527697]
*
WT1
chr11
[  32409321,
32457176]
*
ZRSR2
chrX
[  15808595,
15841383]
*
  -------
  seqinfo: 19 sequences from an unspecified genome; no seqlengths
LeukemiaGenesSeq = getSeq(x = Hsapiens, Regions + 1000)
writeXStringSet(x = LeukemiaGenesSeq,
 file = "Reference_FASTA_files/LeukemiaGenes.fa",
 format = "fasta",
 width = 50)

В данном случае мы загрузили в рабочее пространство R среды библиотеку BSgenome.Hsapiens.UCSC.hg19, тем самым получив доступ к объекту класса BSgenome, где хранится сиквенс всего генома человека сборки hg19. Далее, через стандартную функцию read.table мы загрузили заранее подготовленную таблицу, содержащую информацию об интересующих нас генах, и преобразовали геномные координаты этих генов в объект класса GRanges. После этого с помощью функции getSeq мы извлекли нужные нам интервалы генома (плюс по 1000 фланкирующих нуклеотидов) из объекта класса BSgenome и сохранили их в новом объекте LeukemiaGenesSeq, который является объектом класса DNAStringSet – одном из ключевых классов, используемых функциями, зависимыми от библиотеки Biostrings. В завершение нам осталось сохранить полученный объект LeukemiaGenesSeq в виде файла формата FASTA, для чего мы воспользовались функцией writeXStringSet из библиотеки Biostrings, а затем преобразовать его в объект класса GmapGenome как это мы уже неоднократно делали ранее. Очевидно, что если мы захотим воспользоваться объектом класса GmapGenome, созданным таким путем, то и имеющийся у нас BAM-файл должен быть получен подходящим образом.
Теперь мы можем полностью задать параметры, контролирующие этап сбора количественной информации о вариациях нуклеотидной последовательности ридов в нашем BAM-файле. При этом в качестве примера мы присвоим аргументу genome имя ранее созданного нами объекта hg19Chr12 класса GmapGenome:

genome = "/media/gvv/NGS_Illumina_MiSeq/Reference_Genomes/hg19Сhr12"
Parameters = TallyVariantsParam(genome = genome,
     read_pos_breaks = NULL,
     high_base_quality = 23L,
     minimum_mapq = 13L,
     variant_strand = 1L,
     ignore_query_Ns = TRUE,
     ignore_duplicates = TRUE,
     keep_extra_stats = FALSE)

На наш взгляд, полученный объект Parameters лучше всего сразу же подать не на аргумент tally.param функции callVariants, а использовать его как один из аргументов функции tallyVariants, с помощью которой можно сначала собрать количественную информации о ридах BAM-файла и лишь после этого задействовать функцию callVariants. В конечном итоге такая стратегия позволяет сократить время, затрачиваемое на анализ (благодаря параллелизации вычислений, которые предусмотрены для этой функции в виде отдельного аргумента).
Функция tallyVariants имеет только три аргумента:
q x – индексированный BAM-файл или объект класса BamFile.
q param – параметры, контролирующие этап сбора количественной информации о вариациях нуклеотидной последовательности в ридах. Этому аргументу можно присвоить имя объекта Parameters, созданного нами ранее, или же задать параметры заново с помощью функции TallyVariantsParam.
q BPPARAM – объект виртуального класса BiocParallelParam, создаваемого функциями библиотеки BiocParallel. Объекты данного класса являются хранилищами параметров, контролирующих параллельные вычисления.
Пример кода, где задействована функция tallyVariants, представлен ниже:

library(BiocParallel)
talliedVariants = tallyVariants(x = dataBAM,
          param = Parameters,
          BPPARAM = MulticoreParam(workers = 4))

Мы загрузили в рабочее пространство R среды дополнительную библиотеку BiocParallel, так как функция MulticoreParam, необходимая для параллельных вычислений, является составной частью этой библиотеки. Полученный нами объект talliedVariants принадлежит к классу VRanges, который похож на класс GRanges, но специализирован на хранении данных по нуклеотидным вариациям. Этот объект может напрямую подаваться на аргумент x функции callVariants вместо индексированного BAM-файл или объекта класса BamFile (см. выше).

print(talliedVariants[1:5, ])
VRanges object with 21516 ranges and 0 metadata columns:
seqnames
ranges
strand
ref
alt
totalDepth
refDepth
altDepth
<Rle>
<IRanges>
<Rle>
<character>
<characterOrRle>
<integerOrRle>
<integerOrRle>
<integerOrRle>
[1]
chr12
[  712147,
712147]
+
T
C
1
1
0
[2]
chr12
[3625022,
3625022]
+
G
A
1
0
1
[3]
chr12
[3625071,
3625071]
+
T
G
0
0
0
[4]
chr12
[4051269,
4051269]
+
G
T
1
0
1
[5]
chr12
[4051757,
4051757]
+
C
A
1
0
1
  -------
  seqinfo: 93 sequences from 2 genomes (hg19Chr12, NA)
  hardFilters: NULL

Аргумент calling.filters так же является комплексным и задается с помощью функции низкого уровня VariantCallingFilters. Эта функция имеет три аргумента:
q read.count – пороговое количество высококачественных ридов с альтернативным нуклеотидом, при превышении которого позиция может быть классифицирована как структурный вариант. По умолчанию этот аргумент имеет значение 2, что позволяет отличить реальные структурные варианты от ошибок секвенирования в тех случаях, когда геномная позиция имеет низкое покрытие ридами.
q p.lower – нижняя граница частоты встречаемости альтернативного нуклеотида, выше которой он может быть классифицирован как структурный вариант. По умолчанию этому аргументу присвоено значение 0.2, что позволяет идентифицировать гетерозиготные суб-клоны, встречающиеся в популяции клеток с частотой 50%. Пользователь может изменить эту нижнюю границу по своему усмотрению.
q p.error – вероятность ошибки секвенирования. По умолчанию этому аргументу присвоено значение 1/1000, что вполне приемлемо для данных, полученных с помощью секвенаторов типа Illumina v.1.8 и пороге высококачественного прочтения в 23 балла.
Аргумент post.filters необходим для настраивания фильтров (правил отбора), позволяющих из списка первичных кандидатур отобрать те из них, которые являются реальными структурными вариантами. Этот аргумент задается с помощью функции VariantPostFilters, которая в настоящее время имеет только один аргумент max.nbor.count – максимально допустимое число соседей, представленных альтернативными нуклеотидами. Данный аргумент позволяет оценить насыщенность локального окружения геномной позиции, претендующей на структурный вариант, альтернативными нуклеотидами, и убрать те случаи, для которых этот показатель превышает заданный порог (0.1 по умолчанию).
Следует отметит, что задание аргумента post.filters не является обязательным в контексте использования функции callVariants, так как фильтрация кандидатур может быть произведена после их идентификации и с привлечением других подходов, нежели те, что предоставляются библиотекой VariantTools. Тем не менее, пост-фильтрация является настолько важным этапом в идентификации SNP по данным геномного секвенирования, что мы посвятим этому вопросу отдельный раздел.
Получение первичного списка сайтов SNP.
Теперь мы обобщим все наши рассуждения, изложенные выше, в виде единого R кода. При этом мы опустим этап создания объекта класса GmapGenome, предполагая, что он уже создан по одной из выше описанных стратегий, и используем наш объект hg19Chr12.

rm(list = ls())
setwd("/media/gvv/NGS_Illumina_MiSeq/")
library(VariantTools)
library(BiocParallel)
dataBAM = BamFile(file = "Files_BAM/ExampleBAM.bam",
        index = "Files_BAM/ExampleBAM.bam.bai")
genome = "Reference_Genomes/hg19Chr12"
Parameters = TallyVariantsParam(genome = genome,
     read_pos_breaks = NULL,
     high_base_quality = 23L,
     minimum_mapq = 13L,
     variant_strand = 1L,
     ignore_query_Ns = TRUE,
     ignore_duplicates = TRUE,
     keep_extra_stats = FALSE)
talliedVariants = tallyVariants(x = dataBAM,
          param = Parameters,
          BPPARAM = MulticoreParam(workers = 4))
preFilters = VariantCallingFilters(read.count = 2L, p.lower = 0.2, p.error = 1/1000)
postFilters = VariantPostFilters(max.nbor.count = 0.1)
primaryVariants = callVariants(x = talliedVariants,
           tally.param = Parameters,
           calling.filters = preFilters,
           post.filters = postFilters)

Полученный нами объект primaryVariants содержит первичных список структурных вариаций нуклеотидной последовательности в нашем экспериментальном геноме. Этот объект принадлежит к тому же классу, что и объект talliedVariants (фактически он им и является, но после первичной фильтрации с помощью аргументов calling.filters и post.filters). Он может быть подвергнут дальнейшей фильтрации, преобразованию и/или сохранению. Как уже говорилось, проблеме пост-фильтрации будет посвящен отдельный раздел. А в завершение этого раздела мы покажем лишь один из путей сохранения объектов класса VRanges – через преобразование и сохранение их в файле формата VCF. Код представлен ниже.

sampleNames(primaryVariants) = "Sample 1"
fileVCF = as(primaryVariants, "VCF")
print(fileVCF)
class: ExpandedVCF
dim: 57 1
rowData(vcf):
  GRanges with 4 metadata columns: REF, ALT, QUAL, FILTER
info(vcf):
  DataFrame with 0 columns:
geno(vcf):
  SimpleList of length 3: AD, DP, FT
geno(header(vcf)):
        Number     Type       Description
AD   2               Integer    Allelic depths (number of reads in each observed allele)
DP   1               Integer    Total read depth
FT    1               String      Variant filters
writeVcf(fileVCF, "VCFexample.vcf")

Вот таким простым способом мы сохранили полученные нами данные в формате VCF (версия v.4.1), который на сегодняшний момент является одним из самых оптимальных форматов для дальнейшей работы с нуклеотидными вариациями.


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

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