遺伝子アノテーションに含まれる同一座標を持つtranscriptを探す（あるいはEnsemblアノテーションに関する注意）

##背景
諸事情で古いEnsemblのGTFを使っています．
CufflinksとHTSeqの出力を比較した際，Cufflinksでは出力されない遺伝子IDがいくつか見つかりました．それらの遺伝子IDをみてみますと，別の遺伝子IDにアサインされているtranscriptが全く同一の遺伝子構造を持っていることがわかりました．
つまり，以下のことが分かりました（常識？）．
- Ensemblの（少なくとも）古いバージョンの遺伝子アノテーションには，本来マージされるべきtranscriptが別の物として登録されている
- そのような'重複した'transcriptのエントリをCufflinksとHTSeqでは異なる取り扱いをする
-- Cufflinksでは，非明示的にそれらをマージし，片方のIDだけ出力する
-- HTSeqでは，別々のtranscriptとして扱う


そこで，Rでそのようなtranscriptがどのくらいあるかを調べてみました．

## 準備

```r:アノテーションパッケージのインストール   
source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Mmusculus.UCSC.mm10.ensGene")
```


```r
# 遺伝子アノテーションのロード
library("TxDb.Mmusculus.UCSC.mm10.ensGene")
# オブジェクト名の確認
ls("package:TxDb.Mmusculus.UCSC.mm10.ensGene")
ensGene <- TxDb.Mmusculus.UCSC.mm10.ensGene
```

`findOverlap`がcircular seqではエラーがでるので，ミトコンドリア配列をオフにします．

```r
isActiveSeq(ensGene) <- c(chrM = F)
```

挙動確認用．

```r
length(transcripts(ensGene)) == length(exonsBy(ensGene, by = "tx"))
```

```
## [1] TRUE
```

## transcriptの座標から候補を絞り込む
まず，transcriptのstartとendが同一のtranscriptがあるものを探します．ここではエキソン構造が異なる場合も含まれます．

```r
tx <- transcripts(ensGene)
dup <- countOverlaps(tx, tx, type = "equal") > 1
tx.dup <- tx[dup]
hit <- findOverlaps(tx.dup, tx.dup, type = "equal")
hit <- hit[queryHits(hit) != subjectHits(hit), ]
```


```r
table(dup)
```

```
## dup
## FALSE  TRUE 
## 89589  5907
```

```r
head(tx.dup$tx_name)
```

```
## [1] "ENSMUST00000061280" "ENSMUST00000172291" "ENSMUST00000166384"
## [4] "ENSMUST00000168907" "ENSMUST00000119714" "ENSMUST00000163561"
```

```r
hit
```

```
## Hits of length 9246
## queryLength: 5907
## subjectLength: 5907
##      queryHits subjectHits 
##       <integer>   <integer> 
##  1            1           2 
##  2            2           1 
##  3            3           4 
##  4            4           3 
##  5            5           6 
##  6            6           5 
##  7            7           8 
##  8            8           7 
##  9            9          10 
##  ...        ...         ... 
##  9238      5899        5898 
##  9239      5900        5901 
##  9240      5901        5900 
##  9241      5902        5903 
##  9242      5903        5902 
##  9243      5904        5905 
##  9244      5905        5904 
##  9245      5906        5907 
##  9246      5907        5906
```


## エキソン構造を比較する
次にエキソン構造が一致しているものを取り出します．

```r
ex <- exonsBy(ensGene, by = "tx")
dup.ex <- countOverlaps(ex, ex, type = "within", maxgap = 0) > 1
table(dup.ex)
```

```
## dup.ex
## FALSE  TRUE 
## 79223 16273
```

```r
# ex からtranscript_idでGRangesを取り出せるように
names(ex) = tx$tx_name

# ちょっと高速化
library(multicore)

iden <- mclapply(1:length(hit), function(x) {
    a <- queryHits(hit)[x]
    b <- subjectHits(hit)[x]
    # startとendのデータフレームにした上で，同一かどうか比較
    identical(as.data.frame(ex[[tx.dup$tx_name[a]]])[, 2:3], as.data.frame(ex[[tx.dup$tx_name[b]]])[, 2:3])
})
table(unlist(iden))
```

```
## 
## FALSE  TRUE 
##  9202    44
```


## 結果

```r
hit.identical <- hit[unlist(iden), ]
hit.identical <- hit.identical[queryHits(hit.identical) < subjectHits(hit.identical), 
    ]
pair <- data.frame(a = tx.dup$tx_name[queryHits(hit.identical)], b = tx.dup$tx_name[subjectHits(hit.identical)])
head(pair)
```

```
##                    a                  b
## 1 ENSMUST00000055830 ENSMUST00000159201
## 2 ENSMUST00000126634 ENSMUST00000134468
## 3 ENSMUST00000056136 ENSMUST00000163374
## 4 ENSMUST00000036819 ENSMUST00000163908
## 5 ENSMUST00000097080 ENSMUST00000114295
## 6 ENSMUST00000051845 ENSMUST00000168574
```


試しに一つのペアを表示してググってみます．

```r
pair[3, ]
```

```
##                    a                  b
## 3 ENSMUST00000056136 ENSMUST00000163374
```

```r
ex[[as.character(pair[3, 1])]]
```

```
## GRanges with 2 ranges and 3 metadata columns:
##       seqnames                 ranges strand |   exon_id   exon_name
##          <Rle>              <IRanges>  <Rle> | <integer> <character>
##   [1]     chr1 [172341210, 172341451]      + |      8917        <NA>
##   [2]     chr1 [172368921, 172374085]      + |      8920        <NA>
##       exon_rank
##       <integer>
##   [1]         1
##   [2]         2
##   ---
##   seqlengths:
##                    chr1                 chr2 ...       chrUn_JH584304
##               195471971            182113224 ...               114452
```

```r
ex[[as.character(pair[3, 2])]]
```

```
## GRanges with 2 ranges and 3 metadata columns:
##       seqnames                 ranges strand |   exon_id   exon_name
##          <Rle>              <IRanges>  <Rle> | <integer> <character>
##   [1]     chr1 [172341210, 172341451]      + |      8917        <NA>
##   [2]     chr1 [172368921, 172374085]      + |      8920        <NA>
##       exon_rank
##       <integer>
##   [1]         1
##   [2]         2
##   ---
##   seqlengths:
##                    chr1                 chr2 ...       chrUn_JH584304
##               195471971            182113224 ...               114452
```

## 結果の検証
調べてみると，最新のEnsemblアノテーションではマージされているようです．

![ID history](http://cdn-ak.f.st-hatena.com/images/fotolife/y/yuifu/20130325/20130325001639.png)

Bioconductorのアノテーションパッケージは最新のものではないのでしょう．


では，最新のアノテーションではこのようなことはないのでしょうか．
調べてみましょう．

## EnsemblからTranscripDbオブジェクトを作成

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

```r

isActiveSeq(mmusculusEnsembl) <- c(MT = F)

tx.new <- transcripts(mmusculusEnsembl)
dup.new <- countOverlaps(tx.new, tx.new, type = "equal") > 1
tx.new.dup <- tx.new[dup.new]
hit <- findOverlaps(tx.new.dup, tx.new.dup, type = "equal")
hit <- hit[queryHits(hit) != subjectHits(hit), ]
table(dup.new)
```

```
## dup.new
## FALSE  TRUE 
## 87109  3810
```

```r
head(tx.new.dup$tx_name)
```

```
## [1] "ENSMUST00000166384" "ENSMUST00000168907" "ENSMUST00000142117"
## [4] "ENSMUST00000179089" "ENSMUST00000177608" "ENSMUST00000180062"
```

```r

ex.new <- exonsBy(mmusculusEnsembl, by = "tx")
dup.ex.new <- countOverlaps(ex.new, ex.new, type = "within", maxgap = 0) > 1
table(dup.ex.new)
```

```
## dup.ex.new
## FALSE  TRUE 
## 77370 13549
```

```r
names(ex.new) = tx.new$tx_name

library(multicore)

iden.new <- mclapply(1:length(hit), function(x) {
    a <- queryHits(hit)[x]
    b <- subjectHits(hit)[x]
    identical(as.data.frame(ex.new[[tx.new.dup$tx_name[a]]])[, 2:3], as.data.frame(ex.new[[tx.new.dup$tx_name[b]]])[, 
        2:3])
})
table(unlist(iden.new))
```

```
## 
## FALSE  TRUE 
##  6264    18
```

## 結果
```r
hit.identical.new <- hit[unlist(iden.new), ]
hit.identical.new <- hit.identical.new[queryHits(hit.identical.new) < subjectHits(hit.identical.new), 
    ]

pair.new <- data.frame(a = tx.new.dup$tx_name[queryHits(hit.identical.new)], 
    b = tx.new.dup$tx_name[subjectHits(hit.identical.new)])
head(pair.new)
```

```
##                    a                  b
## 1 ENSMUST00000142117 ENSMUST00000179089
## 2 ENSMUST00000055830 ENSMUST00000159201
## 3 ENSMUST00000126634 ENSMUST00000134468
## 4 ENSMUST00000171570 ENSMUST00000179386
## 5 ENSMUST00000066859 ENSMUST00000111976
## 6 ENSMUST00000110186 ENSMUST00000110188
```

## 結果の検証
ありました… Ensemblのウェブサイトでも確認できました．

![hoge](http://cdn-ak.f.st-hatena.com/images/fotolife/y/yuifu/20130325/20130325012557.png)


とりあえず，皆さまEnsemblアノテーションを使う場合にはご留意を．