Язык R

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

Лекция 14

Арина Никольская, Анастасия Жарикова

6 декабря 2024

Еще раз про

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

Репозитории

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.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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.4   BiocManager_1.30.25 compiler_4.3.3     
 [4] fastmap_1.2.0       cli_3.6.3           tools_4.3.3        
 [7] htmltools_0.5.8.1   rstudioapi_0.17.1   yaml_2.3.10        
[10] rmarkdown_2.29      knitr_1.49          jsonlite_1.8.9     
[13] xfun_0.49           digest_0.6.37       rlang_1.1.4        
[16] evaluate_1.0.1     

Bioconductor version '3.18'

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

create a valid installation with

  BiocManager::install(c(
    "bit", "cli", "colorspace", "cpp11", "ensembldb", "farver", "glue",
    "gtable", "Hmisc", "parallelly", "rlang", "RSQLite", "uwot", "withr"
  ), 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(rep("chr1", 4), "chrK")),
    ranges = IRanges(start=c(10,20,100,10,0), width = 50, 
                     names = head(letters,5)),
    strand = Rle(strand(c("+", "+", "-", "-", "*"))),
    score = 1:5,
    GC = seq(1, 0, length=5))
gr
GRanges object with 5 ranges and 2 metadata columns:
    seqnames    ranges strand |     score        GC
       <Rle> <IRanges>  <Rle> | <integer> <numeric>
  a     chr1     10-59      + |         1      1.00
  b     chr1     20-69      + |         2      0.75
  c     chr1   100-149      - |         3      0.50
  d     chr1     10-59      - |         4      0.25
  e     chrK      0-49      * |         5      0.00
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

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

granges(gr)
GRanges object with 5 ranges and 0 metadata columns:
    seqnames    ranges strand
       <Rle> <IRanges>  <Rle>
  a     chr1     10-59      +
  b     chr1     20-69      +
  c     chr1   100-149      -
  d     chr1     10-59      -
  e     chrK      0-49      *
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

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

mcols(gr)
DataFrame with 5 rows and 2 columns
      score        GC
  <integer> <numeric>
a         1      1.00
b         2      0.75
c         3      0.50
d         4      0.25
e         5      0.00

GenomicRanges

Подвыборка

gr[2:3, "GC"]
GRanges object with 2 ranges and 1 metadata column:
    seqnames    ranges strand |        GC
       <Rle> <IRanges>  <Rle> | <numeric>
  b     chr1     20-69      + |      0.75
  c     chr1   100-149      - |      0.50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

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

gr$len = width(gr)
gr
GRanges object with 5 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1     10-59      + |         1      1.00        50
  b     chr1     20-69      + |         2      0.75        50
  c     chr1   100-149      - |         3      0.50        50
  d     chr1     10-59      - |         4      0.25        50
  e     chrK      0-49      * |         5      0.00        50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges

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

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

GenomicRanges

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

gr[seqnames(gr) == 'chr1',] -> grF
grF
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1     10-59      + |         1      1.00        50
  b     chr1     20-69      + |         2      0.75        50
  c     chr1   100-149      - |         3      0.50        50
  d     chr1     10-59      - |         4      0.25        50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

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

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

Например, изменить границы интервала

  • С набором интервалов внутри одного объекта GRanges

  • С двумя наборами интервалами из разных GRanges

Например, объединить или пересечь интервалы

Отдельные интервалы

Посмотрим на интервалы

GenomicRanges — promoters

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

promoters(grF, upstream = 5, downstream = 10)
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1      5-19      + |         1      1.00        50
  b     chr1     15-29      + |         2      0.75        50
  c     chr1   140-154      - |         3      0.50        50
  d     chr1     50-64      - |         4      0.25        50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges — promoters

Область выделяется с учетом цепи

GenomicRanges — shift

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

shift(grF, 10)
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1     20-69      + |         1      1.00        50
  b     chr1     30-79      + |         2      0.75        50
  c     chr1   110-159      - |         3      0.50        50
  d     chr1     20-69      - |         4      0.25        50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges — shift

Сдвигаем вправо

GenomicRanges — resize

Изменить размеры интервалов

resize(grF, 30)
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1     10-39      + |         1      1.00        50
  b     chr1     20-49      + |         2      0.75        50
  c     chr1   120-149      - |         3      0.50        50
  d     chr1     30-59      - |         4      0.25        50
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges — resize

Изменяем размер от start

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

Вспомним, как выглядели исходные интервалы

GenomicRanges — reduce

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

reduce(grF) -> grR
grR
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1     10-69      +
  [2]     chr1     10-59      -
  [3]     chr1   100-149      -
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

GenomicRanges — reduce

Неперекрывающиеся интервалы с учетом цепи

Два набора интервалов

gr_a <- GRanges(seqnames = c("chr1", "chr1"), strand = c("+", "-"),
               ranges = IRanges(start = c(50, 75), width = 20, names = c("a1", "a2")))
gr_a
GRanges object with 2 ranges and 0 metadata columns:
     seqnames    ranges strand
        <Rle> <IRanges>  <Rle>
  a1     chr1     50-69      +
  a2     chr1     75-94      -
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
gr_b <- GRanges(seqnames = c("chr1", "chr1"), strand = c("+", "+"),
               ranges = IRanges(start = c(60, 100), width = 25, names = c("b1", "b2")))
gr_b
GRanges object with 2 ranges and 0 metadata columns:
     seqnames    ranges strand
        <Rle> <IRanges>  <Rle>
  b1     chr1     60-84      +
  b2     chr1   100-124      +
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

Пересекающиеся интервалы

Создадим интервалы с перекрытием

GenomicRanges — findOverlaps

Поиск пересечений

findOverlaps(gr_a, gr_b, minoverlap = 1)
Hits object with 1 hit and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  -------
  queryLength: 2 / subjectLength: 2
findOverlaps(gr_a, gr_b, minoverlap = 1, ignore.strand=TRUE)
Hits object with 2 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         2           1
  -------
  queryLength: 2 / subjectLength: 2

GenomicRanges — union

Объединение

union(gr_a, gr_b)
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1     50-84      +
  [2]     chr1   100-124      +
  [3]     chr1     75-94      -
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

GenomicRanges — union

Перекрывающиеся интервалы объединяются

GenomicRanges — intersect

Пересечение

intersect(gr_a, gr_b)
GRanges object with 1 range and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1     60-69      +
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

GenomicRanges — setdiff

Вычитание

setdiff(gr_a, gr_b)
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1     50-59      +
  [2]     chr1     75-94      -
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

GenomicRanges — setdiff

Ближайшие соседи

nearest(x = gr_b, #query
        subject = gr_a)
[1] 1 1
distanceToNearest(gr_a, gr_b)
Hits object with 1 hit and 1 metadata column:
      queryHits subjectHits |  distance
      <integer>   <integer> | <integer>
  [1]         1           1 |         0
  -------
  queryLength: 2 / subjectLength: 2

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

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

  • UCSC

  • NCBI

  • Ensembl

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

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

GenomicRanges — хромосомы

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

library(GenomeInfoDb)
keepStandardChromosomes(gr, species = 'Homo_sapiens', pruning.mode="coarse")
GRanges object with 4 ranges and 3 metadata columns:
    seqnames    ranges strand |     score        GC       len
       <Rle> <IRanges>  <Rle> | <integer> <numeric> <integer>
  a     chr1     10-59      + |         1      1.00        50
  b     chr1     20-69      + |         2      0.75        50
  c     chr1   100-149      - |         3      0.50        50
  d     chr1     10-59      - |         4      0.25        50
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

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

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
renameSeqlevels(grR,na.omit(NCBI_name))
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]        1     10-69      +
  [2]        1     10-59      -
  [3]        1   100-149      -
  -------
  seqinfo: 2 sequences from an unspecified genome; no seqlengths

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

Общая идея

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

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

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

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

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

  • …

karyoploteR

library(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=5)

karyoploteR

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

kp <- plotKaryotype(plot.type=2, chromosomes=c("chr1"))
kpDataBackground(kp)
kpDataBackground(kp, data.panel = 2)
kpAxis(kp, data.panel=1)
kpAxis(kp, data.panel=2)
kpText(kp, chr="chr1", x=60e6, y=0.5, labels="data.panel=1", data.panel = 1)
kpText(kp, chr="chr1", x=60e6, y=0.5, labels="data.panel=2", 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(plot.type=2, chromosomes = "chr21")
kpPlotDensity(kp, data=regions)
kpPlotRegions(kp, data=regions, data.panel=2)

Gviz

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

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

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

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

Gviz

Данные по CpG-островкам

library(GenomicRanges)
library(Gviz)
data(cpgIslands)
cpgIslands
GRanges object with 10 ranges and 0 metadata columns:
       seqnames            ranges strand
          <Rle>         <IRanges>  <Rle>
   [1]     chr7 26549019-26550183      *
   [2]     chr7 26564119-26564500      *
   [3]     chr7 26585667-26586158      *
   [4]     chr7 26591772-26593309      *
   [5]     chr7 26594192-26594570      *
   [6]     chr7 26623835-26624150      *
   [7]     chr7 26659284-26660352      *
   [8]     chr7 26721294-26721717      *
   [9]     chr7 26821518-26823297      *
  [10]     chr7 26991322-26991841      *
  -------
  seqinfo: 1 sequence from hg19 genome; no seqlengths

Gviz

Подготовка треков: аннотация

atrack <- AnnotationTrack(cpgIslands, name = "CpG")
plotTracks(atrack)

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

Gviz

Подготовка треков: координаты

gtrack <- GenomeAxisTrack()
plotTracks(list(gtrack, atrack))

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

Gviz

Подготовка треков: идеограмма

gen<-genome(cpgIslands) #hg19
chr <- as.character(unique(seqnames(cpgIslands))) #chr7
itrack <- IdeogramTrack(genome = gen, chromosome = chr)
plotTracks(list(itrack, gtrack, atrack))

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

Gviz

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

data(geneModels)
grtrack <- GeneRegionTrack(geneModels, genome = gen,
                           chromosome = chr, name = "Gene Model",
                           transcriptAnnotation = "symbol",
                           background.panel = "#FFFEDB",
                           background.title = "darkblue")
plotTracks(list(itrack, gtrack, atrack, grtrack))

Gviz — zoom

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

plotTracks(list(itrack, gtrack, atrack, grtrack),
           from = 26570000, to = 26605000, sizes=c(1,3,2,6))

Gviz zoom

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

strack <- SequenceTrack(Hsapiens, chromosome = chr)
plotTracks(list(itrack, gtrack, atrack, grtrack, strack), 
           from = 26591822, to = 26591852, cex = 0.8, sizes=c(1,3,2,6,1))

Gviz — покрытие

Принимает распространенные форматы файлов

bamFile <- system.file("extdata/test.bam", package = "Gviz")
itrack <- IdeogramTrack(genome = gen, chromosome = "chr1")
cov_track <- DataTrack(range = bamFile, genome = "hg19", type = "l", 
                     name = "Coverage", window = -1, chromosome = "chr1")
plotTracks(list(itrack, cov_track), from = 189990000, to = 190000000, sizes=c(1,2))

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 genomic

circos.initializeWithIdeogram(species = "mm10")
mouse_chromInfo = read.chromInfo(species = "mm10")$df # длины хромосом

# получим координаты середин хромосом
mouse_mid = data.frame(
    chr = paste0("chr", 1:19),
    mid = round((mouse_chromInfo[1:19, 2] + mouse_chromInfo[1:19, 3])/2)
)

circos.genomicLink(sample_n(mouse_mid ,5), 
                   sample_n(mouse_mid ,5), col = rand_color(5))

Конец!