Transcriptdbについて

Transcriptdbについて
========================================================
GenomicFeaturesパッケージのTranscriptDbオブジェクトでtranscriptのメタデータを管理することができる．TranscriptDbオブジェクトの裏にはSQLiteデータベースがある．

## TranscriotDbオブジェクトと周辺のクラスとの関係
![alt text](http://cdn-ak.f.st-hatena.com/images/fotolife/y/yuifu/20130321/20130321125011_original.png)


## パッケージの読み込み
```{r}
## source("http://bioconductor.org/biocLite.R")
## biocLite("GenomicFeatures")
library("GenomicFeatures")
```

## TranscriptDbオブジェクトのロード
元々あるデータからTranscriptDbオブジェクトをロードするには２つの方法があります．一つはloadDBを使う方法で，.sqliteファイルを直接読み込みます．もう一つはパッケージとして提供されているものを使う方法です．

```{r}
# Load transcript DB
samplefile <- system.file("extdata", "UCSC_knownGene_sample.sqlite",
                           package="GenomicFeatures")
txdb <- loadDb(samplefile)
txdb

# Load as annotation package
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
ls("package:TxDb.Hsapiens.UCSC.hg19.knownGene") #check object name
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene #shorthand (for convenience)
txdb
```

## TrancsriotDbオブジェクト
```{r}
is(txdb)

metadata(txdb)

seqinfo(txdb)
head(seqinfo(txdb)@seqlengths)
seqlengths(txdb)
seqlevels(txdb)
names(which(isCircular(txdb)))

# Use of as.list
txdump <- as.list(txdb)
is(txdump)
names(txdump)
## txdump
txdb1 <- do.call(makeTranscriptDb, txdump) stopifnot(identical(as.list(txdb1), txdump))

```

## データを染色体でフィルタする
```{r}
# See which chromosomes are currently available
isActiveSeq(txdb)

# List of chromosomes and/or contigs
seqlevels(txdb)  

# Filter by chromosome
isActiveSeq(txdb)[seqlevels(txdb)] <-  FALSE 
# same as: isActiveSeq(txdb) <- rep(FALSE, length(seqinfo(txdb)))

isActiveSeq(txdb) <- c("chr15"=TRUE)
```


## selectでSQLのようにデータをデータフレームとして抽出できる
```{r}
# Which keys are available
keytypes(txdb)

# Which values are available
cols(txdb)

keys <- c("100033416", "100033417", "100033420")

selected <- select(txdb, keys=keys, cols=c("TXNAME", "TXSTRAND"), keytype="GENEID")
selected
is(selected)
```


## 区間をGRangesオブジェクトとして取得する
```{r}
GR <- transcripts(txdb)
is(GR)
GR
GR[1:3]
length(GR)

# Filter by vals
GR <- transcripts(txdb, vals <- list(tx_chrom = "chr15", tx_strand = "+"))
length(GR)
unique(strand(GR))

# Extract Exons as ranges
EX <- exons(txdb)
EX[1:4]
length(EX)
length(GR)
```

## TSSの座標（幅１bpの区間）を抽出する
```{r}
(gr <- GR[3956:3959])
(extractedTSS <- resize(gr, width=1, fix="start"))
is(extractedTSS)
```


## ある区間と重なるものを区間として抽出する
```{r}
GR <- GRanges(seqnames = rep("chr15",2),
ranges = IRanges(start=c(500,70000000), end=c(10000,75000000)), strand = strand(rep("-",2)))
(extractedGR <- transcriptsByOverlaps(txdb, GR))
is(extractedGR)
```

## 区間を何らかのまとまりでまとめて（byで指定），GRangesListオブジェクトとして取得する
```{r}
GRList <- transcriptsBy(txdb, by = "gene")
length(GRList)
names(GRList)[10:13]
GRList[11:12]
is(GRList)

# 
GRList <- exonsBy(txdb, by = "tx")
length(GRList)
# Elements of GRangesList is accessed by [[]]
GRList[[12]]
tx_ids <- names(GRList)
head(select(txdb, keys=tx_ids, cols="TXNAME", keytype="TXID"))

# More specific
length(intronsByTranscript(txdb))
length(fiveUTRsByTranscript(txdb))
length(threeUTRsByTranscript(txdb))
```

## TranscriptDbで指定されている塩基配列をDNAStringSetとして取り出す
```{r}
library(BSgenome.Hsapiens.UCSC.hg19)
ls("package:BSgenome.Hsapiens.UCSC.hg19")
is(Hsapiens)

# Extract transcripts from genome
tx_seqs1 <- extractTranscriptsFromGenome(Hsapiens, txdb)
is(tx_seqs1)
length(tx_seqs1)
tx_seqs2 <- extractTranscriptsFromGenome(Hsapiens, exonsBy(txdb,by="tx"))
is(tx_seqs2)
length(tx_seqs2)
# CDS
cds_seqs <- extractTranscriptsFromGenome(Hsapiens, cdsBy(txdb, by="tx"))
translate(cds_seqs)

# Promoter sequences
## GRanges --> DNAStringSet
promoter_seqs1 <- getPromoterSeq(transcripts(txdb), Hsapiens, upstream=10, downstream=0)
is(promoter_seqs1)
length(promoter_seqs1)
## GRangesList --> DNAStringSetList
promoter_seqs2 <- getPromoterSeq(exonsBy(txdb,by="tx"), Hsapiens, upstream=10, downstream=0)
is(promoter_seqs2)
length(promoter_seqs2)
```

## TranscriptDbオブジェクトをつくる，パッケージ化する
```{r}
# Which tables in UCSC are available
supportedUCSCtables()

## mm9KG <- makeTranscriptDbFromUCSC(genome = "mm9", tablename = "knownGene")

head(getChromInfoFromUCSC("hg19"))

## mmusculusEnsembl <-
##    makeTranscriptDbFromBiomart(biomart = "ensembl",
##                               dataset = "mmusculus_gene_ensembl")
```

## TranscriptDbオブジェクトを保存する
```{r}
## saveDb(mm9KG, file="fileName.sqlite")
## mm9KG <- loadDb("fileName.sqlite")
```


## SessionInfo
```{r}
sessionInfo()
```


```{r timelog}
Sys.time()
```


    
## 参考
- [Making and Utilizing TranscriptDb Objects (pdf)](http://www.bioconductor.org/packages/2.11/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf)
- [Package ‘GenomicFeatures’](http://www.bioconductor.org/packages/2.11/bioc/manuals/GenomicFeatures/man/GenomicFeatures.pdf)