Установите и подгрузите пакет GenomicRanges.
#install.packages("BiocManager")
#BiocManager::install("GenomicRanges")
library(GenomicRanges)
## Загрузка требуемого пакета: stats4
## Загрузка требуемого пакета: BiocGenerics
## Загрузка требуемого пакета: generics
##
## Присоединяю пакет: 'generics'
## Следующие объекты скрыты от 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, union
##
## Присоединяю пакет: 'BiocGenerics'
## Следующие объекты скрыты от 'package:stats':
##
## IQR, mad, sd, var, xtabs
## Следующие объекты скрыты от 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
## unsplit, which.max, which.min
## Загрузка требуемого пакета: S4Vectors
##
## Присоединяю пакет: 'S4Vectors'
## Следующий объект скрыт от 'package:utils':
##
## findMatches
## Следующие объекты скрыты от 'package:base':
##
## expand.grid, I, unname
## Загрузка требуемого пакета: IRanges
##
## Присоединяю пакет: 'IRanges'
## Следующий объект скрыт от 'package:grDevices':
##
## windows
## Загрузка требуемого пакета: GenomeInfoDb
Откройте таблицу gencode_v37_hs_LARGE_exons.tsv и изучите её
содержимое. В таблице содержится информация об экзонах гена
LARGE1.
Какие типы транскриптов представлены в таблице (поле transcript_type) ?
На какой цепи закодированы интервалы? Конвертируйте интервалы в объект
GRanges.
data <- ...
Задавать поля объекта можно сразу столбцами таблицы.
exons_gr <- GRanges(
seqnames = Rle(...),
ranges = IRanges(start = ..., end = ...,
names = ...),
strand = Rle(...))
exons_gr
Это нам пригодится в самостоятельной работе.
Визуализируйте интервалы с помощью Gviz.
#BiocManager::install("Gviz")
library(Gviz)
chr <- ...
genome <- "hg38"
gtrack <- GenomeAxisTrack()
itrack <- IdeogramTrack(genome = ..., chromosome = ...)
atrack <- AnnotationTrack(..., name = "exons")
plotTracks(list(...),
from = 33*10**6, to = 34*10**6)
Заметим, что один и тот же экзон может несколько раз встречаться в разных изоформах транскрипта LARGE1. Чтобы упростить представление экзонов, применим функцию reduce из пакета GenomicRanges.
exons_gr_reduced <- ...
Визуализируем результат
atrack_reduced <- AnnotationTrack(..., name = "exons, reduced")
plotTracks(list(...),
from = 33*10**6, to = 34*10**6)
Теперь найдем координаты интронов. Для этого применим функцию gaps() — она находит участки, комплементарные имеющимся интервалам.
gaps_gr <- gaps(...)
gaps_gr
Функция gaps() также найдет участок за пределами первого экзона. Ограничим область поиска первым экзоном, указав start.
#Способ 1, можно придумать другой
first_exon_start <- min(start(exons_gr))
gaps_gr <- gaps(..., start = first_exon_start)
gaps_gr
Какие области длиннее, экзоны или интроны?
1) Добавьте новый столбец feature: для gaps_gr укажите “intron”, для
exons_gr_reduced “exon”.
2) Объедините два объекта GRanges в один.
3) Добавьте столбец с длиной интервалов и конвертируйте в формат,
подходящий для ggplot.
4) Средствами ggplot визуализируйте распределение длин интронов и
экзонов гена LARGE1.
gaps_gr$feature = ...
exons_gr$feature = ...
combined_gr <- c(,)
combined_gr$len <- ...
combined_df <- as.data.frame(..., row.names = NULL)
Отберите случайным образом по 50 строк из набора исходных интервалов два раза. Положите эти наборы в 2 разных GRanges объекта. Помните, что процедура выбора интервалв должна быть воспроизводима. Далее работайте с этими двумя наборами интервалов.
С помощью функции findOverlaps из пакета GenomicRanges узнайте, есть ли пересечения между наборами интервалов.
С помощью функции nearest из пакета GenomicRanges для интервалов из одного набора найдите ближайших соседей в втором наборе интервалов.
Визуализируйте интервалы с помощью Gviz, добавив их двумя разными треками друг под другом. Это тоже пригодится нам на самостоятельной.