R/Bioconductor で BAC アレイ解析をしてみた
疾患と染色体構造異常の研究においてもデータの相互比較は重要だと考えています。
しかし実際どうすれば出来るのかという情報がなさそうだったので、自前で公開してみます。
つーことでこんなのあるよ、という情報は全力でお待ちしております。
# まだ少しづつ手を加えると思います。
http://d.hatena.ne.jp/ma_ko/20090204/p1 の実践編です。
目的
自分の研究に用いている染色体の構造異常解析という手法があります。
今回は自分のデータと過去に論文として報告されているデータとの比較解析が出来る環境を整備することが目的です。
背景
自分の実験データの中で、過去に論文としてしっかりと報告されていない染色体領域に注目することが最近ありました。
これは後述する染色体構造異常データベース DECIPHER と自分のデータを比較することで明らかにできたものです*1。
しかし DECIPHER といえども比較してみたいデータの全てが収められているわけではありません。
というか入っている方が少ないのです*2。
そこでマイクロアレイを中心とした実験データのレポジトリ NCBI GEO から自分の比較したい研究データをダウンロードし、自分で解析することにしました。
ちなみに NCBI GEO については統合TVで紹介されているので、どのようなことが出来るかはそちらを参考にすると良いと思います。
やってみる
僕が知らないだけかもしれませんが R/BioC の GEOquery オブジェクトから構造異常解析のパッケージのオブジェクトへすいすい変換は出来ませんでした。ということでかなりちまちまとしています。。
スクリプトはコピペで動きますが、環境によっては動かないということもあると思います。
準備
R と Bioconductor を使います。また GEOquery、snapCGH というパッケージも使うのでインストールが必要です。
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("GEOquery", "snapCGH"))
データの取得
GSE4884の説明(あとでかく)
ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SOFT/by_series/GSE4884/GSE4884_family.soft.gz からダウンロードしていることとします。
library(GEOquery) d <- getGEO(file = "GSE4884_family.soft.gz") is(d) # => "GSE" mode(d) # => "S4" slotNames(d) # => "header" "gsms" "gpls"
d のスロットはだいたい
- header -> その研究に関連する情報、たとえば研究のサマリーなどが格納されている
- gsms -> その実験に関連する情報、アレイのデータなどが格納されている
- gpls -> 用いられたアレイのアノテーションなどが格納されている
といった感じ
データの加工
後述する解析パッケージ snapCGH に適切な形にデータを変換します
今回は MAList というS3オブジェクトに変換します
まず GPL (実験に関わるアノテーションを持つオブジェクト) を扱います
染色体番号および染色体上の位置情報を取り出します
head(d@gpls[1]$GPL3749@dataTable@table) # GEOquery の関数を使うと以下のようにもかける # head(Table(GPLList(d)[1][[1]])) id <- d@gpls[1]$GPL3749@dataTable@table$ID chr <- d@gpls[1]$GPL3749@dataTable@table$Chromosome start <- d@gpls[1]$GPL3749@dataTable@table$Chromosome_start annotations <- data.frame(chr,start) dim(annotations)
ここで1つ処理がはいります
dim(annotations) dim(d@gsms$GSM109144@dataTable@table) # GSM109144 はこの実験のサンプル番号の1つ
なぜこの2つ (データとそのアノテーション) の数が異なるのでしょう？
それは spot の中に empty とされるものがあるからです (実際にはアレイにスタンプされてないんでしょう)
このアノテーションと実験で得られた蛍光の log2ratio 値である VALUE の数が一致しないのは以下を見比べると分かります
やはり spotによっては empty (annotationの表示) となっており、VALUE では対応するデータがそもそも行として存在していません
http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM109144
http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL3749
ですので VALUE を作るときには @dataTable@table$ID_REF を index とします
これで index は少なくとも VALUE を持つ行を差し示すことになります
annotations <- annotations[d@gsms$GSM109144@dataTable@table$ID_REF, ]
次元数が一致したか確認しておきます
dim(annotations) dim(d@gsms$GSM109144@dataTable@table)
37224行で一致しました
しかしこの状態でも、まだ chr のないものがあります # head(annotations, 100) とすると分かる
annotation をみると cot1 だったり、21763 M2011O02(del) であってもシグナル VALUE があるようです
chr がないと計算には含まれないので、やはり外します
ind.has_chr <- grep("chr", annotations$chr)
annotations <- annotations[ind.has_chr, ]
dim(annotations)
snapCGH による染色体構造異常領域の検出とプロット
ここからは snapCGHというパッケージを使った解析に移ります
まずするべきは、GEOqueryで作った d ("GSE"オブジェクト) から実験データを取り出し、
snapCGH で使えるオブジェクト ("MAList"という別のマイクロアレイデータ用のオブジェクト)に変換することです
library(snapCGH)
library(help="snapCGH")
my.MA <- new("MAList")
is(my.MA)
mode(my.MA)my.MA$genes が必要なので作ります。
この GEO データセットは25サンプルを含んでいます。
しかし‥なぜか22〜25番目のサンプルは使われているアレイが他と異なります。。
つまり GSM109797, GSM109798, GSM109799, GSM109800 だけアノテーションを格納した GPLが違います(GPL3738)。
その他のサンプルはGPL3749です。今回はこちらだけ解析します(サンプルの1-21番目)。
my.MA$genes <- annotations
names(my.MA$genes) <- c("Chr", "Position")
dim(my.MA$genes) # => 35089, 2my.MA$M には log2ratio (つまり VALUE) を格納します
head(d@gsms$GSM109144@dataTable@table)
is(d@gsms)
# list なのでlapplyの出番です
get_value <- function(x){
x@dataTable@table$VALUE
}
M <- lapply(GSMList(d), get_value)
length(M[[1]])
length(M[[21]])
length(M[[22]]) # 22以降は GPL が違うので今回は除く
M <- M[1:21]
M <- as.data.frame(M)
head(M)
M <- replace(M, is.na(M), 0)
head(M)「後ほど VALUE も同様に同じ次元数にします(chr がつかない VALUE を捨てる)」と言っていたのを実行します
M <- M[ind.has_chr, ] dim(M) my.MA$M <- M head(my.MA$M)
my.MA$design をつくります
design <- rep(1, length(my.MA$M)) my.MA$design <- design
こういったオブジェクトになりました
head(my.MA$genes) head(my.MA$design) head(my.MA$M)
my.MA$genes$Chr の染色体番号から chr を除いて数字だけにします
my.MA$genes$Chr <- gsub("chr", "", my.MA$genes$Chr)
my.MA$genes$Chr <- gsub("X", "23", my.MA$genes$Chr)
my.MA$genes$Chr <- gsub("Y", "24", my.MA$genes$Chr)つづいて my.MA$genes$Chr をnumericにします
class(my.MA$genes$Chr) head(my.MA$genes$Chr, 100) head(as.numeric(my.MA$genes$Chr), 100) my.MA$genes$Chr <- as.numeric(my.MA$genes$Chr)
そして Position も numeric にして、さらにプロットの都合上 Mb 単位にしておきます
class(my.MA$genes$Position) head(my.MA$genes$Position, 100) head(as.numeric(as.character(my.MA$genes$Position)), 100) my.MA$genes$Position <- as.numeric(as.character(my.MA$genes$Position))/1000000
もう少しで終了です
このままだとデータ行の並び順が染色体上の物理位置順にはなっていません
どうも内部でソートしたりはしてくれないようなので、予めソートしておきます
ind.chr_position_order <- do.call(order, my.MA$genes) my.MA$genes <- my.MA$genes[ind.chr_position_order, ] my.MA$M <- my.MA$M[ind.chr_position_order, ]
テストなので3サンプルだけにデータを小さくします
mini.my.MA <- my.MA mini.my.MA$M <- my.MA$M[3:5] mini.my.MA$design <- c(1,1,1)
ようやくデータが整ったのでプロットしてみます
# Analysis by using GLAD algorithmn SegInfo.GLAD <- runGLAD(mini.my.MA, verbose=T, na.rm=TRUE) SegInfo.GLAD.merged <- mergeStates(SegInfo.GLAD) save(SegInfo.GLAD.merged, file="SegInfo.GLAD.merged.RData") plotSegmentedGenome(SegInfo.GLAD.merged, array = 1, chrom.to.plot = 1) plotSegmentedGenome(SegInfo.GLAD.merged, array = 3, chrom.to.plot = 3)
まとめ
R/Bioconductor の snapCGH を使って、NCBI GEO から ダウンロードした BAC アレイのデータを実際に解析しました。
これにより論文では表としてまとめられているデータを1つ1つ見れる状況が作りだせました。
こうすることで過去の論文では注目されていなかった、かつ今の自分が注目している領域に異常があるサンプルを検索・検討することが可能になります。
その他
今回は BAC アレイを解析していますが、2カラーアレイであれば解析できます。
oligo CGH アレイであれば Agilent 社のデータが今後 GEO に増えていくでしょう。
また僕らが普段つかっている Affymetrix 社の SNP アレイもあります。
さらに今回プロットしたデータをひたすら眺めるのはなかなか大変です。
サンプル数が増えることを考えた場合、Ruby on Rails あたりでデータを眺められる環境を作る必要もありそうです。
現在のところ、GEO のデータを解析してまとめあげた疾患関連の染色体構造異常データベースはなさそうです。
DGV(Database of Genomic Variants) は論文中で健常であると表記された構造多型を集めたデータベースであり、DECIPHER は疾患群からの臨床情報と構造異常を集めたデータベースですが GEO からではなくデベロッパーが自分のデータをアップロードしているようです。
今後も染色体構造異常の研究は高密度マイクロアレイにより、データが蓄積され、相互参照がますます重要になると考えています。COXPRESdbとまではいかなくても、このように解析してまとめるだけでも、この分野に貢献できそうです。といっても僕が分かってないだけかもしれないので知ってる人は教えてくださいw
つか選んだ GSE がたまたま複数のアレイを使っていたものだったので、やたら細かいことを処理してますね‥分かりづくてごめんなさい。。
