вторник, 20 октября 2015 г.

Инсталляция библиотек R функций VariantTools и VariantAnnotation в среде операционной системы Linux Ubuntu.

Библиотеки VariantTools и VariantAnnotation являются собраниями R функций, с помощью которых можно проводить разнообразные манипуляции со структурными вариациями генома – от идентификации сайтов структурных вариаций по данным NGS до аннотирования таких сайтов. В виде исходного кода библиотеки могут быть получены с официального сайта проекта Bioconductor, а библиотека VariantAnnotation доступна еще и в виде бинарного пакета под Windows.
При работе в операционной системе Linux Ubuntu (или аналогичной) большинство библиотек R функций устанавливаются одним из стандартных путей: либо с помощью функции biocLite (если библиотека находится под «опекой» Bioconductor), либо с помощью более универсальной функции install.packages (если библиотека находится в каком-нибудь ином, чем Bioconductor, репозитории, например, CRAN), - причем напрямую из рабочего пространства R среды. Однако в случае с библиотеками VariantTools и VariantAnnotation ни один из этих путей напрямую не сработает. Это связано с зависимостью указанных библиотек от библиотеки rtracklayer. В свою очередь библиотека rtracklayer зависит от библиотек XML и RCurl. Обычно такого рода зависимости при установке библиотек R функций учитываются автоматически, но не на этот раз. Поэтому для успешной установки библиотек VariantTools и VariantAnnotation первоначально необходимо инсталлировать на компьютере пакеты r-cran-xml и libcurl4-gnutls-dev:

sudo apt-get install r-cran-xml
sudo apt-get install libcurl4-gnutls-dev

Далее можно запустить R-консоль и установить критические зависимости библиотеки VariantTools:

packages = c(“XML", "RCurl", "rtracklayer")
source("http://bioconductor.org/biocLite.R")
biocLite(packages)

Теперь можно установить и саму библиотеку VariantTools (библиотека VariantAnnotation является зависимостью библиотеки VariantTools, так что она будет установлена автоматически):

source("http://bioconductor.org/biocLite.R")
biocLite("VariantTools")

Вот, собственно говоря, и все. Следует отметить, что описанная выше простота в установке является кажущейся, так как незнание первого этапа может попортить немало нервов. А так теперь можно воспользоваться широкими возможностями, предоставляемыми библиотеками VariantTools и VariantAnnotation, и не изобретать заново велосипед.
В завершение этого сообщения хочу выразить благодарность Петру Назарову (Genomics Research Laboratory of the Department of Oncology, Luxembourg Institute of HealthStrassen, Luxembourg) за консультативную помощь в работе с библиотеками VariantTools и VariantAnnotation.

четверг, 15 октября 2015 г.

Новости из мира РНК.

В майском номере журнала Nature появилось продолжение эпопеи с рекурсивным сплайсингом. На сей раз свой труд на суд научного сообщества предоставила интернациональная группа исследователей из Великобритании, Германии, Саудовской Аравии и Словении и озаглавили они сие дело «Recursive splicing in long vertebrate genes» (PubMed ID 25970246).
Под рекурсивным сплайсингом подразумевается двухстадийное удалении длинных (более 24 kb) интронов, содержащих скрытые рекурсивные сайты сплайсинга (см. рисунок). Во время первой стадии происходит распознавание и удаление 5’-концевой порции интрона, при этом используется канонический 5’-сайт сплайсинга и скрытый рекурсивный сайт сплайсинга в середине самого интрона. Структура скрытого рекурсивного сайта сплайсинга такова, что по завершении первой стадии происходит восстановление 5’-сайта сплайсинга. Как только 5’-сайт сплайсинга восстановился запускается второй раунд удаления интрона, на сей раз его 3’-коцевой части и уже с вовлечением канонического 3’-сайта сплайсинга.

Схема рекурсивного сплайсинга (источник: Sibley C. R. et al. Recursive splicing in long vertebrate genes. // Nature. 2015 May 21;521(7552):371-5. doi: 10.1038/nature14466).

Впервые такой механизм удаления длинных интронов был обнаружен у D. melanogaster, причем аж в 1998 году. Сделала это группа Javier A. Lopez из Университета Карнеги-Меллона (Питтсбург, США). Желающие могут познакомиться с этой пионерской работой, заглянув в 6-й номер 2-го тома журнала Molecular Cell за 1998 год (PubMed ID 9885566). После этого интерес к проблеме удаления очень длинных интронов едва теплился (мне известно лишь три работы, сделанные в период между 1998 годом и настоящим временем). Но вот он снова возник, и теперь в отношении генов млекопитающих. В частности, авторы выше упомянутой работы сосредоточились на генах, работающих в головном мозге. Среди таких генов есть немало с очень длинными интронами, которые как раз таки и удаляются по рекурсивному механизму.

среда, 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.