Домашнее задание №9

Задание 1

Воспользуйтесь пакетом с аннотацией генов человека org.Hs.eg.db. Найдите ENSEMBL и UNIPROT идентификаторы для генов STAT1 и SERPINE1. Что это за гены - посмотрите на GENENAME? Опишите ваши находки.

select(org.Hs.eg.db, c("STAT1", "SERPINE1"),  c("ENSEMBL", "UNIPROT", "GENENAME"), "SYMBOL")

## 'select()' returned 1:many mapping between keys and columns

##     SYMBOL         ENSEMBL    UNIPROT
## 1    STAT1 ENSG00000115415     P42224
## 2 SERPINE1 ENSG00000106366 A0A024QYT5
## 3 SERPINE1 ENSG00000106366     P05121
##                                             GENENAME
## 1 signal transducer and activator of transcription 1
## 2                           serpin family E member 1
## 3                           serpin family E member 1

После многих практикумов по биоинфе я научился делать выводы по табличкам со случайными циферками и буковками. * Для гена STAT1 найдена одна запись в UniProt. Это ген преобразователя сигнала и активатора транскрипции 1 (в душе не чаю, что это значит). * Для гена SERPINE1 найдено две записи в UniProt. Это ген член 1 семейства серпинов Е (не уверен, что слово член стоит на своём месте).

Задание 2

Воспользуйтесь пакетом с аннотацией генов человека org.Hs.eg.db. Переведите все содержащиеся в org.Hs.eg.db ENTREZ-идентификаторы генов в символьные названия SYMBOL и ENSEMBL идентификаторы с помощью функции AnnotationDbi::select().

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

genes <- keys(org.Hs.eg.db)
df <- select(org.Hs.eg.db, genes, c("SYMBOL", "ENSEMBL"), "ENTREZID")

## 'select()' returned 1:many mapping between keys and columns

df %>%
  filter(duplicated(df$ENTREZID)) %>% 
  dplyr::select("ENTREZID") %>% 
  base::unique() %>% 
  nrow()

## [1] 2419

Однозначная конвертация идентификаторов не удалась. У 2419 генов возникли дубликаты.

Задание 3

BioMart отлично справляется с конвертацией идентификаторов генов, например, человека, однако он может быть еще полезен при конвертации идентификатров между разными организмами (определении ортологов).

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

Для решения такой задачи удобнее всего создать два mart объекта - для мыши и человека. А затем с помощью функции getLDS() получить информацию из двух связанных объектов.

Если вы получаете следующую ошибку, попробуйте при создании mart объектов указать параметр host: например, host = "https://jul2023.archive.ensembl.org/".

Error: biomaRt has encountered an unexpected server error.
Consider trying one of the Ensembl mirrors (for more details look at ?useEnsembl)

mouse_genes <- c("Hdac2", "Timeless", "Prkcg", "Hlf", "Sin3a", "Ogt", "Id3", "Csnk1d", "Cdk5", "Ep300", "Cipc", "Relb")
mart_mus <- useMart("ENSEMBL_MART_ENSEMBL", "mmusculus_gene_ensembl", host = "https://dec2021.archive.ensembl.org")
mart_homo <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl", host = "https://dec2021.archive.ensembl.org")

getLDS(attributes = "external_gene_name",
       filters = "external_gene_name",
       values = mouse_genes,
       mart = mart_mus,
       attributesL = "external_gene_name",
       martL = mart_homo)

##    Gene.name Gene.name.1
## 1       Cipc        CIPC
## 2      Ep300       EP300
## 3        Ogt         OGT
## 4   Timeless    TIMELESS
## 5      Hdac2       HDAC2
## 6      Sin3a       SIN3A
## 7        Hlf         HLF
## 8       Relb        RELB
## 9       Cdk5        CDK5
## 10     Prkcg       PRKCG
## 11       Id3         ID3
## 12    Csnk1d      CSNK1D

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

P.S.

АХАХАХАХАХАХ, хост не работает, и другой тоже и третий. На форуме биокондуктора (https://support.bioconductor.org/p/129636/) сказано, что если зеркала не работают, то нужно написать вот такую функцию (она ещё и с вложенными циклами, хотя на парах нам говорили, что использовать циклы в R крайне не прилично). И знаете что? А ничто. Ссылка по которой должна скачиваться табличка с генами человека и мыши просто не работает. А знаете что работает? Википедия. В ней есть статьи про эти гены, более того, в ней даже есть все идентификаторы.

UPD: Если продолжить перебирать все хосты из списка (listEnsemblArchives()), то можно наткнуться на рабочий.

#функция с форума

library(dplyr)
mouse_human_genes = read.csv("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt",sep="\t") #нерабочая ссылка

convert_mouse_to_human <- function(gene_list){

  output = c()

  for(gene in gene_list){
    class_key = (mouse_human_genes %>% filter(Symbol == gene & Common.Organism.Name=="mouse, laboratory"))[['DB.Class.Key']]
    if(!identical(class_key, integer(0)) ){
      human_genes = (mouse_human_genes %>% filter(DB.Class.Key == class_key & Common.Organism.Name=="human"))[,"Symbol"]
      for(human_gene in human_genes){
        output = append(output,human_gene)
      }
    }
  }

  return (output)
}

#P.S. этот чанк не должен исполняться, это просто пример кода

Задание 4

Давайте обратимся к базе данных генов человека Ensembl версии 98 (это можно сделать как с помощью biomaRt, так и с помощью AnnotationHub - выбирайте по своему вкусу) и изучим характеристики всех генов разных биотипов (белок-кодирующие, lncRNA и т.д.). Получите значения длин генов и постройте графики распределения для разных биотипов.

#я не очень понимаю зачем, но этот пакет создал на компе папку размером 125MB (надеюсь R не планирует скачать геном человека)
ah <- AnnotationHub()
query(ah, c("EnsDb", "hsapiens", "98"))

## AnnotationHub with 1 record
## # snapshotDate(): 2023-10-20
## # names(): AH75011
## # $dataprovider: Ensembl
## # $species: Homo sapiens
## # $rdataclass: EnsDb
## # $rdatadateadded: 2019-05-02
## # $title: Ensembl 98 EnsDb for Homo sapiens
## # $description: Gene and protein annotations for Homo sapiens based on Ensem...
## # $taxonomyid: 9606
## # $genome: GRCh38
## # $sourcetype: ensembl
## # $sourceurl: http://www.ensembl.org
## # $sourcesize: NA
## # $tags: c("98", "AHEnsDbs", "Annotation", "EnsDb", "Ensembl", "Gene",
## #   "Protein", "Transcript") 
## # retrieve record with 'object[["AH75011"]]'

#а вот вообще вечность что-то скачивает (видимо действительно геном человека) UPD: оно внатуре сожрало гигабайт
hs_ens98 <- ah[["AH75011"]]

## downloading 1 resources

## retrieving 1 resource

## loading from cache

## require("ensembldb")

df1 <- select(hs_ens98, keys(hs_ens98), c("GENEBIOTYPE"))
df2 <- as.data.frame(ranges(genes(hs_ens98))) %>% 
  transmute(GENEID = names, len = width)
df <- left_join(df1, df2, by = "GENEID")
head(df)

##            GENEID    GENEBIOTYPE    len
## 1 ENSG00000000003 protein_coding  12884
## 2 ENSG00000000005 protein_coding  14950
## 3 ENSG00000000419 protein_coding  23689
## 4 ENSG00000000457 protein_coding  44637
## 5 ENSG00000000460 protein_coding 192074
## 6 ENSG00000000938 protein_coding  23122

#АХТУНГ логарифмическая шкала
df %>% 
  ggplot() +
  aes(x = GENEBIOTYPE, y = len, colour = GENEBIOTYPE) +
  geom_boxplot() +
  labs(title = "Распределение длин генов для разных биотипов", 
       x = "Биотип", y = "Длина") + 
  scale_y_log10(labels = scales::label_number()) +
  theme_bw() +
  theme(legend.position = "none",
        legend.title = element_text(hjust = 0.5),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, hjust = 1))

Домашнее задание №9

ФББ - R - 2023

Муравьёв Георгий

Задание 1

Задание 2

Задание 3

Задание 4