Язык R

и его применение в биоинформатике

Лекция 14

Анастасия Жарикова

1 декабря 2023

Еще раз про

установку пакетов

Репозитории

CRAN

  • Репозиторий CRAN
# установить пакет; делаем один раз
install.packages("tidyverse") 
# подгрузить библиотеку; делаем при загрузке R-сессии
library(tidyverse) 
# проверить установленную версию пакета
package_version("tidyverse")

Пакеты - полезные команды

remove.packages("packagename") - удалить пакет

update.packages() - обновить все пакеты

library() - список доступных пакетов

library("packagename") - загрузить установленный пакет в текущую R сессию

vignette("packagename") - посмотреть “красивый” мануал по пакету, есть не для всех пакетов

sessionInfo() - информация о загруженной сессии R

Пакеты

  • Установить пакет нужно один раз, но подгружать при каждом запуске рабочего сеанса

  • После уcтановки нужно только подгрузить пакет с помощью library("packagename")

Bioconductor

  • www.bioconductor.org
  • репозиторий R пакетов для решения задач биологической направленности:
    • software пакеты со структурами (классами) и функциями для анализа биологических данных
    • annotation пакеты, содержащие базы данных различных аннотаций
    • experiment пакеты с наборами данных, полученных в различных экспериметах
    • workflow пакеты с туториалами анализа биологических данных

Установка Bioconductor пакетов

С каждой новой версией R выходит и новая версия Bioconductor с новыми версиями пакетов. Поэтому определенной версии R соотвествует определенная версия пакета из Bioconductor. Посмотреть соответствие версий можно на сайте.

Пакеты из Bioconductor предпочтительнее устанавливать с помощью специальной системы управления пакетами - BiocManager.

# Установить BiocManager
install.packages("BiocManager")

# Установить нужный пакет из Bioconductor
BiocManager::install("AnnotationHub")

# Удалить уже ненужный пакет
remove.packages("AnnotationHub")

Версия Bioconductor

# Узнать версию Bioconductor
BiocManager::version()
[1] '3.18'
# Все ли Bioconductor пакеты актуальных версий
BiocManager::valid()

* sessionInfo()

R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=Russian_Russia.utf8  LC_CTYPE=Russian_Russia.utf8   
[3] LC_MONETARY=Russian_Russia.utf8 LC_NUMERIC=C                   
[5] LC_TIME=Russian_Russia.utf8    

time zone: Europe/Moscow
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] htmlwidgets_1.6.3   BiocManager_1.30.22 compiler_4.3.2     
 [4] fastmap_1.1.1       cli_3.6.1           tools_4.3.2        
 [7] htmltools_0.5.7     rstudioapi_0.15.0   yaml_2.3.7         
[10] rmarkdown_2.25      knitr_1.45          jsonlite_1.8.7     
[13] xfun_0.41           digest_0.6.33       rlang_1.1.2        
[16] evaluate_0.23      

Bioconductor version '3.18'

  * 1 packages out-of-date
  * 0 packages too new

create a valid installation with

  BiocManager::install("haven", update = TRUE, ask = FALSE, force = TRUE)

more details: BiocManager::valid()$too_new, BiocManager::valid()$out_of_date

Доступные пакеты

BiocManager::available(pattern = "Mmusculus")
 [1] "BSgenome.Mmusculus.UCSC.mm10"        "BSgenome.Mmusculus.UCSC.mm10.masked"
 [3] "BSgenome.Mmusculus.UCSC.mm39"        "BSgenome.Mmusculus.UCSC.mm8"        
 [5] "BSgenome.Mmusculus.UCSC.mm8.masked"  "BSgenome.Mmusculus.UCSC.mm9"        
 [7] "BSgenome.Mmusculus.UCSC.mm9.masked"  "EnsDb.Mmusculus.v75"                
 [9] "EnsDb.Mmusculus.v79"                 "PWMEnrich.Mmusculus.background"     
[11] "TxDb.Mmusculus.UCSC.mm10.ensGene"    "TxDb.Mmusculus.UCSC.mm10.knownGene" 
[13] "TxDb.Mmusculus.UCSC.mm39.knownGene"  "TxDb.Mmusculus.UCSC.mm39.refGene"   
[15] "TxDb.Mmusculus.UCSC.mm9.knownGene"

Установить из архива

Можно скачать архив с пакетом из любого репозитория

Сохранить архив

Установить пакет из архива

install.packages('C:/Users/User/Downloads/abc_2.1.zip', repos=NULL, type='source')`

Работа с геномными интервалами

Геномные интервалы

Локус или набор локусов в конкретном месте генома

Для точного определения геномного интервала нужно знать:

  • хромосому (имя соответствует таковому в референсном геноме)

  • начало интервала

  • конец интервала

  • цепь: “+” или “-” (может не быть)

  • версию референсного генома, относительно которого указаны координаты

Пример: chr1:2000000-3000000

Пример: chr1 2000000 3000000 -

Можно дополнительно задавать столбцы с метаданными интервалов

bedtools

Для работы с геномными интервалами вне R можно использовать

bedtools

bedops

tabix

Пакет GenomicRanges

manual

library(GenomicRanges)

GenomicRanges

Создать набор интервалов

gr <- GRanges(
    seqnames = Rle(c("chrK", "chr2"), c(1, 3)),
    ranges = IRanges(101:104, end = 111:114, names = head(letters, 4)),
    strand = Rle(strand(c("-", "+", "*")), c(1, 2, 1)),
    score = 1:4,
    GC = seq(1, 0, length=4))
gr
GRanges object with 4 ranges and 2 metadata columns:
    seqnames    ranges strand |     score        GC
       <Rle> <IRanges>  <Rle> | <integer> <numeric>
  a     chrK   101-111      - |         1  1.000000
  b     chr2   102-112      + |         2  0.666667
  c     chr2   103-113      + |         3  0.333333
  d     chr2   104-114      * |         4  0.000000
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

Без метаданных

granges(gr)
GRanges object with 4 ranges and 0 metadata columns:
    seqnames    ranges strand
       <Rle> <IRanges>  <Rle>
  a     chrK   101-111      -
  b     chr2   102-112      +
  c     chr2   103-113      +
  d     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

Только метаданные

mcols(gr)
DataFrame with 4 rows and 2 columns
      score        GC
  <integer> <numeric>
a         1  1.000000
b         2  0.666667
c         3  0.333333
d         4  0.000000

GenomicRanges

Подвыборка

gr[2:3, "GC"]
GRanges object with 2 ranges and 1 metadata column:
    seqnames    ranges strand |        GC
       <Rle> <IRanges>  <Rle> | <numeric>
  b     chr2   102-112      + |  0.666667
  c     chr2   103-113      + |  0.333333
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

Добавить метаданные

gr$len = width(gr)
gr
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chrK   101-111      - |         1  1.000000        11
  b     chr2   102-112      + |         2  0.666667        11
  c     chr2   103-113      + |         3  0.333333        11
  d     chr2   104-114      * |         4  0.000000        11
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

Информация о хромосомах

seqnames(gr)
factor-Rle of length 4 with 2 runs
  Lengths:    1    3
  Values : chrK chr2
Levels(2): chrK chr2

GenomicRanges

Подвыборка только канонических хромосом

library(GenomeInfoDb)
keepStandardChromosomes(gr, species = 'Homo_sapiens', pruning.mode="coarse")
GRanges object with 3 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  b     chr2   102-112      + |         2  0.666667        11
  c     chr2   103-113      + |         3  0.333333        11
  d     chr2   104-114      * |         4  0.000000        11
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

GenomicRanges

Подвыборка одной хромосомы

gr[seqnames(gr) == 'chr2',] -> grF
grF
GRanges object with 3 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  b     chr2   102-112      + |         2  0.666667        11
  c     chr2   103-113      + |         3  0.333333        11
  d     chr2   104-114      * |         4  0.000000        11
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Манипуляции с интервалами

  • создать набор интервалов без пересечений

  • пересечь 2 набора интервалов

  • изменить границы набора интервалов

  • исследовать покрытие интервалами

  • найти ближайший участок из другого набора интервалов

  • многое другое …

GenomicRanges

Набор интервалов без пересечений

reduce(grF) -> grR
grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Набор интервалов с учетом направления цепи

Имена хромосом

Имена хромосом в разметках могут быть представлены в разных стилях:

  • UCSC

  • NCBI

  • Ensembl

Разные консорциумы

“chr1” - “1” - “NC_000001.11”

Работы со стилями хромосом

library(GenomeInfoDb)

head(genomeStyles("Homo_sapiens"),5)
  circular auto   sex NCBI UCSC dbSNP Ensembl
1    FALSE TRUE FALSE    1 chr1   ch1       1
2    FALSE TRUE FALSE    2 chr2   ch2       2
3    FALSE TRUE FALSE    3 chr3   ch3       3
4    FALSE TRUE FALSE    4 chr4   ch4       4
5    FALSE TRUE FALSE    5 chr5   ch5       5
head(genomeStyles("Caenorhabditis_elegans"),5)
  circular auto   sex NCBI   UCSC Ensembl
1    FALSE TRUE FALSE    I   chrI       I
2    FALSE TRUE FALSE   II  chrII      II
3    FALSE TRUE FALSE  III chrIII     III
4    FALSE TRUE FALSE   IV  chrIV      IV
5    FALSE TRUE FALSE    V   chrV       V

Группы хромосом

  • auto

  • circular

  • sex

extractSeqlevelsByGroup(species="Homo_sapiens", style="NCBI", group="auto")
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "20" "21" "22"
extractSeqlevelsByGroup(species="Homo_sapiens", style="NCBI", group="circular")
[1] "MT"
extractSeqlevelsByGroup(species="Homo_sapiens", style="NCBI", group="sex")
[1] "X" "Y"

Проверить названия хромосом

seqlevelsStyle(paste0("chr",c(1:30)))
[1] "UCSC"
seqlevelsStyle(c("2L","2R","X","Xhet"))
[1] "NCBI"

Изменить стиль имен хромосом

mapSeqlevels(c("chrII", "chrIII", "chrM"), "NCBI")
 chrII chrIII   chrM 
  "II"  "III"   "MT"

GenomicRanges + стиль

Можно изменить стиль имен хромосом при необходимости

mapSeqlevels(seqlevels(grR), "NCBI") -> NCBI_name
grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
renameSeqlevels(grR,na.omit(NCBI_name))
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        2   102-113      +
  [2]        2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges + flank

Область до начала интервала

grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
flank(grR, 10, both = T)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2    92-111      +
  [2]     chr2    94-113      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges + promoters

Область вокруг начала интервала

grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
promoters(grR, upstream = 10, downstream = 5)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2    92-106      +
  [2]     chr2    94-108      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges + shift

Сдвинуть границы интервалов

grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
shift(grR, 10)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   112-123      +
  [2]     chr2   114-124      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges + resize

Изменить размеры интервалов, начиная с start

grR
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-113      +
  [2]     chr2   104-114      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
resize(grR, 100)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr2   102-201      +
  [2]     chr2   104-203      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Сравнение двух наборов

gr1 <- GRanges(seqnames = c("chr1", "chr2"), strand = "*",
               ranges = IRanges(start = c(2, 6), width = 3))
gr1
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1       2-4      *
  [2]     chr2       6-8      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths
gr2 <- GRanges(seqnames = c("chr1", "chr2"), strand = "*",
               ranges = IRanges(start = c(1, 3), width = 3))
gr2
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1       1-3      *
  [2]     chr2       3-5      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Пересечение

findOverlaps(gr1, gr2)
Hits object with 1 hit and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  -------
  queryLength: 2 / subjectLength: 2

Пересечение

pintersect(gr1, gr2)
GRanges object with 2 ranges and 1 metadata column:
      seqnames    ranges strand |       hit
         <Rle> <IRanges>  <Rle> | <logical>
  [1]     chr1       2-3      * |      TRUE
  [2]     chr2       6-5      * |      TRUE
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Непересечение

psetdiff(gr1, gr2)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1         4      *
  [2]     chr2       6-8      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Объединение

union(gr1, gr2)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1       1-4      *
  [2]     chr2       3-8      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

Ближайшие

nearest(gr1, gr2)
[1] 1 2

Расстояние до ближайшего

distanceToNearest(gr1, gr2)
Hits object with 2 hits and 1 metadata column:
      queryHits subjectHits |  distance
      <integer>   <integer> | <integer>
  [1]         1           1 |         0
  [2]         2           2 |         0
  -------
  queryLength: 2 / subjectLength: 2

Визуализация

Общая идея

Продемонстрировать результаты, привязанные к геномным координатам

  • отрисовать покрытие локуса чтениями

  • сравнить характеристики двух треков по высоте покрытия, количеству вариантов, …

  • посмотреть характеристики любого заданного локуса

  • отрисовать особенность генома (GC-состав, покрытие, генную разметку, …) полногеномно

  • …

karyoploteR

library(karyoploteR)

karyoploteR

kp <- plotKaryotype()

Основная идея

  • выбираем локус, хромосому или несколько хромосом

  • выбираем, что хотим отрисовать (работает с GRanges объектами)

  • рисовать можно сверху или снизу хромосомы

karyoploteR

Выбираем две хромосомы

kp <- plotKaryotype(genome = "hg38", chromosomes=c("chr10", "chr12"))

karyoploteR

Можно добавить опорные координаты

kp <- plotKaryotype(genome = "hg38", chromosomes=c("chr10", "chr12"))
kpAddBaseNumbers(kp)

karyoploteR

Расположить хромосомы в линию

kp <- plotKaryotype(genome = "hg38", chromosomes=c("chr10", "chr12"), plot.type=7)

karyoploteR

В зависимости от выбранного типа графика (plot.type) необходимые треки можно рисовать сверху и/или снизу хромосомы

kp <- plotKaryotype(plot.type=2, chromosomes=c("chr1", "chr2"))
kpDataBackground(kp)
kpDataBackground(kp, data.panel = 2)
kpAxis(kp, data.panel=1)
kpAxis(kp, data.panel=2)

karyoploteR покрытие

library(regioneR)
library(BSgenome.Hsapiens.UCSC.hg19)
regions <- createRandomRegions(nregions=10000, length.mean = 1e6, mask=NA, non.overlapping = FALSE)
kp <- plotKaryotype()
kpPlotDensity(kp, data=regions)

karyoploteR взаимодействия

starts <- sort(createRandomRegions(nregions = 25, length.sd = 8e6))
ends <- sort(createRandomRegions(nregions = 25, length.sd = 8e6))

kp <- plotKaryotype()
kpPlotRegions(kp, starts, r0=0, r1=0.5, col="#ff8d92")
kpPlotRegions(kp, ends, r0=0, r1=0.5, col="#8d9aff")
kpPlotLinks(kp, data=starts, data2=ends, col="#fac7ffaa", r0=0.5)

Gviz

  • Отдельно создаем необходимые треки

  • Отрисовываем в желательном типе и порядке

  • Работает с GRanges объектами

  • Большое количество референсных геномов

Gviz

Подготовка треков

library(GenomicRanges)
library(Gviz)
data(cpgIslands)
data(geneModels)
class(cpgIslands)
[1] "GRanges"
attr(,"package")
[1] "GenomicRanges"

Gviz

chr <- as.character(unique(seqnames(cpgIslands)))
gen <- genome(cpgIslands)
gtrack <- GenomeAxisTrack()
itrack <- IdeogramTrack(genome = gen, chromosome = chr)
atrack <- AnnotationTrack(cpgIslands, name = "CpG")

itrack миниатюра хромосомы, где красным указан выбранный локус

gtrack трек с детальной разметкой выбранного локуса

atrack визуализация CpG островков

Gviz

Отрисовываем все треки через список

plotTracks(list(itrack, gtrack, atrack))

Gviz

Добавляем генную разметку

data(geneModels)
grtrack <- GeneRegionTrack(geneModels, genome = gen,
                           chromosome = chr, name = "Gene Model")
plotTracks(list(itrack, gtrack, atrack, grtrack))

Gviz

Трек с данными (например, варианты)

lim <- c(26700000, 26750000)
coords <- sort(c(lim[1], 
                 sample(seq(from = lim[1], to = lim[2]), 99), 
                 lim[2]))
dat <- runif(100, min = -10, max = 10)
dtrack <- DataTrack(data = dat, start = coords[-length(coords)],
                    end = coords[-1], chromosome = chr, genome = gen, 
                    name = "Uniform")
plotTracks(list(itrack, gtrack, atrack, grtrack, dtrack), 
           from = lim[1], to = lim[2])

Gviz

Варианты отображения - гистограмма

plotTracks(list(itrack, gtrack, atrack, grtrack, dtrack), 
           from = lim[1], to = lim[2], type = "histogram")

Gviz

Варианты отображения

plotTracks(list(itrack, gtrack, atrack, grtrack, dtrack), 
           from = lim[1], to = lim[2], type = "smooth")

Gviz zoom

Удобный способ приближать данные …

plotTracks(list(itrack, gtrack, atrack, grtrack),
           from = 26700000, to = 26750000)

Gviz zoom

… вплоть до нуклеотидов

strack <- SequenceTrack(Hsapiens, chromosome = chr)
plotTracks(list(itrack, gtrack, atrack, grtrack, strack), 
           from = 26591822, to = 26591852, cex = 0.8)

circlize

circlize

library(circlize)
cell_cycle = data.frame(phase = factor(c("G1", "S", "G2", "M"), levels = c("G1", "S", "G2", "M")),
                        hour = c(11, 8, 4, 1))
color = c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3")
circos.par(start.degree = 90)
circos.initialize(cell_cycle$phase, xlim = cbind(rep(0, 4), cell_cycle$hour))
circos.track(ylim = c(0, 1), panel.fun = function(x, y) {
    circos.arrow(CELL_META$xlim[1], CELL_META$xlim[2], 
        arrow.head.width = CELL_META$yrange*0.8, arrow.head.length = cm_x(0.5),
        col = color[CELL_META$sector.numeric.index])
    circos.text(CELL_META$xcenter, CELL_META$ycenter, CELL_META$sector.index, 
        facing = "downward")
    circos.axis(h = 1, major.at = seq(0, round(CELL_META$xlim[2])), minor.ticks = 1,
        labels.cex = 0.6)
}, bg.border = NA, track.height = 0.3)

circlize

Конец!