統計解析環境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