1．1数据的包

有机体级(' org ')包包含中心标识符(例如，Entrez基因id)和其他标识符(例如，GenBank或Uniprot登录号，RefSeq id等)之间的映射。org包的名称总是这样的形式org。< Sp >。< id > .db(如。org.Sc.sgd.db),< Sp >是一种由两个字母组成的有机体的缩写。Sc为酿酒酵母),<标识>是描述中心标识符类型的缩写(小写)。sgd的基因标识符酿酒基因组数据库，或如对于Entrez基因ids)。中的“如何使用' .db '注释包”小插图AnnotationDbiPackage (org包只是" .db "注释包的一种)是一个关键引用。“.db”和大多数其他Bioconductor注释包每6个月更新一次。

注释包通常包含一个以包本身命名的对象。这些对象被统称为AnnotationDb对象，并命名了更具体的类OrgDb，ChipDb或TranscriptDb对象。可以应用于这些对象的方法包括关口()，键()，keytypes ()而且select ()．表中总结了检索注释的常用操作。

类别	函数	描述
发现	列()	列出可以返回的列的类型
	keytypes ()	列出可用作键的列
	键()	列出给定键类型的预期值
	select ()	检索匹配的注释键，keytype而且列
操作	setdiff ()，联盟()，相交()	集合操作
	复制()，独特的()	标记或删除重复项
	%, %，匹配()	找到匹配
	任何()，所有()	任何真正的？都是吗?
	合并()	结合两个不同的基于共享密钥
农庄*	成绩单()，外显子()，cd ()	特征(转录本，外显子，编码序列)为农庄．
	transcriptsBy ()，exonsBy ()	特征按基因、转录本等分组，如GRangesList．
	cdsBy ()

1．2互联网资源

下表是支持基于web查询的选择Bioconductor包的简短摘要。

包	描述
AnnotationHub	Ensembl, Encode, dbSNP, UCSC数据对象
biomaRt	合奏和其他注释
PSICQUIC	蛋白质相互作用
uniprot.ws	蛋白质注释
KEGGREST	KEGG通路
SRAdb	测序实验。
rtracklayer	基因组的痕迹。
GEOquery	数组和其他数据
ArrayExpress	数组和其他数据

1.3练习

练习1:本练习演示注释包的' select '接口的基本用法。

1. 安装及连接org.Hs.eg.db注释包;它包含的“符号映射”信息智人基于NCBI的“Entrez”标识符。

   库(org.Hs.eg.db)

2. 快速查看此包中的数据摘要

   org.Hs.eg.db

   # # OrgDb对象:# # | DBSCHEMAVERSION: 2.1 # # | Db型:OrgDb # # |支持包:AnnotationDbi # # | DBSCHEMA: HUMAN_DB # # |生物:智人# # |物种:人类# # | EGSOURCEDATE: 2019 - apr26 # # | EGSOURCENAME: Entrez基因# # | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA # # | CENTRALID:如# # | TAXID: 9606 # # | GOSOURCENAME:基因本体# # | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/ # # | GOSOURCEDATE: 2019 - apr24 # # | GOEGSOURCEDATE:2019 - apr26 # # | GOEGSOURCENAME: Entrez基因# # | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA # # | KEGGSOURCENAME: KEGG基因组# # | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes # # | KEGGSOURCEDATE: 2011 - mar15 # # | GPSOURCENAME: UCSC基因组生物信息学(智人)# # | GPSOURCEURL: # # | GPSOURCEDATE: 2018 - dec3 # # | ENSOURCEDATE: 2019 - apr08 # # | ENSOURCENAME:运用# # | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta # # | UPSOURCENAME: Uniprot # # | UPSOURCEURL:http://www.UniProt.org/ ## | UPSOURCEDATE:星期五4月26日20:12:58 2019

   ## ##请参阅:help('select')了解使用信息

3. 这个想法是存在的keytypes ()可以映射到不同的列()；键()可用于查看可用的密钥。浏览包以查看可用的信息类型，例如，

   keytypes (org.Hs.eg.db)

   ## [1] " accnum " " alias " " ensembl " " ensemblprot " ## [5] " ensembltrans " " entrezid " " enzyme " " evidence " ## [9] " evidenceall " " genename " " go " " goall " ## [9] " ipi " " map " " omim " " ontology " ## [17] " ontologyall " " path " " pfam " " pmid " ## [25] " prosite " " refseq " " symbol " " ucsckg " ## [25] " unigene " " uniprot "

   列(org.Hs.eg.db)

   ## [1] " accnum " " alias " " ensembl " " ensemblprot " ## [5] " ensembltrans " " entrezid " " enzyme " " evidence " ## [9] " evidenceall " " genename " " go " " goall " ## [9] " ipi " " map " " omim " " ontology " ## [17] " ontologyall " " path " " pfam " " pmid " ## [25] " prosite " " refseq " " symbol " " ucsckg " ## [25] " unigene " " uniprot "

   (钥匙(org.Hs.eg.db,“象征”))

   ##[1]“a1bg”“a2m”“a2mp1”“nat1”“nat2”“natp”

4. 从数据库中提取数据有两种基本方法org . *包-mapIds ()要在键和单列之间创建1:1映射，和select ()(通常有必要直接指定这个函数，以避免与dplyr发生冲突，例如AnnotationDbi: select ())。探索这些功能，例如，

   set.seed(123) egid <- sample(keys(org.Hs.eg.db)， 6) mapIds(org. hs .db, egid， "SYMBOL"， "ENTREZID")

   ## 'select()'返回键和列之间的1:1映射

   ##“rnu7-186p”“loc110120927”“irak1”“loc110485085”“pebp1p3”## 100129013 ##“actr3p2”

   注释dbi::select(org. hs . e.g. .db, egid, c("SYMBOL"， "ENSEMBL"， "GENENAME")， "ENTREZID")

   ## 'select()'返回键和列之间的1:1映射

   # # 1 # # ENTREZID象征运用106480363 rnu7 - 186 p < NA > # # 2 110120927 LOC110120927 < NA > # # 3 3654 IRAK1 ENSG00000184216 # # 4 110485085 LOC110485085 < NA > # # 5 100129570 PEBP1P3 < NA > # # 6 100129013 ACTR3P2 < NA > # # 1 # # GENENAME RNA,得以小型核186假基因# # 2 VISTA强化剂hs721 # # 3白介素- 1受体相关激酶1 # # 4 BRCA1P1基因间的重组地区# # 5磷脂酰乙醇胺结合蛋白1假基因3 # # 6 ACTR3假基因2

5. 有些键列映射是1:many，例如，Entrez ID“3812”映射到44个ensemble id。什么mapIds ()映射Entrez ID时返回“3812”到集合id ?使用附加参数multiVals = "CharacterList"进一步探索。将结果与返回的结果进行比较select ()．

   egid <- "3812" mapIds(org.Hs.eg.db, egid， "ENSEMBL"， "ENTREZID")

   ## 'select()'返回1:多个键和列之间的映射

   ## 3812 ## " ensg00000240403 "

   mapIds(org.Hs.eg.db, egid， "ENSEMBL"， "ENTREZID"， multiVals = "CharacterList")

   ## 'select()'返回1:多个键和列之间的映射

   ##长度为1的字符列表## [["3812"]]ENSG00000240403 ENSG00000278850…ENSG00000284381

   注释dbi::select(org. hs . e.g. .db, egid, c("SYMBOL"， "ENSEMBL")， multiVals = "CharacterList")

   ## 'select()'返回1:多个键和列之间的映射

   # # 1 # # ENTREZID象征运用3812 KIR3DL2 ENSG00000240403 # # 2 3812 KIR3DL2 ENSG00000278850 # # 3 3812 KIR3DL2 ENSG00000276424 # # 4 3812 KIR3DL2 ENSG00000275566 # # 5 3812 KIR3DL2 ENSG00000276004 # # 6 3812 KIR3DL2 ENSG00000278809 # # 7 3812 KIR3DL2 ENSG00000278726 # # 8 3812 KIR3DL2 ENSG00000278403 # # 9 3812 KIR3DL2 ENSG00000275629 # # 3812 KIR3DL2 ENSG00000276882 # # 11 3812 KIR3DL2 ENSG00000278710 # # 3812 KIR3DL2 ENSG00000275083 # # 13 3812 KIR3DL2 ENSG00000277181 # # 14 3812 KIR3DL2ENSG00000278656 # # 15 3812 KIR3DL2 ENSG00000274722 # # 16 3812 KIR3DL2 ENSG00000275416 # # 17 3812 KIR3DL2 ENSG00000278474 # # 18 3812 KIR3DL2 ENSG00000278707 # # 19 3812 KIR3DL2 ENSG00000275838 # # 3812 KIR3DL2 ENSG00000276739 # # 21 3812 KIR3DL2 ENSG00000277709 # # 22 3812 KIR3DL2 ENSG00000273911 # # 23 3812 KIR3DL2 ENSG00000278361 # # 24 3812 KIR3DL2 ENSG00000277982 # # 25 3812 KIR3DL2 ENSG00000275626 # # 26 3812 KIR3DL2 ENSG00000276357 # # 27 3812 KIR3DL2 ENSG00000275511 # # 28 3812 KIR3DL2ENSG00000273735 # # 29 3812 KIR3DL2 ENSG00000275262 # # 3812 KIR3DL2 ENSG00000278442 # # 31 3812 KIR3DL2 ENSG00000278758 # # 32 3812 KIR3DL2 ENSG00000284384 # # 33 3812 KIR3DL2 ENSG00000283975 # # 34 3812 KIR3DL2 ENSG00000284213 # # 35 3812 KIR3DL2 ENSG00000284295 # # 36 3812 KIR3DL2 ENSG00000284063 # # 37 3812 KIR3DL2 ENSG00000284466 # # 38 3812 KIR3DL2 ENSG00000284101 # # 39 3812 KIR3DL2 ENSG00000284046 # # 40 3812 KIR3DL2 ENSG00000284053 # # 41 3812 KIR3DL2 ENSG00000283951 # # 42 3812 KIR3DL2Ensg00000284192 ## 43 3812 kir3dl2 ensg00000284528 ## 44 3812 kir3dl2 ensg00000284381

6. 在返回值上使用tidyverse似乎经常是有用的mapIds ()而且select ()；探索这种用法

   library(tidyverse) egid <- keys(org. hs . exe .db) # all ENTREZIDs mapIds(org. hs .db, egid， "SYMBOL"， "ENTREZID") %>% as_tibble() %>% rownames_to_column("ENTREZID")

   ##警告:不鼓励在vector上调用' as_tibble() '，因为该行为在未来可能会改变。使用' tibble::enframe(name = NULL) '代替。##此警告在每个会话中显示一次

   ## A tibble: 61521 x 2 ## ENTREZID值##   ## 11 A1BG ## 2 2 A2M ## 3 3 A2MP1 ## 4 4 NAT1 ## 5 5 NAT2 ## 6 6 NATP ## 7 7 SERPINA3 ## 8 8 AADAC ## 9 9 AAMP ## 10 10 AANAT ## #

   注释dbi::select(org.Hs.eg.db, egid, c("SYMBOL"， "GO"， "GENENAME")， "ENTREZID") %>% as_tibble()

   ## #小猫咪:309097 x 6 # # ENTREZID象征本体GENENAME证据去# # <空空的> <空空的> <空空的> <空空的> <空空的> < >从而向# # 1 A1BG: 0002576助教BP alpha-1-B糖蛋白# 1 # 2 A1BG: 0003674 ND MF alpha-1-B糖蛋白# # 3 1 A1BG: 0005576 CC alpha-1-B糖蛋白# # 4注重科技进步1 A1BG: 0005576 IDA CC alpha-1-B糖蛋白5 # # 1 A1BG: 0005576助教CC alpha-1-B糖蛋白# # 6 1 A1BG: 0005615 CC注重科技进步alpha-1-B糖蛋白1 # # 7 A1BG: 0008150 ND BP alpha-1-B糖蛋白# # 8 1 A1BG: 0031093助教CCα -1- b糖蛋白## 9 1 A1BG GO:0034774 TAS CC α -1- b糖蛋白## 10 1 A1BG GO:0043312 TAS BP α -1- b糖蛋白## #

练习2:biomaRt．

本练习需要Internet访问

1. 浏览Biomart网站https://www.ensembl.org/biomart用于检索各种基因组注释。

   首先选择一个数据库(例如，' Ensembl Genes 92 ')、数据集(例如，' Human Genes (GRCh38.p12) ')、过滤器(例如，' GENE ' / ' Input external reference ' / ' GENE stable id '并输入' ENSG00000000003 ')、属性(默认是可以的)，然后按' Results '将Ensembl标识符映射到转录本标识符。

2. 安装(如有必要)并加载biomaRt包中。使用listMarts ()要查看可用的数据库，useMart ()选择您感兴趣的超市。

   库(biomaRt)头(listMarts ())

   ##生物艺术版本## 1集合_mart_集合基因97 ## 2集合_mart_小鼠小鼠品系97 ## 3集合_mart_snp集合变异97 ## 4集合_mart_funcgen集合调控97

   mart <- useMart("ENSEMBL_MART_ENSEMBL")

3. 使用listDatasets ()而且useDataset ()若要选择智人基因数据集。

   头(listDatasets(集市))

   # # # # 1 abrachyrhynchus_gene_ensembl粉足雁基因数据集描述(ASM259213v1) # # 2 acalliptera_gene_ensembl东部快乐基因(fAstCal1.2) # # 3 acarolinensis_gene_ensembl变色龙蜥蜴基因(AnoCar2.0) # # 4 acitrinellus_gene_ensembl Midas丽鱼科鱼的基因(Midas_v5) # # 5 ahaastii_gene_ensembl伟大发现猕猴桃基因(aptHaa1) # # 6 amelanoleuca_gene_ensembl熊猫基因(ailMel1) # # # #版本1 ASM259213v1 # # 2 fAstCal1.2 # # 3 AnoCar2.0 # # 4 Midas_v5 # # 5 aptHaa1 # # 6 ailMel1

   dataset <- usdataset ("hsapiens_gene_ensembl"， mart)

4. 使用listFilters ()查看可用的筛选器。筛选器是您要查询的数据类型。选择一个。

   头(listFilters(数据集)

   ##名称描述## 1 chromosome_name染色体/支架名称## 2 start开始## 3 end结束## 4 band_start Band开始## 5 band_end Band结束## 6 marker_start标记开始

   filters <- "ensembl_gene_id" #参见' listFilters() '

5. 使用listAttrbutes ()查看可用的属性。属性表示要检索的信息。选择一些!

   头(listAttributes(数据集)

   1 ensembl_gene_id基因稳定ID feature_page 2 ensembl_gene_id_version基因稳定ID版本feature_page 3 ensembl_transcript_id转录稳定ID feature_page 4 ensembl_transcript_id_version转录稳定ID版本feature_page 5 ensembl_peptide_id蛋白质稳定ID feature_page 6 ensembl_peptide_id_version蛋白质稳定ID版本feature_page

   attrs <- c("ensembl_gene_id"， "hgnc_symbol") #参见' listAttributes() '

6. 创建一个ensemble基因id的字符向量，组成并执行查询，将结果转换为tibble。

   ids <- c("ENSG00000000003"， "ENSG00000000005"， "ENSG00000000419"， "ENSG00000000457"， "ENSG00000000460"， "ENSG00000000938") tbl <- getBM(attrs, filters, ids, dataset) %>% as_tibble()

   ##发现缓存

   资源描述

   6 x 2 ## ensembl_gene_id hgnc_symbol ##   ## 1 ENSG00000000003 TSPAN6 ## 2 ENSG00000000005 TNMD ## 3 ENSG00000000419 DPM1 ## 4 ENSG00000000457 SCYL3 ## 5 ENSG00000000460 C1orf112 ## 6 ENSG00000000938 FGR

练习3:KEGGREST

本练习需要Internet访问

1. 探索KEGG网站https://www.genome.jp/kegg/KEGG是一个路径信息数据库。

2. 加载KEGGREST打包并发现可用的数据库

   库(KEGGREST) KEGGREST: listDatabases ()

   ##[13]“途径”“brite”“模块”“ko”“基因组”“vg”##[7]“ag”“化合物”“聚糖”“反应”“rclass”“酶”##[13]“疾病”“药物”“dgroup”“environ”“基因”“配体”##[19]“kegg”

3. 使用keggList ()在路径数据库中查询人类路径;将结果显示为tibble

   hsa_paths <- keggList("path "， "hsa") %>% tibble(path = names(.)， description = .) hsa_paths . sh <- keggList("path "， "hsa") %>% tibble(path = names(.)， description = .

   ## #小猫咪:##   ## 1路径:hsa00010糖酵解/糖异生-智人(人类)## 2路径:hsa00020柠檬酸循环(TCA循环)-智人(人类)## 3路径:hsa00030磷酸戊糖途径-智人(人类)## 4路径:hsa00040戊糖和葡萄糖醛酸盐相互转化-智人(人类)## 5路径:hsa00051果糖和甘露糖代谢-智人(人类)## 6路径:hsa00052半乳糖代谢-智人(人类)## 7路径:hsa00053抗坏血酸盐和aldarate代谢-智人(人类)## 8路径:hsa00061脂肪酸生物合成-智人(人类)## 9路径:hsa00062脂肪酸延伸-智人(人类)## 10路径:hsa00071脂肪酸降解-智人(人类)## #…多323行

4. 使用keggLink ()恢复每个通路中的基因。

   hsa_path_eg < - keggLink(“通路”,“保险公司”)% > %宠物猫(egid通路=,=子(“hsa:“,”“名称()))hsa_path_eg

   ##一个tibble: 30,967 x 2 ##路径egid ##   ## 1路径:hsa00010 10327 ## 2路径:hsa00010 124 ## 3路径:hsa00010 125 ## 4路径:hsa00010 126 ## 5路径:hsa00010 127 ## 6路径:hsa00010 128 ## 7路径:hsa00010 130 ## 8路径:hsa00010 131 ## 10路径:hsa00010 160287 ## #

   Hsa_path_eg %>% group_by(pathway) %>% summary (genes = list(egid))

   ## A tibble: 333 x 2 ##路径基因##   ## 1路径:hsa00010 [68]> ## 2路径:hsa00020  ## 4路径:hsa00040  ## 5路径:hsa00051  ## 7路径:hsa00053  ## 8路径:hsa00061  ## 9路径:hsa00062  ## #…with 323 more rows

5. 更新hsa_path_eg表中包含基因符号和集合id的信息org.Hs.eg.db包中。检索相关信息mapIds ()．如何处理映射到多个ensemble id的entrez基因id ?

   hsa_kegg_anno <- hsa_path_eg %>% mutate(symbol = mapIds(org.Hs.eg.db, egid， " symbol "， "ENTREZID")， ensembl = mapIds(org.Hs.eg.db, egid， " ensembl "， "ENTREZID"))

   ## 'select()'返回键和列之间的1:1映射

   ## 'select()'返回1:多个键和列之间的映射

6. 使用left_join ()控件的路径描述hsa_kegg_anno表格

   left_join (hsa_kegg_anno hsa_pathways)

   ##连接，by = "path "

   ## #小猫咪:30967 x 5 # #通路egid象征运用描述# # <空空的> <空空的> <空空的> <空空的> < >从而向# # 1路径:hsa00…10327年AKR1A1 ENSG000001糖酵解/糖质新生-坎…# # 2路径:hsa00…124年ADH1A ENSG000001糖酵解/糖质新生-坎…# # 3路径:hsa00……125 ADH1B ENSG000001糖酵解/糖质新生-坎……# # 4路径:hsa00…126年ADH1C ENSG000002糖酵解/糖质新生-坎…# # 5路径:hsa00…127年ADH4 ENSG000001糖酵解/糖质新生-坎…# # 6路径:128 ADH5 hsa00……130 ADH6 ENSG000001…糖酵解/糖异生- Hom…## 8路径:hsa00…130589 GALM ENSG000001.糖酵解/糖异生- Hom…## 9路径:hsa00…131 ADH7 ENSG000001.糖酵解/糖异生- Hom…## 10路径:hsa00…160287 LDHAL6A ENSG000001.糖酵解/糖异生- Hom…## #，共30,957行

2.1成绩单注释包

基因组中心包对于涉及基因组坐标的注释非常有用。例如，在感兴趣的区域发现编码序列的坐标，并从这些检索相应的DNA或蛋白质编码序列是直接的。容易使用基因组中心注释执行的操作类型的其他示例包括定义感兴趣的区域，以便在RNA-seq实验中计数对齐的读取，以及在ChIP-seq分析中检索感兴趣区域的DNA序列，例如，用于motif表征。

2．2rtracklayer

的rtracklayer软件包允许我们查询UCSC基因组浏览器，以及提供进口()而且export ()GFF、GTF和BED等常用注释文件格式的函数。下面的练习演示了的一些功能rtracklayer．

2．3练习

练习4:TxDb。*包

1. 安装及连接TxDb.Hsapiens.UCSC.hg38.knownGene包中。这包含了基因模型智人基于人类基因组的“hg38”构建，使用UCSC“knownGene”注释轨道中的基因注释;TxDb的最新版本和不同的注释轨道可用。查看包的摘要，并为便于输入创建别名

   图书馆(TxDb.Hsapiens.UCSC.hg38.knownGene) TxDb.Hsapiens.UCSC.hg38.knownGene

   # # TxDb对象:# # # Db型:TxDb支持包:# # # # # # GenomicFeatures数据来源:UCSC基因组:# # # # # # hg38生物:智人# # #分类ID: 9606 # # # UCSC的表:knownGene # # # UCSC的跟踪:GENCODE v29 # # #资源URL: http://genome.ucsc.edu/ # # #的基因类型ID: Entrez基因ID # # #完整数据集:是的# # # miRBase构建ID: NA # # # transcript_nrow: 226811 # # # exon_nrow: 647025 # # # cds_nrow: 301804 # # #创建的数据库:来自Bioconductor的GenomicFeatures包## #创建时间:2019-04-25 15:16:59 +0000(2019年4月25日星期四)## #基因组特性创建时的版本:1.35.11 ## RSQLite创建时的版本:2.1.1 ## DBSCHEMAVERSION: 1.2

   txdb <- txdb . hsapiens . ucsc .hg38. knowngene

2. 这个包的主要目的是提供基因组特征的基因组坐标，如外显子()，编码序列(cd ())，成绩单()而且基因()．例如，探索，

   Ex <- exons(txdb

   ## seqnames ranges strand | exon_id ##    |  ## [1] chr1 11869-12227 + | 1 ## [2] chr1 12010-12057 + b| 2 ## [3] chr1 12179-12227 + | 3 ## [4] chr1 12613-12697 + | 4 ## [5] chr1 12613-12721 + | 5 ## ... ... ... ... . ...[647021] chrUn_GL000220v1 155997-156149 + | 647021 ## [647022] chrUn_KI270442v1 380608-380726 + | 647022 ## [647023] chrUn_KI270442v1 217250-217401 - | 647023 ## [647024] chrUn_KI270744v1 51009-51114 - | 647024 ## [647025] chrUn_KI270750v1 148668-148843 + | 647025 ## ------- ## seqinfo:来自hg38基因组的595个序列(1个循环)

   库(ggplot2) qplot (log10(宽度(特异)))

   ## ' stat_bin() '使用' bin = 30 '。用' binwidth '选择更好的值。

   

   Ex [where .max(width(Ex))]

   ## seqnames ranges strand | exon_id ##    |  ## [1] chr12 102197585-102402596 + | 350923 ## ------- ## seqinfo:来自hg38基因组的595个序列(1个循环)

3. 提取所有基因，然后只保留“标准”染色体1:22,X, Y和m表()的seqnames ()来确定每条染色体上有多少基因。也可以用dplyr的方式来做;注意seqnames (gn)需要被胁迫as.factor ()．

   -基因(txdb)长度(Gn)

   ## [1] 26034

   性病< - paste0(“空空”,c(1:22说,“X”,“Y”,“M”))seqlevels (gn,修剪。Mode = "粗")<- STD长度(gn)

   ## [1] 26012

   seqlevels (gn)

   ##[1]“chr1”“chr2”“chr3”“chr4”“chr5”“chr6”“chr7”“chr8”“chr9”##[10]“chr10”“chr11”“chr12”“chr13”“chr14”“chr15”“chr16”“chr17”“chr18”##[19]“chr19”“chr20”“chr21”“chr22”“chrX”“chrry”“chrM”

   表(seqnames(gn))

   ## ## chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 ## 2625 1711 1497 1043 1206 1146 1231 922 1047 1039 1577 1315 ## chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr21 chr22 chrX chrY ## 592 869 892 1036 1407 421 1655 750 347 582 1038 62 ## chrM ## 2

   Tibble (chr = as.factor(seqnames(gn))) %>% group_by(chr) %>% summary (n = n())

   ## #一个tibble: 25 x 2 ## chr n ##   ## 1 chr1 2625 ## 2 chr2 1711 ## 3 chr3 1497 ## 4 chr4 1043 ## 5 chr5 1206 ## 6 chr6 1146 ## 7 chr7 1231 ## 8 chr8 922 ## 9 chr9 1047 ## 10 chr10 1039 ## ##

4. exonsBy ()按基因或转录子分组外显子;提取按基因分组的外显子。(挑战!)你能识别不同染色体上的外显子基因吗?标准染色体上有这些基因吗?

   exByGn <- exonsBy(txdb， "gene") ## trans <- length (unique(seqnames(exByGn))) table(trans)

   ## trans ## 1 2 34 5 6 7 8 9 10 13 16 ## 26175 629 111 27 34 39 79 72 2 17 11 ## 22 24 27 28 30 33 ## 111 111 3 1

   seqnames(exByGn[trans > 1])

   ## RleList长度为1020 ## $ ' 10000 ' ## factor-Rle长度为56，运行2次##长度:30 26 ##值:chr1 chr1_KI270763v1_alt ##级别(595):chr1 chr2 chr3…chrUn_KI270756v1 chrUn_KI270757v1 ## ## $ ' 100037417 ' ## factor-Rle长度为6，有2次运行##长度:3 3 ##值:chr22 chr22_KI270879v1_alt ##级别(595):chr1 chr2 chr3…chrUn_KI270756v1 chrUn_KI270757v1 ## ## $ ' 100049076 ' ## factor-Rle长度为20，有2次运行##长度:11 9 ##值:chr5 chr5_GL339449v2_alt ##级别(595):chr1 chr2 chr3…chrUn_KI270756v1 chrUn_KI270757v1 ## ## $ ' 100128260 ' ## factor-Rle长度为4，有2次运行##长度:2 2 ##值:chrX chrY ##级别(595):chr1 chr2 chr3…chrUn_KI270756v1 chrUn_KI270757v1 ## ## $ ' 100128292 ' ## factor-Rle长度为2，运行2次##长度:1 1 ##值:chr10 chr10_KI270825v1_alt ##级别(595):chr1 chr2 chr3…chrUn_KI270757v1 ## ##…## <1015个元素>

   ## std <- paste0("chr"， c(1:22， "X"， "Y"， "M")) unames <- unique(seqnames(exByGn[trans > 1])) transstd <- all(unames %in% std) unames[transstd]

   ##长度因子列表40 ## [["100128260"]]chrX chrY ## [["100132062"]] chr1 chr5 ## [["100302278"]] chr1 chr15 ## [["100359394"]] chrX chrY ## [["100422831"]] chr9 chr19 ## [["100422872"]] chr2 chr5 ## [["100422872"]] chr2 chr5 ## [["100422872"]] chrX chrY ## [["100422872"]] chr1 chr16 ## [["101929127"]] chr2 chr9 ## [["101954271"]] chr1 chr2 chr2 chr3 chr7 chr10 chr15 chr19 ##…## <30多个元素>

5. 前面的练习表明了这种基因“22947”在4号和10号染色体上都有外显子。了解更多关于这个基因的信息org.Hs.eg.db并在NCBI网站上搜索该基因符号。

   egid <- "22947" AnnotationDbi::select(org.Hs.eg.db, egid, c("SYMBOL"， "GENENAME")， "ENTREZID")

   ## 'select()'返回键和列之间的1:1映射

   DUX4L1双同源体4 like 1(假基因)

   url <- paste0("https://www.ncbi.nlm.nih.gov/gene/"， egid)

6. 注意TxDb。*包也支持keytypes ()，列(),select ()用于外显子、cds、转录本和基因标识符之间的映射。

练习5:BSgenome。*包

1. 安装(如有必要)并加载BSgenome.Hsapiens.UCSC.hg38包，包含hg38构建的整个序列智人．检查它的内容，并创建一个简单的别名。

   图书馆(BSgenome.Hsapiens.UCSC.hg38) BSgenome.Hsapiens.UCSC.hg38

   ##人类基因组:## #生物:智人(人类)## #提供者:UCSC ## #提供者版本:hg38 ## #发布日期:2013年12月## #发布名称:基因组参考联盟GRCh38 ## # 455序列:## # chr1 chr2 chr3 ## # chr4 chr5 chr6 ## # chr7 chr8 chr9 ## # chr10 chr11 chr12 ## # chr13 chr14 chr15 ## # ... ... ...chrUn_KI270751v1 chrUn_KI270751v1 chrUn_KI270752v1 ## #(使用'seqnames()'查看所有序列名称，使用'$'或'[[]' ##操作符访问给定序列)

   hg38 < - BSgenome.Hsapiens.UCSC.hg38

2. 基因组序列可以通过染色体检索，例如:hg38 [[" chr1 "]]，或按基因组范围，例如:getSeq (hg38农庄(chr1:1000000 - 2000000))．检索您最喜欢的DNA块并计算GC含量。

   dna <- getSeq(hg38, GRanges("chr1:1000000-2000000")) letterFrequency(dna， "GC"， as.prob=TRUE)

   [1，] 0.5728534

3. 使用org . *，TxDb。*,BSgenome。*包装来检索BRCA1外显子DNA序列。

   brca1_egid <- mapIds(org. hs . e.g. .db， "BRCA1"， "ENTREZID"， "SYMBOL")

   ## 'select()'返回键和列之间的1:1映射

   brca1_exons <- exonsBy(txdb， "gene")[[brca1_egid]] getSeq(hg38, brca1_exons)

   一个DNAStringSet实例，长度为80，宽度为1508Ctgcaaataaacttggtagcaaacacttcca ## [2] 998 caattgggcagatgtgtgaggcacctgtggtg…agactgtggctcaaaaaaaaaaaaaaaaaaaaaaa ## [3] 717 caattgggcagatgtgtgaggcacctgtggtg…AAGGAAACTTGAAACCTGGGCATGGTGGCTC # # [4] 240 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG…174 caattgggcagatgtgtgtgaggcacctgtggtg…accagagccacaggaccccaagaatgagctta ## ... ... ...## [76] 94 gagctcgctgagacttcctggacgggggacag…Gcgctcaggaggccttcaccctctgctctgg ## [77] 126 acagataaattaaaactgcgactgcgcggcgt…gcgctcaggaggccttcacccttgctctgg ## [78] 174 ttagcggtagccccttggtttccgcaacg…[79] 175 cttagcggtagccccttggtttccgcaac…gcgctcaggaggccttcaccctcctctgg ## [80] 120 aaagcgtgggaattacagataataaaactg

练习6

本练习使用注释资源从基因符号“BRCA1”到与该基因相关的每个转录本的基因组坐标，最后到转录本的DNA序列。这可以使用EnsDb与BSgenome包装，或与之结合TxDb，Homo.sapiens而且BSgenome包。这里我们将重点介绍前一种方法。

0. 使用AnnotationHub发现和检索当前的Ensembl注释(' EnsDb ')智人．

1. 使用cdsBy ()功能检索基因' BRCA1 '的所有编码序列的基因组坐标EnsDb.Hsapiens.v75包中。若要仅检索指定基因的数据，请提交GenenameFilter或函数的筛选公式/表达式过滤器参数。这避免了提取所有基因的编码区，这需要很长时间。

2. 在基因组坐标中可视化转录本Gviz包来构造GeneRegionTrack，并用plotTracks ()．

3. 使用Bsgenome.Hsapiens.UCSC.hg38包和extractTranscriptSeqs ()提取每个转录本的DNA序列。

解决方案

检索感兴趣基因的转录本分组编码序列，验证每个编码序列都是3的倍数。

library(EnsDb.Hsapiens.v75) edb <- EnsDb.Hsapiens.v75)v75 brca1cds <- cdsBy(edb, by = "tx"， filter = ~ genename == "BRCA1") class(brca1cds)

## [1] "CompressedGRangesList" ## attr(，"package") ## [1] "GenomicRanges"

长度(brca1cds)

29

Brca1cds[[1]] #外显子在CDS

14 # #农庄对象范围和3元数据列:# # seqnames范围链| gene_name exon_id # # < Rle > < IRanges > < Rle > | <人物> <人物> # # [1]17 41243452 - 41246659 | BRCA1 ENSE00003549219 # # 17 [2] 41242961 - 41243049 | BRCA1 ENSE00003547126 # # [3] 17 41234421 - 41234592 | BRCA1 ENSE00003527960 # # 17 [4] 41228505 - 41228631 | BRCA1 ENSE00003484886 # # [5] 17 41226348 - 41226538 | BRCA1 ENSE00003537850  ## ... ... ... ... . ... ...# # [10] 17 41209069 - 41209152 | BRCA1 ENSE00003458468 # # [11] 17 41203080 - 41203134 | BRCA1 ENSE00003477922 # # [12] 17 41201138 - 41201211 | BRCA1 ENSE00003628864 # # [13] 17 41199660 - 41199720 | BRCA1 ENSE00003687053 # # [14] 17 41197695 - 41197819 | BRCA1 ENSE00001312675 # # exon_rank # # <整数> # # 9 # # [1][2]10 # # 11 # # [3][4]12 # # 13 [5]  ## ... ...## [10] 18 ## [11] 19 ## [12] 20 ## [13] 21 ## [14] 22 ## ------- ## seqinfo:来自GRCh37基因组的1个序列

cdwidth <- width(brca1cds) #每个外显子的宽度all((sum(cdwidth) %% 3) == 0) # sum within cds，模量3

## [1] false

为什么有些成绩单的cd没有达到预期的长度?获取第一个cd大小错误的转录本的转录本ID，并在Ensembl基因组浏览器中查找该转录本(http://www.ensembl.org)。

Tx_cds_fail <- names(brca1cds)[(sum(cdwidth) %% 3) != 0] length(Tx_cds_fail)

## [1]

tx_cds_fail [1]

## [1] " enst00000412061 "

在记录的描述里说CDS 5 '不完整．因此，除了已知的蛋白质编码转录本外，Ensembl还为已知的无义介导的mRNA衰变靶标转录本或具有不完整CDS的转录本提供注释。然而，该等成绩单不会列于例如TxDb.Hsapiens.UCSC.hg38.knownGene包中。

接下来我们用Gviz(这个包有一个很棒的小插图，装饰图案(“Gviz”)）

使用ensembldb包中的函数以适合Gviz grt的##格式提取数据<- getGeneRegionTrackForGviz(edb, filter = ~genename == "BRCA1") plotTracks(list(GenomeAxisTrack()， GeneRegionTrack(grt)))

提取每个转录本的编码序列。EnsDb数据库提供来自Ensembl的注释，因此使用Ensembl风格的染色体名称(如“Y”)BSgenome包基于UCSC注释，该注释使用在每个染色体名称前加上“chr”的命名风格(例如“chrY”)。因此改变seqlevelsStyle从默认的UCSC染色体命名到Ensembl命名风格。

库(BSgenome.Hsapiens.UCSC.hg19)

## ##附加包:' bsgenome . hspapiens . ucsc .hg19'

bsgenome . hspiens . ucsc . bsgenome . hspiens . ucschg38': ## ## Hsapiens

基因组< - BSgenome.Hsapiens.UCSC。hg19 ##将seqlevelsStyle从UCSC更改为Ensembl seqlevelsStyle(genome) <- "Ensembl" tx_seq <- extractTranscriptSeqs(genome, brca1cds) tx_seq

##长度为29的DNAStringSet实例Atcccccacagccactactga enst00000309486 ## [2] 4875 atggatttatctgctcttcgc…Atcccccacagccactactga enst00000346315 ## [3] 2043 atggatttatctgctcttcgc…Atcccccacagccactactga enst00000351666 ## [4] 2166 atggatttatctgctcttcgc…Atcccccacagccactactga enst00000352993 ## [5] 4797 atggatttatctgctcttcgc…Atcccccacagccactactga enst00000354071 ## ... ... ...## [25] 1419 atggatttatctgctcttcgc…Aagaaggcaagcctccccaac enst00000494123 ## [26] 531 atgaatgtagaaaaggctgaa…Aagaaggcaagcctccccaac enst00000497488 ## [27] 522 atggatgctgagtttgtgtgt…Atcccccacagccactactga enst00000586385 ## [28] 1065 atgcacagttgctctgggagt…Atcccccacagccactactga enst00000591534 ## [29] 291 atgagtgacagcaagaaaacc…ATCCCCCACAGCCACTACTGA ENST00000591849

我们也可以检查CDS序列的成绩单与不完整的CDS。它们中的许多都没有开始密码子，因此表明CDS的5 '端是不完整的。

tx_seq [tx_cds_fail]

##一个长度为8的DNAStringSet实例## width seq names ## [1] 1312 GTTTGGATTCTGCAAAAAAGGCtgaagagataaagaaaaaaaaaa enst00000412061 ## [2] 958 ttcagcttgacacaggtttgga…Tcactccaaatcagtagagag enst00000473961 ## [3] 667 atggatttatctgctcttcgcg…gttttttctgcaaaaaagg enst00000476777 ## [4] 1867 atggatttatctgctcttcgcg…Atagttgttctagcagtgaag enst00000477152 ## [5] 1870 atggattttatctgctcttcgcg…Agtctattaaagaaagaaaaa enst00000478531 ## [6] 1495 gagctattgaaaatcatttgtg…Agtctattaaagaaagaaaaa enst00000484087 ## [7] 800 gagctattgaaaatcatttgtg…Taaagaaccaggagtggaaag enst00000487825 ## [8] 296 atggatttatctgctcttcgcg…AAAAGATGAAGTTTCTATCAT ENST00000489037

内含子坐标可以通过计算每个转录本覆盖的基因组范围(从第一个外显子的开始到最后一个外显子的结束)来识别，然后计算该范围与每个外显子覆盖的基因组坐标之间的(代数)集差

内含子<- psetdiff(unlist(range(brca1cds))， brca1cds)

检索内含子序列getSeq ()(这些都是不这样组装extractTranscriptSeqs ()外显子序列组装成成熟转录本);请注意，内含子以适当的受体和供体位点序列开始和结束。不幸的是，UCSC和Ensembl也使用了不同的基因组组装名称。的基因组名称更改内含子对象中的对象进行匹配基因组对象。

独特的(基因组(基因组))

##[1]“hg19”

基因组(内含子)

## 17 ##“GRCh37”

##更改内含子上的基因组名称以匹配BSgenome包基因组(内含子)<- c(' 17 ' = unique(基因组(基因组)))seq <- getSeq(基因组，内含子)names(seq)

# #[1]“ENST00000309486”“ENST00000346315”“ENST00000351666”“ENST00000352993”# #[5]“ENST00000354071”“ENST00000357654”“ENST00000412061”“ENST00000461221”# #[9]“ENST00000461574”“ENST00000461798”“ENST00000468300”“ENST00000470026”# #[13]“ENST00000471181”“ENST00000473961”“ENST00000476777”“ENST00000477152”# #[17]“ENST00000478531”“ENST00000484087”“ENST00000487825”“ENST00000489037”# #[21]“ENST00000491747”“ENST00000492859”“ENST00000493795”“ENST00000493919”# # [25]" enst00000494123 " " enst00000497488 " " enst00000586385 " " enst00000591534 " ## [29] " enst00000591849 "

seq[["ENST00000352993"]] # 20内含子

一个DNAStringSet实例，长度为20 ## width seq ## [1] 1840 GTAAGGTGCCTGCATGTACCTGTGCTATATGG…1417 gtaagtattgggtgccctgtcagagagggg…tttgaatgctctttttccttcctggggatccag ## [3] 1868 gtaagagcctgggagaaccccagttccagc…cagtgattttacatctaaatgtccattttag# # [4] 5934 gtaaagctccctccctcaagttgacaaaaatc…6197 gtaagtacttgatgttacaaactaaccagaga…atcctgatgggttgtttggtttctttcag ## ... ... ...## [16] 4241 gtaaaaccatttttttttttcttcttcttcttcttctt…Ttgcttgactgttctttaccatactgtttag ## [17] 606 gtaagtgttgaatatcccaagaatgacactca…caaacataatgttttttcccttgtattttacag ## [18] 1499 gtatataatttttggtaatgggtaggttggaa…tgagtgttttttctcaaacaatttaatttcag ## [19] 9192 gtaagttttttgaatgtgttatgtggctccattat…aaattgttttttttttttttttttttttttttttttttttttttttttttttttttttccccccctaccctgctag ## [20] 8237 gtaagtcagcacaagagtgtattaatttggga

练习7

本练习需要Internet访问

这里我们使用rtracklayer以检索ENCODE项目中跨细胞系鉴定的雌激素受体结合位点。我们关注的是一个特别有趣区域附近的结合位点。

1. 定义感兴趣的区域农庄实例具有适当的基因组坐标。我们的区域对应于一个特定基因的上游和下游10Mb。
2. 为UCSC基因组浏览器创建一个会话
3. 查询UCSC基因组浏览器编码雌激素受体ERalpha_一个转录标记;确定合适的基因轨道、表和转录因子需要生物学知识和检测工作。
4. 可视化结合位点的位置和它们的分数;注释感兴趣区域的中点。

解决方案

定义感兴趣的区域

library(GenomicRanges) roi <- GRanges("chr10"， IRanges(92106877, 112106876, names="ENSG00000099194"))

创建会话

library(rtracklayer) session <- browserSession()

在我们感兴趣的区域中查询UCSC以获得特定的轨道、表和转录因子

trackName <- " wgencoderegtbsclusteredv2 " tableName <- "wgEncodeRegTfbsClusteredV2" trFactor <- "ERalpha_a" ucscTable <- getTable(session, track=trackName, range=roi, table=tableName, name=trFactor))

可视化结果

plot(score ~ chromStart, ucscTable, pch="+") abline(v=start(roi) + (end(roi) - start(roi) + 1) / 2, col="blue")

5.1会话信息

sessionInfo ()

## R version 3.6.1 Patched (2019-07-16 r76845) ##平台:x86_64-apple-darwin17.7.0(64位)##运行在macOS High Sierra 10.13.6 ## ##矩阵产品:default ## BLAS: /Users/ma38727/bin/R-3-6-branch/lib/libRblas。/Users/ma38727/bin/R-3-6-branch/lib/libRlapack。dylib # # # #语言环境:# # [1]en_US.UTF-8 / en_US.UTF-8 en_US.UTF-8 / C / en_US.UTF-8 / en_US。UTF-8 ## ##附加的基本包:##[1]网格并行stats4统计图形grDevices utils ##[8]数据集方法基础## ##其他附加包:# # # # [1] BSgenome.Hsapiens.UCSC.hg19_1.4.0 [2] KEGGREST_1.25.0 # # [3] forcats_0.4.0 # # [4] stringr_1.4.0 # # [5] dplyr_0.8.3 # # [6] purrr_0.3.2 # # [7] readr_1.3.1 # # [8] tidyr_0.8.3 # # [9] tibble_2.1.3 # # [10] ggplot2_3.2.0 # # [11] tidyverse_1.2.1 # # [12] AnnotationHub_2.17.3 # # [13] BiocFileCache_1.9.1 # # [14] dbplyr_1.4.2 # # [15] Gviz_1.29.0 # # [16] biomaRt_2.41.6 # # [17] BSgenome.Hsapiens.UCSC.hg38_1.4.1 # # [18] BSgenome_1.53.0 # # [19] rtracklayer_1.45.1 # # [20] Biostrings_2.53.2 # # [21]XVector_0.25.0 # # [22] EnsDb.Hsapiens.v75_2.99.0 # # [23] ensembldb_2.9.2 # # [24] AnnotationFilter_1.9.0 # # [25] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.6 # # [26] GenomicFeatures_1.37.4 # # [27] GenomicRanges_1.37.14 # # [28] GenomeInfoDb_1.21.1 # # [29] org.Hs.eg.db_3.8.2 # # [30] AnnotationDbi_1.47.0 # # [31] IRanges_2.19.10 # # [32] S4Vectors_0.23.17 # # [33] Biobase_2.45.0 # # [34] BiocGenerics_0.31.5 # # [35] BiocStyle_2.13.2 # # # #通过加载一个名称空间(而不是附加):## [1] readxl_1.3.1 backports_1.1.4 ## [3] Hmisc_4.2-0 lazyeval_0.2.2 ## [5] splines_3.6.1 biocparallel_19.0 ## [5] checkmate_1.9.4 memoise_1.1.0 ## [13] cluster_2.1.0 modelr_0.1.4 ## [15] matrixStats_0.54.0 askpass_1.1 ## [17] prettyunits_1.0.2 colorspace_1.4-1 ## [19] blob_1.2.0 rvest_0.3.4 ## [23] rappdirs_0.3.1 haven_2.1.1 ## [25] RCurl_1.95-4.12 jsonlite_1.6 ## [27] zeallot_0.1.0[37] xtable_1. 4 progress_1.2.2 ## [39] htmlTable_1.13.1 foreign_0.8-71 ## [41] bit_1.1-14 Formula_1.2-3 ## [43] htmlwidgets_3 . 1.3 httr_1.4.0 ## [45] rcolorbrewer_1 . 1.1-2 acepack_1.4.1 ## [47] pkgconfig_2.0.2 XML_3.98-1.20 ## [51] tidyselect_0.2.5 labeling_0.3 ## [53] rlang_0.4.0 later_0.8.0 ## [55]munsell_0.5.0 cellranger_1.1.0 ## [57] tools_3.6.1 cli_1.1.0 ## [59] generics_0.0.2 RSQLite_2.1.1 ## [61] broom_0.5.2 evaluate_0.14 ## [63] yaml_2.2.0 knitr_1.23 ## [65] bit64_0.9-7 nlme_3.1-140 ## [67] mime_0.7 xml2_1.2.0 ## [69] compiler_3.6.1 rstudioapi_0.10 ## [71] curl_3.3 png_0.1-7 ## [73] interactiveDisplayBase_1.23.0 stringi_1.4.3 ## [75] lattice_0.20-38 ProtGenerics_1.17.2 ## [77] Matrix_1.2-17 vctrs_0.2.0 ## [79] pillar_1.4.2 BiocManager_1.30.4 ## [81] data.table_1.12.2 bitops_1.0-6 ## [83] httpuv_1.5.1 R6_2.4.0 ## [85] latticeExtra_0.6-28 bookdown_0.12 ## [87] promises_1.0.1 gridExtra_2.3 ## [89] codetools_0.2-16 dichromat_2.0-0 ## [91] assertthat_0.2.1 SummarizedExperiment_1.15.5 ## [93] openssl_1.4 withr_2.1.2 ## [95] GenomicAlignments_1.21.4 Rsamtools_2.1.3 ## [97] GenomeInfoDbData_1.2.1 hms_0.5.0 ## [99] rpart_4.1-15 rmarkdown_1.14 ## [101] biovizBase_1.33.0 shiny_1.3.2 ## [103] lubridate_1.7.4 base64enc_0.1-3

5.2确认

本教程中报告的研究得到了国家人类基因组研究所和国家卫生研究院国家癌症研究所的支持，奖励编号为U41HG004059和U24CA180996。

该项目已获得欧盟地平线2020研究和创新计划(资助协议编号633974)下欧洲研究委员会(ERC)的资助。