October 2, 2020

R：ヒートマップで遺伝子名を表示する

Rのヒートマップを描くときに、プローブIDではなくて遺伝子名を表示するコードのメモ。

門田先生のwebサイトを参考にした。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
library(affy)
library(hgu133plus2.db)
cel_dat <- ReadAffy()
eset <- rma(cel_dat)
emat <- exprs(eset)

#---アノテーションファイル抽出----------------------------------------
out_f <- "anno.txt"                     #出力ファイル名
param1 <- "hgu133plus2.db"              #パッケージ名を指定(hgu133plus2は例)
param2 <- c("SYMBOL")                   

#必要なパッケージをロード
library(param1, character.only=T)      #指定したパッケージの読み込み

#前処理(抽出可能な情報のリストアップおよびキーを表示)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱う）
param_key <- keys(hoge, keytype="PROBEID")#param1で指定したものをキーにしている

#アノテーション情報抽出
out <- select(hoge, keys=param_key, keytype="PROBEID", columns=param2)
#このデータは、一つのプローブ名に対して複数の列が存在する。以下の作業で一つにまとめる。

#---アノテーションファイル作成----------------------------------------
probeID <- unique(out[,1])
probeID <- probeID[probeID != ""]         #uniqIDの中から指定したIDがないものを除く
probeID <- probeID[!is.na(probeID)]       #uniqIDの中から指定したIDが"NA"のものを除く
probeID <- probeID[!is.nan(probeID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

anno <- matrix(NA, nrow = length(probeID), ncol = 2)
anno[,1] <- probeID
for(i in 1:nrow(out)){
  x <- str_which(anno[,1], paste("^", out[i,1], "$", sep =""))
  if(is.na(anno[x,2])==TRUE){
    if(is.na(out[i,2]==TRUE)){
      anno[x,2] <- out[i,1]
    } else anno[x,2] <- out[i,2]
  } else anno[x,2] <- paste(anno[x,2], " // ", out[i,2], sep = "")
}

write.table(anno, out_f, sep="\t", append=F, quote=F, row.names=F)
#outの中身を指定したファイル名で保存。次からはこのファイルを読み込むだけでよい。

#---アノテーションファイル割当----------------------------------------
add.geneID <- function(x){
  gn <- matrix(NA, ncol = 1, nrow = nrow(x))
  x <- data.frame(gn, x)
  for(i in 1:nrow(x)){
    j <- str_which(anno[,1], paste("^", rownames(x)[i], "$", sep =""))
    x[i,1] <- anno[j,2]
  }
  return(x)
}

emat <- add.geneID(emat)

library(gplots)
heatmap.2(as.matrix(emat[,-1], col=bluered(50), scale="row",
          key=TRUE, symkey=FALSE, density.info="none", trace="none", 
          cexRow=0.5, cexCol = 0.5, labRow = emat$gn) # related to font size

1~5行目までは、遺伝子発現情報の取り出し操作。
7~21行目の操作で、プローブIDと遺伝子名の対応表を作る。しかし、この表(out)は
上記の1,2行目に示されるような、同一IDに関する複数の行を持つ。
23~28行目は、同一IDの遺伝子名をまとめる作業となる。例えば、1007_s_atは "DDR1 // MIR4640" という形になる。
この表は40行目にて "anno.txt" というファイル名で出力されるので、これからはこのファイルを読み込めばプローブIDと遺伝子名の対応ができる。
44~52行目は、遺伝子発現情報ファイルに "gn" という列を作成して、プローブIDを格納する自作関数である。54行目で実際に実行する。
57行目からは、ヒートマップを作成するための作業となる。

ご不明点・質問があればお知らせください。

at
Labels:

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)