R言語を使ったデータ解析

統計解析環境RのWindows版をこちらからダウンロードする。

マイクロアレイの発現データをこちら(data.txt)からダウンロードする。
(マイクロアレイのデータ解析の実習については、門田幸二氏の作成した例題から抜粋している。)

1. 平均発現量を求める

data <- read.table("data.txt", header=TRUE, row.names=1, sep="\t")       # ファイルの読み込み
apply(data, 1, mean)                                                     # 各行に対して、平均(mean)を計算
tmp <- cbind(rownames(data), data, apply(data, 1, mean))                 # 「行の名前(遺伝子名)」と「data」と「平均値」を列ベクトル単位で結合し、結果をtmpに格納
write.table(tmp, "out1.txt", sep = "\t", append=F, quote=F, row.names=F) # tmpの中身をout1.txtというファイル名で保存。

2. 最大発現量を示す組織を調べる

data <- read.table("data.txt", header=TRUE, row.names=1, sep="\t")       # ファイルの読み込み
max.col(data)                                                            # 最大発現量を示す組織の番号を表示
colnames(data)[max.col(data)]                                            # 最大発現量を示す組織名を表示
tmp <- cbind(rownames(data), max.col(data), colnames(data)[max.col(data)])
                                                                         # 「行の名前」と「最大発現量を示す組織の番号」と「組織名」を列ベクトル単位で結合し、
                                                                         # 結果をtmpに格納
write.table(tmp, "out2.txt", sep = "\t", append=F, quote=F, row.names=F) # tmpの中身をout2.txtというファイル名で保存。

3. 遺伝子間の階層型クラスタリングを行う

library(stats)  # statsパッケージをロードする
data <- read.table("data.txt", header=TRUE, row.names=1, sep="\t")  # ファイルの読み込み
data.dist <- dist(data, method = "euclidean")                       # 遺伝子間の距離を計算。デフォルトはユークリッド距離
data.hclust <- hclust(data.dist, method = "average")                # 平均距離法を適用
plot(data.hclust, main="hclust")                                    # 結果を表示
以下のようなクラスタリングの結果が出力される。

method="single"(最小距離法)、method="complete"(最大距離法)が適用される。これらのクラスタリングの結果を求め、上の結果と合わせて比較せよ。

4. 遺伝子間の非階層型クラスタリング(k-means法)を行う

install.packages("cclust")         # パッケージのインストール
library(cclust)
data <- read.table("data.txt", header=TRUE, row.names=1, sep="\t")   #ファイルの読み込み
k2 <- cclust(data, 2, 20, verbose=TRUE, method="kmeans")      # k = 2のとき
k3 <- cclust(data, 3, 20, verbose=TRUE, method="kmeans")      # k = 3のとき
.....
k2$cluster            # クラスタ分けの結果
k2$centers            # 各クラスタの中心座標

バイオ系データの解析

source("http://www.bioconductor.org/biocLite.R")         # パッケージのインストール
biocLite("GeneR")                                        # 生物配列の各種操作、データベースアクセス
library(GeneR)

ls("package:GeneR")

s <- "cgtagtagctagctagctagctagctag"
strComp(s)
strCompoSeq(s, wsize=1)
strCompoSeq(s, wsize=2)

biocLite("Biostrings")                                   # 生物配列の表記と検索アルゴリズム
library(Biostrings)

ls("package:Biostrings")

GENETIC_CODE
AMINO_ACID_CODE