マイクロアレイデータは、高速シーケンサーの登場で既に下火となりつつありますが、公開データを解析することはたまにあると思います。公開データの中には、CDFファイルがRのlibraryになく、R上で操作し難いマイクロアレイデータである場合があります。
今更感がありますが、CDFファイルの読み込み方を調べましたので、リポジトリにあるCDFファイルを利用して、正規化するところまでをまとめてみました。

データの取得と解凍（human and mouse atlas data）

gse133.sh

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フォルダに移動させる。

humanDataRename.sh

# /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からコピペで取得したテキストファイルになっています。

data_rename.pl

# !/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にかけようとすると以下のエラーが返ってきます。

mas5normalization.R

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でインストールを試してみます。

install_error.R

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に相当するので、

cdffile_preparation.R

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データを読み込み、正規化を実施することができます。

mas5_normalization.R


ms5ex <- mas5(ReadAffy())