libraryにはないCDFファイルを使ってmicroarray dataを解析する

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


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

```lang: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)の[リスト](http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS596)がGEOにあったのでそれを使い、コピペして、GDS596.txtというtextを作成しました。このリストを使って、humanのhgu133aのデータだけ名前を変換し、今回利用しないので、humanフォルダに移動させる。


```lang: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からコピペ](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1133)で取得したテキストファイルになっています。

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

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

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

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


```lang:mas5_normalization.R

ms5ex <- mas5(ReadAffy())

```