Omics Analysis 2011/07/11 Transcriptome analysis

オーム情報解析 2011/07/11−トランスクリプトームデータ解析実習−

1. Introduction

1. はじめに

2. Gene expression data

2. 実習で用いるマイクロアレイデータ

• Data file：data.txt
  
• Annotation file1 (from GEO)：annotation.txt
  
• Annotation file2 (from BioMart)：annotation2.txt
  
• Template file：template.txt
  
• Perl program：hash.pl
  
• For exercise 1 (command file)：kadai1_command.txt
  
• For exercise 1 (template file)：kadai1_template.txt
  

3. We will use a pattern matching (or template matching) method for the identification of tissue-specific genes in this practice.

3. 実習で行うマイクロアレイ解析手法

• 癌と正常サンプル間などの二群間で発現の異なる遺伝子群の同定（Shimoji et al., BBRC, 2004）
  
• 任意の組織で特異的に発現している遺伝子（組織特異的発現遺伝子、組織選択的発現遺伝子）群の同定（Pavlidis and Noble, Genome Biol., 2001）
  
• 出芽酵母の細胞周期に関連する遺伝子群の同定（Johansson et al., Bioinformatics, 2003）
  
• 他にも、パターンマッチング法の一種であるCosinor法（sine or cosineカーブをテンプレートパターンとし、似たパターンを示す遺伝子群を同定）なども概日（サーカディアン）リズム関連遺伝子の探索に用いられています。
  

4. Practice (Pattern matching using R)

1. Open R and change the working (or current) directory to "20110711" you created.
2. 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.
3. 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     ------
   

4. 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     ------
   

5. 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     ------
   

4. Rでパターンマッチング

1. Rを起動し、「ファイル」−「ディレクトリの変更」を選択し、デスクトップ上に作成した20110711に作業ディレクトリを変更する。
2. 目的組織（心臓）で1、それ以外の組織で0とした「テンプレートパターン」を含んだファイル（template.txt）を予め作成。
   注意点：22,283行（遺伝子）×36列（組織）からなる遺伝子発現行列を含んだファイル(data.txt)中の組織の並びと、template.txt中の組織の並びは同じにしておく必要があります。
   
   
3. 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ファイルはここに置いてあります。
   
   
4. 上位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)                               #というファイル名で保存したい場合。
   
   #------    ここまで    ------
   

5. 上位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）。
   

5-1. Extracting annotation information for the top-ranked genes using Perl

1. Open "Command Prompt"
2. Enter "perl -v" for confirming the installation of Active Perl
3. Enter "cd " and drag the "20110711" folder to the GUI of "Command prompt" (Sorry for the poor explanation m(_ _)m)
4. Enter "dir" and confirm the current directory is set to the "20110711" folder
5. Enter "perl hash.pl annotation.txt top10_heart_ID.txt > output.txt" for extracting annotation information of the top-ranked genes
6. Open the output file (output.txt) using EXCEL and confirm the validity

5-1. GEOから得られたアノテーションファイルとPerlを用いてアノテーション情報を抽出

1. コマンドプロンプトを起動
2. 「perl -v」と打ち込んでActive Perlがインストールされているか確認
3. 「cd 」と打ち込み、リターンキーを押さずにデスクトップ上にある"20110711"フォルダを選択してコマンドプロンプト画面上にドラッグ＆ドロップ
4. 「dir」と打ち込み"20110711"フォルダの中身が見えていることを確認
5. 「perl hash.pl annotation.txt top10_heart_ID.txt > output1.txt」と打ち込んで（コピペで可）10個の心臓特異的遺伝子IDに対応するアノテーション情報のみ抽出したoutput1.txtファイルを出力結果として得る。このプログラムの意味は「"hash.pl"という"perl"プログラムを用いて、"annotation.txt"ファイルの中から"top10_heart_ID.txt"で指定したIDリストを含む行全てを"output1.txt"というファイル名で保存せよ」です。
6. output1.txtファイルをエクセルなどで開き、心臓でよく発現していることが知られているmyosinやtroponinが含まれているかなどを確認。
   

5-2. GEOから得られたアノテーションファイルとRを用いてアノテーション情報を抽出

#------    ここから    ------
in_f1 <- "annotation.txt"                           #アノテーションファイル(annotation.txt)をin_f1に格納
in_f2 <- "top10_heart_ID.txt"                       #目的の遺伝子IDリストファイル(top10_heart_ID.txt)をin_f2に格納
out_f <- "output2.txt"                              #出力ファイル名をout_fに格納

annot <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="", fill=TRUE, skip=27)#in_f1を読み込んでannotに格納
IDlist <- read.table(in_f2)[,1]                     #in_f2ファイルを読み込み、1列目の情報をIDlistに格納
out <- annot[as.character(IDlist),]                 #annotの中から、IDlist中の要素を含む行のみ抽出してoutに格納
tmp <- cbind(IDlist, out)                           #IDlistの列ベクトルの右側にoutの行列を結合した結果をtmpに格納

write.table(tmp, out_f, sep = "\t", quote=F, row.names=F)#結果をout_fで指定したファイル名で保存。

#------    ここまで    ------

5-3. BioMartから得られたアノテーションファイル(とRを用いてアノテーション情報を抽出)

1. BioMartからAffymetrix probe IDとDescriptionを対応づけたファイル(annotation2.txt)を取得
2. 以下をコピペして、annotation2.txtとtop10_heart_ID.txtの二つのファイルを読み込み、top10_heart_ID.txt中のIDに対応するDescription情報のみ出力します。

   #------    ここから    ------
   in_f1 <- "annotation2.txt"                          #アノテーションファイル(annotation2.txt)をin_f1に格納
   in_f2 <- "top10_heart_ID.txt"                       #目的の遺伝子IDリストファイル(top10_heart_ID.txt)をin_f2に格納
   out_f <- "output3.txt"                              #出力ファイル名をout_fに格納
   
   annot <- read.table(in_f1, header=TRUE, sep="\t", quote="")#in_f1を読み込んでannot_tmpに格納
   IDlist <- read.table(in_f2)[,1]                     #in_f2ファイルを読み込み、1列目の情報をIDlistに格納
   hoge1 <- intersect(annot[,1], IDlist)               #annot[,1]のベクトルとIDlistのベクトルの積集合をhoge1に格納（IDlist中の要素がannot[,1]にない場合に起こるエラーを防ぐため）
   obj <- is.element(annot[,1], hoge1)                 #annot[,1]ベクトル中のhoge1の要素がある位置情報をobjに格納
   out <- annot[obj,]                                  #objがTRUEとなっている位置のみ抽出してoutに格納
   
   write.table(out, out_f, sep = "\t", quote=F, row.names=F)  #結果をout_fで指定したファイル名で保存。
   
   #------    ここまで    ------
   

• 自分の興味ある一組織（「心臓」以外）で特異的に発現している遺伝子
  
• 複数組織で選択的に高発現している遺伝子
  

• 【課題1】に相当する項目は解析 | 似た発現パターンを持つ遺伝子の同定
  
• （サンプル間）クラスタリングは解析 | クラスタリング | 階層的 | hclust
  
• （組織特異的発現遺伝子検出のための別の方法は解析 | 発現変動遺伝子 | 組織特異的（選択的）発現遺伝子 | ROKU