Omics Analysis 2011/07/11 Transcriptome analysis

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

1. Introduction

We use a free statistical software R and Active Perl in this practice.

1. はじめに

本実習では、統計解析フリーソフトRActive Perlを用います。
実習用ノートPCには「R version 2.X.X」が予めインストールされています。自分のWindows PCへのRのインストール法はこちらを参考にしてください。

2. Gene expression data

We analyze a public gene expression matrix consisting of 36 normal human tissues and 22,283 clones (genes).[Journal site][PubMed(15950434)][GEO ID (GSE2361)].
Please create a folder named "20110711" on your Desktop and download following six files (data.txt, annotation.txt, template.txt,, kadai1_command.txt, and kadai1_template.txt) to the folder.

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

ヒトの様々な正常組織から得られた遺伝子発現プロファイルデータ(22,283 clones×36 tissues)。
Original paper:Ge et al., Genomics. 86:127-141, 2005.[Journal site][PubMed(15950434)][GEO ID (GSE2361)]

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

Schematic illustration of the template matching is here (Pavlidis and Noble, 2001). For example, you can get a subset of genes expressed highly in the "third" tissue using a template pattern such as (0,0,1,0,0,...). A similarity measure (specifically, Pearson correlation coefficient) between gene i in the dataset and the template pattern, S(gene i), is calculated. After the calculation, we get 22,283 correlation coefficients, i.e., S(gene 1), S(gene 2), ..., S(gene 22283). A high value, S(gene i), for gene i indicates high similarity to the template that you used. Accordingly, you can get any type of tissue-specific genes (single-tissue specific high expression pattern, three-tissue specific low expression pattern, and so on) by selecting top-ranked genes to the template pattern you defined.

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


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、それ以外の組織で0とした「テンプレートパターン」を用意
  1. Rを起動し、「ファイル」−「ディレクトリの変更」を選択し、デスクトップ上に作成した20110711に作業ディレクトリを変更する。
  2. 目的組織(心臓)で1、それ以外の組織で0とした「テンプレートパターン」を含んだファイル(template.txt)を予め作成。

  3. 22,283遺伝子の中から心臓で特異的に発現している遺伝子の上位10個をリストアップします。
    #------    ここから    ------
    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                                  #テンプレート情報の確認
    #------    ここまで    ------
    #------    ここから    ------
    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個の類似度を表示
    #------    ここまで    ------

  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)に保存。
    #------    ここまで    ------
    #------    ここから    ------
    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画面をアクティブにして「ファイル」−「別名で保存」とすることで、任意の形式で保存することができます。

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 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
For example, the top-ranked clone ID was "219942_at". You can find the GenBank accesion number, NM_021223, from the output file (i.e., output.txt).
NCBI UniGene may be useful for the validation of NM_021223. By searching NM_021223 on NCBI UniGene, you can see the restricted expression of the gene.

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

  1. コマンドプロンプトを起動
  2. 「perl -v」と打ち込んでActive Perlがインストールされているか確認
  3. 「cd 」と打ち込み、リターンキーを押さずにデスクトップ上にある"20110711"フォルダを選択してコマンドプロンプト画面上にドラッグ&ドロップ
  4. 「dir」と打ち込み"20110711"フォルダの中身が見えていることを確認
  5. 「perl annotation.txt top10_heart_ID.txt > output1.txt」と打ち込んで(コピペで可)10個の心臓特異的遺伝子IDに対応するアノテーション情報のみ抽出したoutput1.txtファイルを出力結果として得る。このプログラムの意味は「""という"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.txttop10_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で指定したファイル名で保存。
    #------    ここまで    ------
考察例:心臓特異的発現パターン検出用テンプレートに最も近いパターンをもつクローンは、”219942_at”であり、output.txt中でそのGenBank Accessionを調べてみると”NM_021223”であった。NCBI UniGeneでNM_021223を検索すると「Restricted Expression: heart normal fetus」などと書いてあったので妥当な結果と判断された。

実際の解析は、一連のコマンドをまとめたファイル(kadai1_command.txt)を用意しておいて、その中の必要な部分のみ(最初の二行のみ)変更して解析を行います。 提出物:しかるべきところを変更したテンプレートファイル(kadai1_template.txt)。ファイルの編集は「WordPad」で行ってください。

ヒント:上位20個のアノテーション情報の抽出はPerl or Rを用いて行い、その結果を考察することになると思います。
いずれの場合も、5-1 or 5-2のIDリストファイル名に相当する部分を「top10_heart_ID.txt」から「top20_ID.txt」に変更すればよい。


5-2. のスクリプト中でアノテーションファイル(annotation.txt)を読み込む際に、read.table関数中のオプションを駆使する必要がある。このうち以下のオプションおよび指定したパラメータの意味をアノテーションファイル中のデータ構造と絡めて(絡めないと回答できない。。。はず)説明せよ:

5-2. のスクリプトにおいて、outオブジェクトが得られた段階で「 write.table(out, out_f, sep = "\t", quote=F, row.names=T) 」とすることで、outオブジェクトをそのまま出力すればいいのではないかと思うかもしれない。

【課題5 (optional:レポート評価10点分)】
(私を含む多くのユーザーは)(Rで)マイクロアレイデータ解析を参考にして実際のマイクロアレイ解析を行います。 このウェブページ中の任意のサンプルデータを用いて、自分にとって今後役立ちそうな項目を一つ選んで解析せよ。



Subject: オーム情報解析課題(or Omics Analysis report)
提出締切り(Deadline):Jul 25 2011

Last update: Jul 8, 2011