- Open R and change the working (or current) directory to "20110711" you created.
- Prepare a file (template.txt) including template pattern for detecting heart-specific genes. In this practice, you don't have to prepare the file because you already have the file (template.txt) in "20110711" folder.
- Reading two files (data.txt and template.txt) with the following commands (It's OK by just doing "copy and paste".).
#------ Start ------
in_f1 <- "data.txt" #we use the gene expression data file (data.txt) as in_f1
in_f2 <- "template.txt" #we use the template file (template.txt) as in_f2
library(genefilter) #Open "genefilter" library to use a function, genefinder, in the library for pattern-matching
data_tmp <- read.table(in_f1, header=TRUE, row.names=1)#Reading in_f1 (data.txt) and we use the data as "data_tmp"
data <- as.matrix(data_tmp) #Exchanging data format (from "data.frame") to "matrix" for object "data_tmp" using the function "as.matrix" and we use the "data_tmp" after formatting change as "data"
template_tmp <- read.table(in_f2) #Reading "in_f2" (template.txt) and we use the template data as "template_tmp"
template <- template_tmp[,2] #Extracting the information of the second column in the "template_tmp" (i.e., template.txt) and we will use the information as "template"
template #Checking the information of "template". The gray color commands indicate "not essential."
#------ End ------
Extracting the top 10 genes having highest similarities to the template pattern.
#------ Start ------
gene_num <- 10 #We use the number of top-ranked genes (i.e., 10) as gene_num
dim(data) #Showing the numbers of rows and coloumns of an object "data"
tmp <- rbind(data, template) #Combining two objects ("data" and "template") by row. The "template" is added below "data" and we use the combined one as "tmp". You can refer the detailed usage of the function, rbind, by typing "?rbind".
dim(tmp) #Showing the numbers of rows and coloumns of an object "tmp". You can see the row number of "tmp" becomes the row number of "data" plus 1. This is because one row is added for "data".
template_posi <- which(rownames(tmp) == "template")#We use a row number whose row name is "template" as template_posi. Practically, the "template_posi" can be replaced by "nrow(tmp)".
closeg <- genefinder(tmp, template_posi, gene_num, method="correlation")#Extracting the top gene_num genes from the gene expression matrix "tmp" having highest similarities to the template pattern defined as that in the "template_posi"th line (or row) in the matrix "tmp". We use the results as "closeg"
closeg[[1]]$indices #Row numbers for the top-ranked genes
rownames(tmp[closeg[[1]]$indices,]) #Gene IDs for the genes. The underline part indicates those expression data.
closeg[[1]]$dists #Similarity values defined as "1 - correlation coefficient". These values can be changed if the corresponding parameter in the "genefinder" function is changed. For example, method="euclidean". For details, type "?genefinder".
#------ End ------
- Save the expression data for the top-ranked genes as the output file name of "top10_heart.txt".
#------ Start ------
out_f1 <- "top10_heart.txt" #we use the output filename (top10_heart.txt) as out_f1
tmp2 <- cbind(rownames(tmp[closeg[[1]]$indices,]), tmp[closeg[[1]]$indices,])#Two objects (Gene IDs and the expression data for the "gene_num" top-ranked genes) are combined by column and we use the combined one as "tmp2".
write.table(tmp2, out_f1, sep = "\t", quote=F, row.names=F)
#You may think "write.table(tmp[closeg[[1]]$indices,], "top10_heart_hoge.txt", sep = "\t", row.names=T)" is more simple compared to the above two commands.
#It is of course OK. But... Open the "top10_heart_hoge.txt" using EXCEL and compare it with "top10_heart.txt". I believe you can understand why I recommend above somewhat complicated commands.
#------ End ------
If you would like to save only the IDs, type as follows:
#------ Start ------
out_f2 <- "top10_heart_ID.txt" #we use the output filename (top10_heart_ID.txt) as out_f2
write.table(rownames(tmp[closeg[[1]]$indices,]), out_f2, quote=F, row.names=F, col.names=F)
#------ End ------
- Making pseudo-color image
#------ Start ------
heatmap(as.matrix(tmp[closeg[[1]]$indices,]), Rowv =NA, Colv=NA, scale="row", col = heat.colors(20), main="Best 10", xlab="Tissue", ylab="Clone ID", margin=c(8,10))
#------ End ------
- Rを起動し、「ファイル」−「ディレクトリの変更」を選択し、デスクトップ上に作成した20110711に作業ディレクトリを変更する。
- 目的組織(心臓)で1、それ以外の組織で0とした「テンプレートパターン」を含んだファイル(template.txt)を予め作成。
注意点:22,283行(遺伝子)×36列(組織)からなる遺伝子発現行列を含んだファイル(data.txt)中の組織の並びと、template.txt中の組織の並びは同じにしておく必要があります。
- 22,283遺伝子の中から心臓で特異的に発現している遺伝子の上位10個をリストアップします。
まず、以下をコピー&ペーストして、パターンマッチングを行って候補遺伝子群を抽出したいファイル(data.txt)と「テンプレートパターン」を含んだファイル(template.txt)を読み込みます。
#------ ここから ------
in_f1 <- "data.txt" #解析したい発現データファイル(data.txt)をin_f1に格納
in_f2 <- "template.txt" #解析したいテンプレートファイル(template.txt)をin_f2に格納
library(genefilter) #パッケージの読み込み
data_tmp <- read.table(in_f1, header=TRUE, row.names=1)#in_f1を読み込んでdata_tmpに格納(これをそのまま使うわけではないので"_tmp"を付加している)
data <- as.matrix(data_tmp) #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。(read.tableで読み込んで得られたdata_tmpのデータの型は"データフレーム"なので、そのままではこの場合は使用不可なのでやる必要があります)
template_tmp <- read.table(in_f2) #in_f2を読み込んでtemplateに格納
template <- template_tmp[,2] #テンプレート情報(2列目)のみ抽出し、templateに格納
template #テンプレート情報の確認
#------ ここまで ------
genefinder()関数を使ってテンプレートパターンに最も似た発現パターンをもつ上位10個を抽出します。この際、類似度は1−相関係数(method=”correlation”)で評価することにします。
#------ ここから ------
gene_num <- 10 #上位10個という情報をgene_numに格納
dim(data) #行列dataの行数と列数を表示
tmp <- rbind(data, template) #templateというテンプレートパターンを行列dataの後に追加
dim(tmp) #行列dataにテンプレートパターンが追加された行列tmpの行数と列数を表示
template_posi <- which(rownames(tmp) == "template")#行のラベル情報が"template"の行番号をtemplate_posiとして格納
closeg <- genefinder(tmp, template_posi, gene_num, method="correlation")#行列tmpの中で、template_posi行目の遺伝子発現プロファイルと近い上位gene_num個の情報をclosegに格納
closeg[[1]]$indices #上位gene_num個の行番号を表示
rownames(tmp[closeg[[1]]$indices,]) #上位gene_num個の遺伝子ID(AFFX_ID)を表示。underline部分が上位gene_num個の遺伝子発現データに相当する。
closeg[[1]]$dists #上位gene_num個の類似度を表示
#------ ここまで ------
参考:一連の計算はExcelでもできます。計算を行ったExcelファイルはここに置いてあります。
- 上位10個の遺伝子発現プロファイル情報を取り出し、ファイル(ファイル名:top10_heart.txt)に保存します:
#------ ここから ------
out_f1 <- "top10_heart.txt" #出力ファイル名(top10_heart.txt)をout_f1に格納
tmp2 <- cbind(rownames(tmp[closeg[[1]]$indices,]), tmp[closeg[[1]]$indices,])#上位10個の遺伝子ID(AFFX_ID)情報とその遺伝子発現情報を列方向で結合(column方向でbind)し、結果をtmp2に格納
write.table(tmp2, out_f1, sep = "\t", quote=F, row.names=F) #結果をファイル(ファイル名:top10_heart.txt)に保存。
#------ ここまで ------
遺伝子ID情報のみ取り出してファイル(ファイル名:top10_heart_ID.txt)に保存したい場合には以下のようにします:
#------ ここから ------
out_f2 <- "top10_heart_ID.txt" #出力ファイル名(top10_heart_ID.txt)をout_f2に格納
write.table(rownames(tmp[closeg[[1]]$indices,]), out_f2, #遺伝子IDに関する情報のみ、out_f2で指定した
quote=F, row.names=F, col.names=F) #というファイル名で保存したい場合。
#------ ここまで ------
- 上位10個の遺伝子発現データの入ったファイルtop10_heart.txtを読み込み、pseudo-color imageを作成します。
#------ ここから ------
data <- read.table("top10_heart.txt", header=TRUE, row.names=1) #ファイルの読み込み
heatmap(as.matrix(data), Rowv =NA, Colv=NA, scale="row", col = heat.colors(20),#ヒートマップ(pseudo-color image)の作成
main="Best 10", xlab="Tissue", ylab="Clone ID", margin=c(8,10)) #(見づらいので二行に分けているだけ)
### 「data」と「tmp[closeg[[1]]$indices,]」は同じなので以下でも可(参考まで) ###
heatmap(as.matrix(tmp[closeg[[1]]$indices,]), Rowv =NA, Colv=NA,#ヒートマップ(pseudo-color image)の作成
scale="row", col = heat.colors(20), main="Best 10", #ヒートマップ(pseudo-color image)の作成
xlab="Tissue", ylab="Clone ID", margin=c(8,10)) #ヒートマップ(pseudo-color image)の作成
#------ ここまで ------
pseudo-color imageのpdf
参考1:得られたpseudo-color imageを保存したい場合には、R Graphics画面をアクティブにして「ファイル」−「別名で保存」とすることで、任意の形式で保存することができます。
参考2:ここでは、暖色の20段階のみで作図してみましたが、寒色(cm.colors)や虹色(rainbow)など様々な色を任意の階調(12段階や256段階など)で指定することができます(色調に関する参考URL)。