出芽酵母のいちばん長い遺伝子

Galaxyなどを使っても良いのですが、Bioconductorのパッケージを使ってやってみます。

準備

モデル生物の中には、遺伝子を含む転写領域の情報がパッケージとして用意されているものがあります。
ここでは出芽酵母（Saccharomyces cerevisiae）の遺伝子領域リストを取得するために、BioconductorからTxDb.Scerevisiae.UCSC.sacCer3.sgdGeneパッケージをインストールします。

パッケージのインストール

source("http://bioconductor.org/biocLite.R")
biocLite("TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")

遺伝子領域のリストを得る

# パッケージの読み込み
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)

# 遺伝子名とその開始位置・終了位置などのリストをGRangeList形式で読み込む
sacCer3.gene.GRangesList <- transcriptsBy(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene,
                                          by="gene")

# データフレーム形式への変換
sacCer3.gene <- as.data.frame(sacCer3.gene.GRangesList)

今回使うパッケージでは、遺伝子領域のリストが表形式ではないGRangesList形式で提供されていますが、
表形式であるデータフレームに変換したほうが解析の見通しがつきやすいでしょう。
head()関数を使ってデータフレームの一部を確認できます。

result

> head(sacCer3.gene)
  element seqnames start   end width strand tx_id tx_name
1   Q0010     chrM  3952  4338   387      +  6665   Q0010
2   Q0010     chrM  4254  4415   162      +  6666   Q0017
3   Q0032     chrM 11667 11957   291      +  6667   Q0032
4   Q0055     chrM 13818 16322  2505      +  6668   Q0050
5   Q0055     chrM 13818 18830  5013      +  6669   Q0055
6   Q0055     chrM 13818 19996  6179      +  6670   Q0060

##一番長い遺伝子を探す

データフレーム内の目的行を抽出

# 遺伝子の長さが最大の行番号を取得
sacCer3.gene.maxRow <- which.max(sacCer3.gene$end-sacCer3.gene$start)

# 上で取得した行番号をつかってデータフレームの一部にアクセス
sacCer3.gene[sacCer3.gene.maxRow,]

上のコードを実行すると、

result

> sacCer3.gene[sacCer3.gene.maxRow,]
     element seqnames  start    end width strand tx_id tx_name
4209 YLR106C   chrXII 349006 363738 14733      -  4390 YLR106C

こんな結果が出てきます。YLR106CというのはMDN1の別名です。

いやぁ、Rって本当にいいもんですね

「勝つか負けるかはもう問題ではない。日本の国民を生かすか殺すかなのです」

そもそもなんで酵母の最長遺伝子長なんてものが知りたくなったのかはこちらに書きました。