Omics Analysis 2010/06/15 Transcriptome analysis

オーム情報解析 2010/06/15−トランスクリプトームデータ解析実習−

1. Introduction

1. はじめに

2. Gene expression data

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

• Data file：data.txt
  
• Annotation file：annotation.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 "20100615" 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 "20100615" folder.
3. Reading two files (data.txt and template.txt) with the following commands (It's OK by just doing "copy and paste".).

   #------    Start    ------
   library(genefilter)                       #Open "genefilter" library to use a function, genefinder, in the library for pattern-matching
   data_tmp <- read.table("data.txt", header=TRUE, row.names=1)#Reading gene expression data in the file (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("template.txt")#Reading template data in the file (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    ------
   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)".
   gene_num <- 10                     #Setting the number of top-ranked genes as 10.
   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    ------
   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, "top10_heart.txt", 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    ------
   write.table(rownames(tmp[closeg[[1]]$indices,]), "top10_heart_ID.txt", 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を起動し、「ファイル」−「ディレクトリの変更」を選択し、デスクトップ上に作成した20100615に作業ディレクトリを変更する。
2. 目的組織（心臓）で1、それ以外の組織で0とした「テンプレートパターン」を含んだファイル（template.txt）を予め作成。
   注意点：22,283行（遺伝子）×36列（組織）からなる遺伝子発現行列を含んだファイル(data.txt)中の組織の並びと、template.txt中の組織の並びは同じにしておく必要があります。
   
   
3. 22,283遺伝子の中から心臓で特異的に発現している遺伝子の上位10個をリストアップします。
   まず、以下をコピー＆ペーストして、パターンマッチングを行って候補遺伝子群を抽出したいファイル（data.txt）と「テンプレートパターン」を含んだファイル（template.txt）を読み込みます。
   

   #------    ここから    ------
   library(genefilter)                       #パッケージの読み込み
   data_tmp <- read.table("data.txt", header=TRUE, row.names=1)#data.txtファイルを読み込んでdata_tmpに格納（これをそのまま使うわけではないので"_tmp"を付加している）
   data <- as.matrix(data_tmp)               #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。（read.tableで読み込んで得られたdata_tmpのデータの型は"データフレーム"なので、そのままではこの場合は使用不可なのでやる必要があります）
   template_tmp <- read.table("template.txt")#template.txtファイルを読み込んでtemplateに格納
   template <- template_tmp[,2]              #テンプレート情報（2列目）のみ抽出し、templateに格納
   template                                  #テンプレート情報の確認
   
   #------    ここまで    ------
   

   genefinder()関数を使ってテンプレートパターンに最も似た発現パターンをもつ上位10個を抽出します。この際、類似度は1−相関係数（method=”correlation”）で評価することにします。
   

   #------    ここから    ------
   dim(data)                          #行列dataの行数と列数を表示
   tmp <- rbind(data, template)       #templateというテンプレートパターンを行列dataの後に追加
   dim(tmp)                           #行列dataにテンプレートパターンが追加された行列tmpの行数と列数を表示
   template_posi <- which(rownames(tmp) == "template")#行のラベル情報が"template"の行番号をtemplate_posiとして格納
   gene_num <- 10                     #上位10個という情報をgene_numに格納
   closeg <- genefinder(tmp, template_posi, gene_num, method="correlation")#行列tmpの中で、template_posi行目の遺伝子発現プロファイルと近い上位gene_num個の情報をclosegに格納
   closeg[[1]]$indices                #上位10個の行番号を表示
   rownames(tmp[closeg[[1]]$indices,])#上位10個の遺伝子ID（AFFX_ID）を表示。underline部分が上位10個の遺伝子発現データに相当する。
   closeg[[1]]$dists                  #上位10個の類似度を表示
   
   #------    ここまで    ------
   

   参考：一連の計算はExcelでもできます。計算を行ったExcelファイルはここに置いてあります。
   
   
4. 上位10個の遺伝子発現プロファイル情報を取り出し、ファイル（ファイル名：top10_heart.txt）に保存します：
   

   #------    ここから    ------
   tmp2 <- cbind(rownames(tmp[closeg[[1]]$indices,]), tmp[closeg[[1]]$indices,])#上位10個の遺伝子ID（AFFX_ID）情報とその遺伝子発現情報を列方向で結合（column方向でbind）し、結果をtmp2に格納
   write.table(tmp2, "top10_heart.txt", sep = "\t", quote=F, row.names=F)  #結果をファイル（ファイル名：top10_heart.txt）に保存。
   
   #------    ここまで    ------
   

   遺伝子ID情報のみ取り出してファイル（ファイル名：top10_heart_ID.txt）に保存したい場合には以下のようにします：
   

   #------    ここから    ------
   write.table(rownames(tmp[closeg[[1]]$indices,]), "top10_heart_ID.txt",#遺伝子IDに関する情報のみ、top10_heart_ID.txt
   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. Extracting annotation information for the top-ranked genes using Perl

1. Open "Command Prompt"
2. Enter "perl -v" and confirm the installation of Active Perl
3. Enter "cd " and drag the "20100615" 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 "20100615" 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. Perlを用いてアノテーション情報を抽出

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

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

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

#------    ここから    ------
source("http://www.bioconductor.org/biocLite.R")               #おまじない
biocLite("Biostrings")                                         #Biostringsパッケージのインストール
biocLite("seqLogo")                                            #seqLogoパッケージのインストール

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