Базовая аннотация сплайсинговых событий.

В этом разделе мы рассмотрим основные подходы в аннотации сплайсинговых событий, идентифицированных с помощью элайнера Subjunc. В качестве примера мы возьмем BED файл, использованный нами ранее, и загрузим его в рабочее пространство R среды с помощью функции import, как это было описано в разделе “Загрузка файла формата BED в рабочее пространство R среды”.

rm(list = ls())

library(rtracklayer)

splEvents = import(con = "D:/Transcriptome/SRR1145838.junction.bed",

                  format = "bed",

                  trackLine = FALSE,

                  genome = "hg38",

                  colnames = NULL,

                  which = NULL,

                  seqinfo = NULL)

После загрузки данных в созданном объекте splEvents мы преобразуем слот ranges – заменим геномные координаты блоков на координаты собственно сплайсинговых событий (так, что бы левая координата start соответствовала последнему нуклеотиду предшествующего экзона, а правая координата end – первому нуклеотиду следующего экзона). Кроме того, в преобразованном объекте мы оставим только два столбца метаданных – name и score, - так как остальные столбцы нам не понадобятся.

ranges(splEvents) = IRanges(start = ranges(splEvents)@start + width(unlist(splEvents$blocks)[c(TRUE, FALSE)]) - 1,

                                               end = end(splEvents) - width(unlist(splEvents$blocks)[c(TRUE, FALSE)]) + 1)

splEvents = splEvents[, 1:2]

splEvents

GRanges object with 245354 ranges and 2 metadata columns:

	seqnames		ranges	strand	|	name	score
	<Rle>		<IRanges>	<Rle>	|	<character>	<numeric>
[1]	chr1	[15038,	15796]	-	|	JUNC00202937	15
[2]	chr1	[15038,	186317]	-	|	JUNC00167075	30
[3]	chr1	[16310,	16607]	-	|	JUNC00140845	1
[4]	chr1	[16765,	16858]	-	|	JUNC00100248	13
[5]	chr1	[17055,	17233]	-	|	JUNC00214170	6

-------

seqinfo: 455 sequences from hg38 genome

Генная принадлежность идентифицированных сплайсинговых событий.

Теперь мы можем приступить к аннотации сплайсинговых событий. На первом этапе мы определим генную принадлежность сплайсинговых событий. Для этого в рабочее пространство R среды необходимо загрузить координаты известных генов (а точнее единиц транскрипции, без внешних промоторных регионов) человека, взятые из автоматической системы аннотирования Ensembl (81-й выпуск аннотаций, основанный на референс-сборке GRCh38.p3 генома человека, июнь 2015). Мы воспользуемся заранее подготовленным текстовым файлом, содержащим пять полей, разделенных табуляцией: идентификационный номер гена по Ensembl (Ensembl_Gene_ID), название хромосомы (Chromosome), начало (Start) и конец (End) гена в хромосоме, а также цепь (Strand), - и после загрузки данных преобразуем их в объект класса GRanges.

ensGenes = read.table("D:/Transcriptome/Ensembl genes, genomic coordinates, GRCh38.p3 human release 81.txt",

sep = "\t",

header = TRUE,

quote = "\"",

as.is = TRUE)

ensGenes = makeGRangesFromDataFrame(ensGenes, keep.extra.columns = TRUE)

ensGenes

GRanges object with 60411 ranges and 1 metadata column:

	seqnames		ranges	strand	|	Ensembl_Gene_ID
	<Rle>		<IRanges>	<Rle>	|	<character>
[1]	chrX	[100627109,	100639991]	-	|	ENSG00000000003
[2]	chrX	[100584802,	100599885]	+	|	ENSG00000000005
[3]	chr20	[50934867,	50958555]	-	|	ENSG00000000419
[4]	chr1	[169849631,	169894267]	-	|	ENSG00000000457
[5]	chr1	[169662007,	169854080]	+	|	ENSG00000000460

-------
seqinfo: 24 sequences from an unspecified genome; no seqlengths

Теперь можно сопоставить геномные интервалы в двух объектах. Для этого мы воспользуемся функцией findOverlaps из R/Bioconductor библиотеки IRanges (саму библиотеку дополнительно загружать не надо, так как она загружается в рабочее пространство R среды вместе с библиотекой rtracklayer). Эта функция имеет семь аргументов:

q  query – объект класса GRanges, содержащий список интервалов, которые необходимо картировать. В нашем случае это объект splEvents. В целом, на этот аргумент может быть подан любой другой объект близкого класса: Ranges, Views, RangesList, ViewsList или RangedData;

q  subject – объект класса GRanges, содержащий список интервалов, по которым будет идти картирование интервалов из query. В нашем случае это объект ensGenes. На этот аргумент может быть подан любой другой объект, но того же класса, что и объект, присвоенный аргументу query;

q  maxgap – максимальное расстояние, на котором могут находиться два интервала, что бы считаться перекрывающимися. Значение этого аргумента может быть только целочисленным положительным;

q  minoverlap – минимально необходимое перекрытие для двух интервалов, что бы классифицировать их как перекрывающиеся. Значение этого аргумента может быть только целочисленным положительным;

q  type – тип перекрытия. Этот аргумент может принимать одно из следующих значений: “any” (значение является значением по умолчанию и допускает любой вариант перекрытия), “start” (два интервала считаются перекрывающимися, если полностью совпадают их начала), “end” (два интервала считаются перекрывающимися, если полностью совпадают их концы), “within” (интервал из query должен полностью уместиться в интервал из subject) и “equal” (два интервала должны полностью совпасть);

q  select – метод выбора перекрывающегося интервала, если интервал из query перекрывается с несколькими интервалами из subject. Этот аргумент может принимать одно из следующих значений: “all” (значение является значением по умолчанию и возвращает все обнаруженные перекрытия), “first” (возвращает первый интервал из всего списка перекрывающихся интервалов), “last” (возвращает последний интервал из всего списка перекрывающихся интервалов) и “arbitrary” (возвращает произвольно выбранный интервал из всего списка перекрывающихся интервалов);

q  ignore.strand – логический аргумент, определяющий, будет ли учитываться цепь при поиске перекрытий (FALSE), или же принадлежность интервалов к той или иной цепи будет игнорироваться (TRUE). Аргумент добавлен к функции findOverlaps в R/Bioconductor библиотеке GenomicRanges, но отсутствует в библиотеке IRanges.

hits = findOverlaps(query = splEvents,

      subject = ensGenes,

      maxgap = 0L,

      minoverlap = 1L,

      type = "within",

      select = "all",

      ignore.strand = FALSE)

hits

Hits object with 229669 hits and 0 metadata columns:

	queryHits	subjectHits
	<integer>	<integer>
[1]	1	28054
[2]	3	28054
[3]	4	28054
[4]	5	28054
[5]	6	28054
...	...	...
[229665]	245311	15278
[229666]	245312	15278
[229667]	245313	15278
[229668]	245314	15278
[229669]	245343	35434

-------

queryLength: 245354
subjectLength: 60411

Полученный объект hits принадлежит к специальному классу объектов Hits, реализованному в R/Bioconductor библиотеке S4Vectors. Объекты такого класса хранят информацию о пересечениях (hits) двух векторов, один из которых называется запросом (query), а второй – предметом запроса (subject). В нашем объекте hits в столбце queryHits стоят номера строк (индексы) из объекта splEvents, а в столбце subjectHits номера строк из объекта ensGenes, с которыми строки из объекта splEvents перекрываются.
Следует отметить, что не все сплайсинговые события перекрылись с аннотированными в Ensembl генами: запрос length(unique(queryHits(hits))) дает цифру 217944, в то время как исходные данные включают 245354 записи. Кроме того, часть сплайсинговых событий может быть приписана не к одному, а к нескольким генам, что обусловлено перекрытием самих генов. В связи с этим мы сделаем две вещи: во-первых, создадим упорядоченный список events только тех сплайсинговых событий, которые перекрываются с известными генами, причем вместе с названиями этих генов; во-вторых, упростим объект splEvents, удалив из него те сплайсинговые события, которые не перекрываются с известными генами, отсортируем этот объект по полю name и добавим, воспользовавшись списком events, новое поле gene_ID.

events = data.frame(cbind(splEvents[queryHits(hits), ]$name,

        ensGenes[subjectHits(hits), ]$Ensembl_Gene_ID),

        stringsAsFactors = FALSE)

events = aggregate(events$X2 ~ events$X1, paste, collapse = ",", data = events)

colnames(events) = c("junction", "gene_ID")

splEvents = splEvents[unique(queryHits(hits)), ]

splEvents = splEvents[order(splEvents$name), ]

splEvents$gene_ID = events$gene_ID

splEvents = sort(splEvents)

splEvents

GRanges object with 217944 ranges and 3 metadata columns:

	seqnames		ranges	strand	|	name	score	gene_ID
	<Rle>		<IRanges>	<Rle>	|	<character>	<numeric>	<character>
[1]	chr1	[629963,	630048]	+	|	JUNC00197520	65	ENSG00000225630
[2]	chr1	[634018,	634319]	+	|	JUNC00241998	1	ENSG00000248527
[3]	chr1	[634550,	634610]	+	|	JUNC00102726	14	ENSG00000198744
[4]	chr1	[827775,	829003]	+	|	JUNC00075428	110	ENSG00000228794
[5]	chr1	[95151419,	95173816]	+	|	JUNC00245349	1	ENSG00000152078, ENSG00000271092

-------

seqinfo: 455 sequences from hg38 genome

При необходимости полученные результаты могут быть сохранены в виде, например, обычного текстового файла со столбцами, разделенными табуляцией:

res = as.data.frame(splEvents)

write.table(res,

     file = "D:/Transcriptome/Splicing events vs genes.txt",

     sep = "\t",

     quote = FALSE,

     col.names = TRUE,

     row.names = FALSE)

Идентифицированные сплайсинговые события vs известные сплайсинговые события.

Теперь мы выясним, какая часть из идентифицированных с помощью элайнера Subjunc сплайсинговых событий содержит уже описанные ранее стыки экзонов, а какая включает что-то новое, ранее не описанное. Что бы реализовать задумку мы воспользуемся транскрипционными моделями генов человека, построенными автоматической системой аннотирования Ensembl (81-й выпуск аннотаций, основанный на референс-сборке GRCh38.p3 генома человека, июнь 2015), и извлечем из них уже известные варианты стыков экзонов. Конечно, данные, обобщенные в Ensembl, не являются всеобъемлющими, но для наших целей вполне подходят.

Для этого нам понадобятся функциональные возможности двух дополнительных R/Bioconductor библиотек – GenomicFeatures и SGSeq. Возможности первой библиотеки предоставят нам доступ к локальной базе SQLite, содержащей транскрипционные модели генов человека (более подробно об этом написано в разделе “Организация и создание объектовкласса TxDb”). Возможности же библиотеки SGSeq позволят нам извлечь из транскрипционных моделей только экзонные стыки, причем полученные данные о стыках будут организованы в точно такой же объект класса GRanges, как и наш объект splEvents. Консолидированный R код представлен ниже.

library(GenomicFeatures)

library(SGSeq)

TxDb = loadDb("D:/Transcriptome/Ensembl.GRCh38.p3.release81.sqlite")

TxFe = convertToTxFeatures(TxDb)

refSpl = TxFe[TxFe@type == "J", ]

refSpl = as.data.frame(refSpl)

refSpl$seqnames = paste("chr", refSpl$seqnames, sep = "")

chr = c(paste("chr", c(1:22), sep = ""), "chrX", "chrY")

refSpl = refSpl[refSpl$seqnames %in% chr, ][, -6:-8]

rownames(refSpl) = NULL

refSpl = refSpl[!duplicated(refSpl), ]

refSpl = makeGRangesFromDataFrame(refSpl, keep.extra.columns = TRUE)

refSpl

GRanges object with 347559 ranges and 0 metadata columns:

	seqnames		ranges	strand
	<Rle>		<IRanges>	<Rle>
[1]	сhr1	[12057,	12179]	+
[2]	сhr1	[12057,	12613]	+
[3]	сhr1	[12697,	12975]	+
[4]	chr1	[12721,	13221]	+
[5]	chr1	[13052,	13221]	+

-------

seqinfo: 24 sequences from an unspecified genome; no seqlengths

Теперь мы можем сравнить два объекта – splEvents и refSpl, - и выяснить насколько полно они перекрываются. Для этого мы воспользуемся все той же функцией findOverlaps, но поменяем значение аргумента type на “equal” (равный).

hits = findOverlaps(query = splEvents,

      subject = refSpl,

      maxgap = 0L,

      minoverlap = 1L,

      type = "equal",

      select = "all",

      ignore.strand = FALSE)

hits

Hits object with 144931 hits and 0 metadata columns:

	queryHits	subjectHits
	<integer>	<integer>
[1]	4	46
[2]	6	48
[3]	7	49
[4]	8	50
[5]	13	55
...	...	...
[144927]	217934	345767
[144928]	217935	345768
[144929]	217938	345769
[144930]	217943	345771
[144931]	217944	347553

-------

queryLength: 217944

subjectLength: 347559

Видно, что из 217944 сплайсинговых событий, содержащихся в объекте splEvents, 144931 событие может быть классифицировано как уже известное, а остальные являются новыми. В действительности же, как это будет показано далее, большинство из идентифицированных “новых” событий подтверждается небольшим количеством ридов и, скорее всего, является шумом.

В завершение мы добавим к объекту splEvents новый столбец status, где будет содержаться информация о статусе идентифицированных сплайсинговых событий – являются ли они уже известными или же новыми.

splEvents$status = "known"

splEvents[-queryHits(hits), ]$status = "new"

При необходимости модифицированный объект splEvents может быть сохранен в виде обычного текстового файла, как это уже описывалось выше.

Количественная характеристика идентифицированных сплайсинговых событий.

Рассмотрим несколько вариантов количественной характеристики идентифицированных сплайсинговых событий. Первой из них будет дистанция сплайсинга (splicing distance) – фактически размер интрона, который удаляется во время объединения экзонов. Консолидированный код для решения этой задачи приведен ниже, а результат отображен на рисунке 1. На рисунке 1 отчетливо видны два пика – минорный пик так называемых минимальных интронов (с размером между 50 и 150 нуклеотидами) и основной пик стандартных интронов длиной более 150 нуклеотидов.

empSplDistance = log(width(splEvents) - 1, 10)

refSplDistance = log(width(refSpl) - 1, 10)

par(omi = c(0, 0.1, 0, 0), mai = c(1, 1, 0.1, 0.1), lwd = 3)

plot(density(empSplDistance), col = "red", ann = FALSE, cex.axis = 2, xlim = c(0, 8), ylim = c(0, 0.6))

lines(density(refSplDistance), col = "blue")

title(ylab = "Density", xlab = "Splicing distance, log10(nucleotides)", cex.lab = 2.5)

legend(x = "topright", legend = c("experimental events", "Ensembl"), col = c("red", "blue"), lwd = 3, bty = "n", cex = 2)

Рисунок 1. Распределение дистанций сплайсинга среди экспериментально идентифицированных и аннотированных в Ensembl сплайсинговых событий.

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

res = cbind(density(empSplDistance)$x, density(empSplDistance)$y,

                   density(refSplDistance)$x, density(refSplDistance)$y)

colnames(res) = c("Experimental events", "Density", "Ensembl events", "Density")

write.table(res, file = "Splicing distances.txt", sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)

Второй характеристикой будет количество ридов, подтверждающих сплайсинговое событие. В англоязычной литературе для обозначения этого параметра используется несколько терминов: number of supporting reads, sequencing depth или coverage. В нашем объекте splEvents такого рода информация содержится в поле score, так что нам остается только извлечь ее и привести в надлежащий вид.

library(edgeR)

splEvents$cpm = log(cpm(splEvents$score), 2)

par(omi = c(0, 0.05, 0, 0), mai = c(1, 1, 0.1, 0.1), lwd = 3)

plot(density(splEvents$cpm), col = "red", ann = FALSE, cex.axis = 2, xlim = c(-5, 10), ylim = c(0, 1.2))

lines(density(splEvents[splEvents$status == "new", ]$cpm), col = "blue")

lines(density(splEvents[splEvents$status == "known", ]$cpm), col = "darkgreen")

abline(v = 0, col = "black", lty = 2)

title(ylab = "Density", xlab = "Sequencing depth, log2(CPM)", cex.lab = 2.5)

legend(x = "topright", legend = c("all events", "new events", "known events"),

            col = c("red", "blue", "darkgreen"), lwd = 3, bty = "n", cex = 2)

Код, представленный выше, включает CPM-нормализацию и log₂-трансформацию первичных данных. Под CPM-нормализацией (counts per million) подразумевается пересчет количества ридов, подтверждающих сплайсинговые события, на 1 млн. ридов, содержащихся в RNA-Seq библиотеке. Для проведения такого рода нормализации мы воспользовались готовой функцией cpm из R/Bioconductor библиотеки edgeR.

На обобщающем рисунке 2 видно, что многие сплайсинговые события подтверждаются лишь небольшим количеством ридов. Особенно это заметно для “новых” событий, так что если выставить порог в 1 рид на миллион (пунктирная линия), то из этой группы событий сохранится лишь незначительная часть.

Рисунок 2. Подтвержденность ридами идентифицированных сплайсинговых событий.

Оценка силы сайтов сплайсинга.

Что бы оценить силу сайтов, задействованных в сплайсинге (score of splice sites, strength of splice sites), нам необходимо получить их сиквенс. Для этого нам понадобятся функциональные возможности трех R/Bioconductor библиотек: GenomicRanges (для создания объектов класса GRanges), Biostrings (необходима для извлечения сиквенсов и их записи в формате FASTA) и BSgenome.Hsapiens.UCSC.hg38 (как источник последовательности референс-сборки GRCh38 генома человека). Первые две уже загружены в рабочее пространство R среды, так как они связаны с библиотекой rtracklayer, а третью мы загрузим дополнительно. Далее последовательно будут созданы объекты fiveSS и threeSS с координатами 5’ и 3’ сайтов сплайсинга, из референс-сборки GRCh38 генома человека по полученным координатам будут извлечены сиквенсы этих сайтов сплайсинга и сохранены в виде отдельных файлов в FASTA формате. Консолидированный код, а также фрагмент одного из итоговых объектов, приведены ниже. Этот код, конечно, не претендует быть самым эффективным и универсальным, но может рассматриваться как один из вариантов решения поставленной задачи.

library(BSgenome.Hsapiens.UCSC.hg38)

splEvents_Plus = splEvents[splEvents@strand == "+", ]

splEvents_Minus = splEvents[splEvents@strand == "-", ]

fiveSS_Plus = as.data.frame(splEvents_Plus)

fiveSS_Plus$start = fiveSS_Plus$start - 2

fiveSS_Plus$end = fiveSS_Plus$start + 8

fiveSS_Minus = as.data.frame(splEvents_Minus)

fiveSS_Minus$end = fiveSS_Minus$end + 2

fiveSS_Minus$start = fiveSS_Minus$end - 8

fiveSS = c(makeGRangesFromDataFrame(fiveSS_Plus, keep.extra.columns = TRUE),

     makeGRangesFromDataFrame(fiveSS_Minus, keep.extra.columns = TRUE))

fiveSS_Seq = getSeq(Hsapiens, fiveSS)

names(fiveSS_Seq) = fiveSS$name

fiveSS_Seq = fiveSS_Seq[order(names(fiveSS_Seq)), ]

writeXStringSet(x = fiveSS_Seq, filepath = "Splicing events, 5' splice sites, FASTA.fasta", format = "fasta", width = 50)

threeSS_Plus = as.data.frame(splEvents_Plus)

threeSS_Plus$start = threeSS_Plus$end - 20

threeSS_Plus$end = threeSS_Plus$end + 2

threeSS_Minus = as.data.frame(splEvents_Minus)

threeSS_Minus$start = threeSS_Minus$start - 2

threeSS_Minus$end = threeSS_Minus$start + 22

threeSS = c(makeGRangesFromDataFrame(threeSS_Plus, keep.extra.columns = TRUE),

       makeGRangesFromDataFrame(threeSS_Minus, keep.extra.columns = TRUE))

threeSS_Seq = getSeq(Hsapiens, threeSS)

names(threeSS_Seq) = threeSS$name

threeSS_Seq = threeSS_Seq[order(names(threeSS_Seq)), ]

writeXStringSet(x = threeSS_Seq, filepath = "Splicing events, 3' splice sites, FASTA.fasta", format = "fasta", width = 50)

head(fiveSS_Seq)

A DNAStringSet instance of length 6

	width	seq	names
[1]	9	GTGGTGGGT	JUNC00000001
[2]	9	GAGGTATGT	JUNC00000002
[3]	9	AGGGTCAGT	JUNC00000003
[4]	9	CTGGTGAGT	JUNC00000004
[5]	9	AAGGTAAGG	JUNC00000005
[6]	9	AAGGTATCA	JUNC00000006

Теперь мы можем выяснить, насколько сильны интересующие нас сайты сплайсинга. Сделать это можно в рамках одной или нескольких математических моделей, связывающих нуклеотидную последовательность сайта сплайсинга и его силу: модель максимальной энтропии (maximum entropy model), модель максимально зависимой декомпозиции (maximum dependence decomposition model), марковская модель первой степени (first-order Markov model) и взвешенная модель (weight matrix model). Об этих моделях можно подробнее почитать, например, в работе Yeo G. и Burge C. B. "Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals”, опубликованной в 2004 году в 11-м томе журнала Journal of Computational Biology. 

Мы применим только первую из указанных моделей. Что бы не создавать ее заново, мы воспользуемся возможностями R/Bioconductor библиотеки spliceSites, разработанной профессором Wolfgang Kaisers из Дюссельдорфского университета им. Генриха Гейне.

source("https://bioconductor.org/biocLite.R")

biocLite("spliceSites")

library(spliceSites)

maxEnt = load.maxEnt()

Далее модель максимальной энтропии maxEnt будет последовательно подгоняться под наши эмпирические данные и будет рассчитываться интересующий нас параметр этой модели – сила сайтов сплайсинга. Для этого мы используем две функции из библиотеки spliceSites – score5 и score3, - работающие с 5’ и 3’ сайтами сплайсинга, соответственно. Эти функции имеют три одинаковых аргумента:

q  x – модель;

q  seq – нуклеотидная последовательность сайта сплайсинга;

q  pos – позиция в нуклеотидной последовательности, для которой будет вестись расчет.

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

scores = matrix(ncol = 2, nrow = length(splEvents))

for (i in 1:length(splEvents)){

            scores[i, 1] = round(score5(x = maxEnt, seq = toString(fiveSS_Seq[i]), pos = 3), digits = 2)

            scores[i, 2] = round(score3(x = maxEnt, seq = toString(threeSS_Seq[i]), pos = 20), digits = 2)

            if(i%%100 == 0){

                        cat("Junction", i, "of the", length(splEvents), "has been processed...\n")

                        flush.console()

            }

}

scores = data.frame(cbind(names(fiveSS_Seq), scores))

colnames(scores) = c("name", "fiveSS", "threeSS")

head(scores)

	name	fiveSS	threeSS
1	JUNC00000001	4.29	10.7
2	JUNC00000002	9.81	5.4
3	JUNC00000003	5.85	6.53
4	JUNC00000004	10.1	7.23
5	JUNC00000005	10.51	9.18
6	JUNC00000006	7.22	6.92

Полученные результаты мы можем добавить к объекту splEvents и, например, сохранить в виде обычного текстового файла со столбцами, разделенными табуляцией:

splEvents = splEvents[order(splEvents$name), ]

splEvents$fiveSS = as.numeric(levels(scores$fiveSS)[scores$fiveSS])

splEvents$threeSS = as.numeric(levels(scores$threeSS)[scores$threeSS])

res = as.data.frame(splEvents)

write.table(res,

     file = "D:/Transcriptome/Splicing events, strength of splice sites.txt",

     sep = "\t",

     quote = FALSE,

     col.names = TRUE,

     row.names = FALSE)

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