#这里是下载人类的已知基因的信息 (35.4 MB)
source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
然后查看我们下载的这个包里面所包含的信息
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
> txdb
TxDb object:
| Db type: TxDb
| Supporting package: GenomicFeatures
| Data source: UCSC
| Genome: hg19
| Organism: Homo sapiens
| UCSC Table: knownGene
| Resource URL: http://genome.ucsc.edu/
| Type of Gene ID: Entrez Gene ID
| Full dataset: yes
| miRBase build ID: GRCh37
| transcript_nrow: 82960
| exon_nrow: 289969
| cds_nrow: 237533
| Db created by: GenomicFeatures package from Bioconductor
| Creation time: 2014-09-26 11:16:12 -0700 (Fri, 26 Sep 2014)
| GenomicFeatures version at creation time: 1.17.17
| RSQLite version at creation time: 0.11.4
| DBSCHEMAVERSION: 1.0
可以看到这个UCSC数据库发布的经典的hg19版本基因组所有的基因信息,共有237533个CDS,共有289969个外显子
然后我们可以用几个简单的函数提取信息
>genes(txdb)
可以看到genes函数可以提取23056个基因信息,还是一个Granges对象
>exons(txdb)
而用exons函数可以提取这个txdb对象的exons信息,共289969个exon
同理还有 transcripts函数可以提取转录本信息,共82960个转录本
还有cds函数,提取到237533个cds信息
#我们可以提取外显子的GRanges对象来具体看看,也可以用genes,transcripts,cds等函数
exon_txdb=exons(txdb)
seqnames(exon_txdb)返回一个class 'Rle' [package "S4Vectors"] with 4 slots,有93个染色体信息,以及每条染色体上面有多少个外显子信息
ranges(exon_txdb)返回外显子的起始终止位点,长度,以及其它信息,也是一个对象class 'IRanges' [package "IRanges"] with 6 slots
还有很多函数
strand(exon_txdb)返回外显子的正负链信息,要么在正链要么在负链
mcols(exon_txdb)返回exon的id编号,1到27750个
seqlengths(exon_txdb)返回每条染色体的长度信息
names,length
GRanges对象还有很多其它类型的操作,非常好玩的,split,shift,resize,flank,reduce,gaps,disjoin,coverage
其它求交集并集和都可以用,union,intersect,setdiff,pintersect,psetdiff
然后我们再回头看看我们的重点,txdb对象
transcriptsBy(txdb,by="gene")
这个是把我们对象按照gene来对转录本分组,可以看到,分成了23459个元素的list,其中第一个基因有两个转录本,也有一些基因只有一个转录本,甚至有些基因会有非常多的转录本,也可以用exonsBy,cdsBy来对它进行处理
每一个元素都是GRangesList对象,就是前面的GRanges对象,
seqnames(x)
ranges(x)
strand(x)
mcols(x, use.names=FALSE)
elementMetadata(x)
values(x)
seqinfo(x)
seqlevels(x)
seqlengths(x)
isCircular(x)
genome(x)
seqnameStyle(x)
seqnames(x)
score(x)
还有很多关于它的介绍
http://web.mit.edu/r_v3.0.1/lib/R/library/GenomicRanges/html/GRangesList-class.html