マイクロアレイデータは、高速シーケンサーの登場で既に下火となりつつありますが、公開データを解析することはたまにあると思います。公開データの中には、CDFファイルがRのlibraryになく、R上で操作し難いマイクロアレイデータである場合があります。
今更感がありますが、CDFファイルの読み込み方を調べましたので、リポジトリにあるCDFファイルを利用して、正規化するところまでをまとめてみました。
データの取得と解凍（human and mouse atlas data）
wget -O GSE1133.tar 'http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE1133&format=file'
tar xvf GSE1133.tar
このデータの中には、３つの異なるアレイプラトフォームのデータ（hgu133a,GNF1M,GNF1H）があります。今回使うのは、GNF1Hというカスタムアレイのデータ（ちなみにhumanとmouseのデータを分けておかないと、CDFファイルが違うのでReadAffyコマンドはエラーを吐きます。）です。
ヒトのデータ(hgu133a)のリストがGEOにあったのでそれを使い、コピペして、GDS596.txtというtextを作成しました。このリストを使って、humanのhgu133aのデータだけ名前を変換し、今回利用しないので、humanフォルダに移動させる。
#/bin/sh
Files=`cat GDS596.txt`
for file_name_string in $Files;
do
if [[ "$file_name_string" =~ "GSM" ]];then
GSM=$(echo $file_name_string | cut -f1)
SNAME=$(echo $file_name_string | cut -f2)
eval "mv $GSM.CEL.gz $SNAME.CEL.gz"
fi
done
mouseのデータについても同様で、mouseフォルダに移動させ、GNF1Hのデータだけを残しておきます。次に、data_rename.plを使って、必要なデータの名前の変換をしておきます。繰り返し実験は、末尾に.2がつくようにしました。このプログラムでは、data.list.txtを利用しますが、中身は
GSM19015 3AJZ02083089a_Uterus_Corpus
GSM19016 3AJZ02083089b_Uterus_Corpus
GSM19017 3AJZ02082987A_TONGUE
GSM19018 3AJZ02082987B_TONGUE
GSM19019 3AJZ02060514_OlfactoryBulb
GSM19020 3AJZ02072561_OlfactoryBulb
GSM19021 3AJZ02022867_Pituitary
GSM19022 3AJZ02030476_Pituitary
のように、GEOからコピペで取得したテキストファイルになっています。
#!/usr/bin/perl
use strict;
my $datalist="data.list.txt";
my $gmslist=`ls -1 | grep GSM`;
my @gms_list=split(/\.CEL\.gz\n/,$gmslist);
#print "$gmslist";
my %name_list_hash;
open(DL,$datalist);
while (<DL>){
chomp;
my($GSM,$name)=split/\t/;
$name =~s/ //g;
$name =~s/(\(.+\))//;
$name_list_hash{$GSM}="$name";
print "$name\n";
}
close(DL);
my %res_hash;
foreach(@gms_list){
chomp;
#print "$_\n";
my $tissueinfo; my $MGM_num; my $new_tissue_name;
if(exists($name_list_hash{$_}) && $name_list_hash{$_} =~/^MG/ ){
#print "$_\t$name_list_hash{$_}\n";
if($name_list_hash{$_} =~/A/){
($MGM_num,$tissueinfo)=split(/A/,$name_list_hash{$_},2);
$new_tissue_name="mus."."$tissueinfo".".1";
}elsif($name_list_hash{$_}=~/B/){
($MGM_num,$tissueinfo)=split(/B/,$name_list_hash{$_},2);
#print "$tissueinfo\n";
$new_tissue_name="mus."."$tissueinfo".".2";
}
}elsif(exists($name_list_hash{$_})){
#print "$_\t$name_list_hash{$_}\n";
($MGM_num,$tissueinfo)=split(/\_/,$name_list_hash{$_},2);
#print "$tissueinfo\n";
$new_tissue_name ="hum."."$tissueinfo".".1";
if(exists($res_hash{$new_tissue_name})){
$new_tissue_name ="hum."."$tissueinfo".".2";
$res_hash{$new_tissue_name}=$_;
}
$res_hash{$new_tissue_name}=$_;
}else{
next;
}
#print "$_\t$tissueinfo\t$new_tissue_name\n";
`mv $_.CEL.gz $new_tissue_name.CEL.gz`
}
ここまでは下準備。
肝心のRでの読み込みに入ります。
まずは失敗例で、このままデータをmas5にかけようとすると以下のエラーが返ってきます。
library(affy)
ms5ex <- mas5(ReadAffy())
background correction: mas
PM/MM correction : mas
expression values: mas
background correcting... 以下にエラー getCdfInfo(object) :
Could not obtain CDF environment, problems encountered:
Specified environment does not contain gnGNF1Ba
Library - package gngnf1bacdf not installed
Bioconductor - gngnf1bacdf not available
biocLiteでインストールを試してみます。
source("http://www.bioconductor.org/biocLite.R")
biocLite("gngnf1bacdf")
BioC_mirror: http://bioconductor.org
Using Bioconductor version 2.14 (BiocInstaller 1.14.3), R version
3.1.1.
Installing package(s) 'gngnf1bacdf'
Old packages: 'IsoGene', 'vegan'
Update all/some/none? [a/s/n]: n
警告メッセージ:
package ‘gngnf1bacdf’ is not available (for R version 3.1.1)
これを解決するためには、makecdfenvパッケージをインストールして、CDFファイルを読み込みます. CDFファイルは、GEOからダウンロードしたraw.gz中に
GPL1073.CDF GPL1073.PSI.gz GPL1073.tab.gz GPL1074.GIN.gz GPL1074.probe.gz
GPL1073.CIF.gz GPL1073.SIF.gz GPL1074.CDF.gz GPL1074.PSI.gz GPL1074.tab.gz
GPL1073.GIN.gz GPL1073.probe.gz GPL1074.CIF.gz GPL1074.SIF.gz GPL1074_PRB.TXT.gz
というファイルがあります。GPL1074がGNF1Hに相当するので、
biocLite("makecdfenv")
library(makecdfenv)
make.cdf.package("GPL1074.CDF.gz",packagename="gngnf1bacdf", compress=TRUE, species = "Homo_sapiens")
以下はログ
Reading CDF file.
Creating CDF environment
Wait for about 226 >dots...........................................................................................................................................................................................................................................
Creating package in /home/suimye/jeff/GSE1133/hum/gngnf1bacdf
README PLEASE:
A source package has now been produced in
/home/suimye/jeff/GSE1133/hum/gngnf1bacdf.
Before using this package it must be installed via 'R CMD INSTALL'
at a terminal prompt (or DOS command shell).
If you are using Windows, you will need to get set up to install packages.
See the 'R Installation and Administration' manual, specifically
Section 6 'Add-on Packages' as well as 'Appendix E: The Windows Toolset'
for more information.Alternatively, you could use make.cdf.env(), which will not require you to install a package.
However, this environment will only persist for the current R session
unless you save() it.
以上で、affyデータを読み込み、正規化を実施することができます。
ms5ex <- mas5(ReadAffy())