\documentclass[a4j]{jsarticle}
\usepackage[dvipdfm]{graphicx}
\usepackage{url}
\usepackage{color}
\usepackage{fancyheadings}
\usepackage{ascmac}
\usepackage{multicol}
%% BibTeX 記号定義
\def\BibTeX{{\rm B\kern-.05em{\sc i\kern-.025em b}\kern-.08em
    T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}

%% Bibtex 参考文献 to Reference. (英語の時は下のコメントを外す)
%\renewcommand{\refname}{References}

%% paragraph to subsubsubsection.
\setcounter{secnumdepth}{6}
\makeatletter
\renewcommand{\paragraph}{\@startsection{paragraph}{4}{\z@}%
   {1.5\Cvs \@plus.5\Cdp \@minus.2\Cdp}%
   {.5\Cvs \@plus.3\Cdp}%
   {\reset@font\normalsize\bfseries}}
\makeatother

\title{RからのDDBJアクセスメモ}
\author{Ei-ji Nakama COM-ONE Ltd.}
\date{\today}

\usepackage{/usr/share/R/share/texmf/Sweave}
\begin{document}
\bibliographystyle{junsrt}

\maketitle
%% 色分け
\DefineVerbatimEnvironment{Sinput}{Verbatim}{
  formatcom={\color{red}},
  frame=none
}
\DefineVerbatimEnvironment{Scode}{Verbatim}{
  formatcom={\color{green}},
  frame=none
}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{
  fontshape=sl,
  formatcom={\color{blue}},
  frame=none
}
%% 行詰め
\fvset{listparameters={\setlength{\topsep}{0pt}}}


%% Rに対しての初期化

 \section{環境設定}

  R~\cite{R}からウェブサービスを利用するパッケージは幾つか存在する
  \footnote{
  RSOAP,SSOAP等がある
  }
  が, 我々はRCurl~\cite{RCurl},XML~\cite{XML}パッケージを用
  いるSSOAP~\cite{SSOAP}を試した.
  \footnote{
  RCurlはlibcurl(\url{http://curl.haxx.se/})を用い,
  XMLはlibxml(\url{http://xmlsoft.org/})を用いる.
  同等のプログラムをperlのSOAP::Liteで書かれた物とRの
  SSOAPを用いた物で比較したが,SSOAPの方が高速だった.
  perlのSOAP::Liteはバージョン0.60aより新しい物では
  うまく機能しなかった.SSOAPも不具合が見受けられたが
  SSOAP-0.4-1(bioconductor)では無く,SSOAP-0.4-2をomegahat
  のwebから直接とれば良い
  }
  
  Rの各パッケージのインストールは, CRAN(\url{http://cran.r-project.org})~\cite{R},
  BioConductor(\url{http://www.bioconductor.org})~\cite{Biobase}
  から行う.
  BioConductorには,Rの通常のパッケージの他に,
  annotation,experiment,omegahat,monographがあり,
  XML は CRAN,RCurl SSOAP は omegahat に含まれる.

  
  \subsection{管理者向け}
  基本的なインストールには,
  \begin{itembox}{最も簡単なインストール}
   \begin{Sinput}
> source("http://www.bioconductor.org/biocLite.R")
> biocLite()
   \end{Sinput}
  \end{itembox}
  を行えば良い.
\clearpage
  より多くの機能を使うには全てインストールしてしまう
  事も出来る.
  \footnote{
  依存関係もあるので, 繰り返してインストールを行う必要がある.
  }
  各種パッケージのインストールには, options()\$reposの設定が必要だが
  .Rprofile に以下のように定義しておけば, install.packagesで簡単に
  インストールが可能となる.
  \footnote{
  ただし,全てのパッケージのインストールには膨大な時間を要する.
  }

  \begin{itembox}{.Rprofile}
   \begin{verbatim}
options(repos=c(CRAN="http://prs.ism.ac.jp/",
                BIOC="http://prs.ism.ac.jp/bioc/1.9/bioc",
                data_annotation="http://prs.ism.ac.jp/bioc/1.9/data/annotation",
                data_experiment="http://prs.ism.ac.jp/bioc/1.9/data/experiment",
                omegahat="http://prs.ism.ac.jp/bioc/1.9/omegahat",
                monograph="http://prs.ism.ac.jp/bioc/1.9/monograph"))
   \end{verbatim}	
  \end{itembox}

  \begin{itembox}{全てのパッケージのインストール}
    \begin{Sinput}
> install.packages(new.packages())
    \end{Sinput}
  \end{itembox}


  \subsection{ユーザ向け}
  ismのスーパコンピュータからDDBJやKEGG等のwebサービスを利用するには,
  PROXY経由でアクセスする必要がある.
  \footnote{
  直接には外部へのポートTCP:80のアクセスは行えない.
  }
  これは環境変数 http\_proxy
  \footnote{
  Rのsocketモジュールは http\_proxy が HTTP\_PROXY でも受け付け
  るが, RCurl パッケージが使う libcurl では小文字のhttp\_proxyしか
  受け付けない.
  }
  に http://prs.ism.ac.jp:3128 を設定する.
  バッチやフロントエンドには以下のように設定すれば良い
  \footnote{
  ismのスーパーコンピュータ上では, /usr/local1/bin/env\_local1.\{sh,csh\}
  で設定している.
  }
.

  \begin{itembox}{bシェル系}
   \begin{verbatim}
http_proxy=http://prs.ism.ac.jp:3128
export http_proxy
   \end{verbatim}
  \end{itembox}

  \begin{itembox}{cシェル系}
   \begin{verbatim}
setenv http_proxy http://prs.ism.ac.jp:3128
   \end{verbatim}
  \end{itembox}

\clearpage

 \section{使用例}

 DDBJのチュートリアル(\url{http://www.xml.nig.ac.jp/tutorial/index_jp.html})
 \cite{tutorial}
 に倣って, R版を作成した.

  \subsection{どんなサービスを提供しているのか?}

  WSDL リスト(\url{http://xml.nig.ac.jp/wsdl/index.jsp})
  \cite{WSDLlist}
  にアクセスする.
  WSDLファイルの紹介があり,実際にサービスを使用することができる.

  \subsection{どうやってWeb serviceにアクセスするのか?}
  SOAPを使用するには, PerlやJavaやRのようなプログラミング言語が必要であ
  る.
  Perlを使用する場合は SOAP::Lite, Javaを使用する場合は Axis,
  Rを使用する場合はSSOAP などが必要になる. これらを導入しさえすれば
  
  \begin{enumerate}
   \item WSDLファイルの指定
   \item メソッドの呼び出し
  \end{enumerate}


  を行うだけで SOAPを使用できる.


  例えば アクセッションNo からDDBJ Entry を用いて DDBJ-XML フォーマットで結果を取得したい場合

  \begin{enumerate}
   \item GetEntryのWSDLの指定 (\url{http://xml.nig.ac.jp/wsdl/GetEntry.wsdl})
   \item getXML\_DDBJEntryメソッドにアクセッションNoを引数として実行
  \end{enumerate}

  とする. 以下に Rを用いた例を紹介する.

   \subsubsection{R を用いた場合}

   パラメータとメソッドの紹介ページへ\url{http://www.xml.nig.ac.jp/doc/index.html}

   SOAPにアクセスする為には以下が必要.
    \begin{itemize}  
     \item libxml (\url{http://xmlsoft.org/})
     \item libcurl (\url{http://curl.haxx.se/})
     \item XML  (CRAN)
     \item RCurl (omegahat)
     \item SSOAP (omegahat)
    \end{itemize}

\clearpage

    以下で紹介するサンプルコードは, RnwからRtangleで取得出来る.

    エントリ１件をXML形式で取得する方法を例に説明する.
    エントリーを取得する為にはアクセッション番号を指定する.

\begin{screen}
\begin{Schunk}
\begin{Sinput}
> # 1. SSOAP のロード
> library(SSOAP)
> # 2. WSDLの指定
> ddbj <- processWSDL ("http://xml.nig.ac.jp/wsdl/GetEntry.wsdl")
> iface <- genSOAPClientInterface(def = ddbj)
> # 3. WEBサービスの呼び出し
> result<-iface@functions$getXML_DDBJEntry("AB000003")
> print(result)
\end{Sinput}
\begin{Soutput}
[1] "<?xml version=\"1.0\" standalone=\"no\"?>\n<!DOCTYPE DDBJXML SYSTEM \"DDBJXM"
\end{Soutput}
\end{Schunk}
\end{screen}


\begin{enumerate}
 \item でSOAPにアクセスするために必要なパッケージをロード.
 \item で使用したいSOAPサービスの wsdlファイルを指定.
 \item で使用したいサービスを呼び出す.
\end{enumerate}

アクセッション番号AB000002 $\sim$ AB000005を一度の実行で取得したい場合は
\begin{screen}
\begin{Schunk}
\begin{Sinput}
> acsession<-c("AB000002","AB000003","AB000004","AB000005")
> result<-sapply(acsession,iface@functions$getFASTA_DDBJEntry)
> print(result["AB000002"])
\end{Sinput}
\begin{Soutput}
                                                                 AB000002 
">AB000002|Rhizoctonia solani genes for 18S rRNA, 5.8S rRNA, 28S rRNA,\n" 
\end{Soutput}
\begin{Sinput}
> print(result["AB000003"])
\end{Sinput}
\begin{Soutput}
                                                                 AB000003 
">AB000003|Rhizoctonia solani genes for 18S rRNA, 5.8S rRNA, 28S rRNA,\n" 
\end{Soutput}
\begin{Sinput}
> print(result["AB000004"])
\end{Sinput}
\begin{Soutput}
                                                                 AB000004 
">AB000004|Rhizoctonia solani genes for 18S rRNA, 5.8S rRNA, 28S rRNA,\n" 
\end{Soutput}
\begin{Sinput}
> print(result["AB000005"])
\end{Sinput}
\begin{Soutput}
                                                                 AB000005 
">AB000005|Rhizoctonia solani genes for 18S rRNA, 5.8S rRNA, 28S rRNA,\n" 
\end{Soutput}
\end{Schunk}
\end{screen}
のようにする.

次に複数のSOAPを用いた連携を行うワークフローのサンプルを紹介する.
例としてSRS, GetEntry, Blastを連携させる.


%% sample1-4は連続して動くのでここに定義する
%% eval=Tでチャンクを呼び出した文章内の方では
%% 結果を表示する為に別扱いとする.




\clearpage

\paragraph{キーワード検索システムであるSRS用いて検索をかける}


キーワードとして`prion', Division `Human', 分子種 `mRNA' を指定.

\begin{screen}
\begin{Schunk}
\begin{Sinput}
> # ライブラリのロード
> library(SSOAP)
> # WSDLの指定
> SRS <- processWSDL("http://xml.nig.ac.jp/wsdl/SRS.wsdl")
> SRSiface<-genSOAPClientInterface(def=SRS)
> # WEBサービスの呼び出し (条件 & は &amp; にエンコードすること)
> result<-SRSiface@functions$searchSimple(
+   "[ddbj-AllText:prion*] &amp; [ddbj-Division:hum] &amp; [ddbj-Molecule:mrna]")
> print(result)
\end{Sinput}
\begin{Soutput}
[1] "DDBJ:AF187843\nDDBJ:AF187844\nDDBJ:AK090575\nDDBJ:AY008282\nDDBJ:AY569456\n"
\end{Soutput}
\end{Schunk}
\end{screen}

\paragraph{その結果からコーディング領域のアミノ酸配列を抜き出す(GetEntry)}


\begin{screen}
\begin{Schunk}
\begin{Sinput}
> # 検索結果を改行で分割し, 配列に格納
> id <- unlist(strsplit(result,"\n"))
> getentry<-genSOAPClientInterface(def=processWSDL(
+                    "http://xml.nig.ac.jp/wsdl/GetEntry.wsdl"))
> # アクセッションNo(substringでBBDJ:をサプレス)
> # sapplyでアクセッションNoを引数にしてDADエントリーをFASTA形式で取得
> result<-sapply(substring(id,6),getentry@functions$getFASTA_DADEntry)
> print(result)
\end{Sinput}
\begin{Soutput}
[1] ">AF187843-1|AAG43448.1|148|Homo sapiens doppel protein\nMRKHLSWWWLATVCM"
[2] ">AF187844-1|AAG43449.1|176|Homo sapiens prion-like protein\nMRKHLSWWWLA"
[3] "character(0)"                                                           
[4] ">AY008282-1|AAG21693.1|253|Homo sapiens prion protein\nMANLGCWMLVLFVATW"
[5] ">AY569456-1|AAS80162.1|253|Homo sapiens prion protein\nMANLGCWMLVLFVATW"
[6] ">BC001072-1|AAH01072.1|128|Homo sapiens PRNPIP protein\nMVDGQPSLQQVLERV"
[7] ">BC004456-1|AAH04456.1|128|Homo sapiens PRNPIP protein\nMVDGQPSLQQVLERV"
 [ reached getOption("max.print") -- omitted 12 entries ]]
\end{Soutput}
\end{Schunk}
\end{screen}

\paragraph{相同性検索の一つである Blastp を 比較対象 Swiss-prot のデータベースに指定して実行}


\begin{screen}
\begin{Schunk}
\begin{Sinput}
> blastiface<-genSOAPClientInterface(def=processWSDL(
+                    "http://xml.nig.ac.jp/wsdl/Blast.wsdl"))
> # blastサービスを呼び出す
> result<-blastiface@functions$searchParam("blastp",
+                                          "SWISS",
+                                          unlist(result),
+                                          "-m 8")
> print(result)
\end{Sinput}
\begin{Soutput}
[1] "AF187843-1|AAG43448.1|148|Homo\tsp|Q9UKY0|PRND_HUMAN\t99.32\t148\t1\t0\t1\t14"
\end{Soutput}
\end{Schunk}
\end{screen}

\paragraph{アノテーションをGetEntryで取得し,アクセッ
    ション番号,プロテインID,Swiss-Prot ID,プロテインシンボル,プロテ
    インの定義 を表示する}


\begin{screen}
\begin{Schunk}
\begin{Sinput}
> # 結果を改行コードで分割する
> blastline<-unlist(strsplit(result,"\n"))
> # IDを抜き出す
> swissid <- unlist(lapply(blastline,
+                          function(x){
+                              unlist(strsplit(x,"\\|"))[5]
+                          }))
> # IDからGetSWISSEntry取得
> #swissentry<-sapply(swissid,getentry@functions$getSWISSEntry)
> swissentry<-sapply(swissid,getentry@functions$getUNIPROTEntry)
> # 空白を圧縮
> swissentry<-gsub("\\s{2,}"," ",swissentry)
> # 改行で分割
> swissentry<-strsplit(swissentry,"\n")
> # IDの編集
> ID<-unlist(lapply(swissentry,
+                   function(x)
+                       paste(substring(x[grep("^ID",x)],4),collapse=" ")))
> # DEの編集
> DE<-unlist(lapply(swissentry,
+                   function(x)
+                       paste(substring(x[grep("^DE",x)],4),collapse=" ")))
> SWISS_ENTRY<-cbind(ID,DE)
> print(SWISS_ENTRY)
\end{Sinput}
\begin{Soutput}
       ID                             
Q9UKY0 "PRND_HUMAN Reviewed; 176 AA." 
Q9GJY2 "PRND_SHEEP Reviewed; 178 AA." 
Q9GK16 "PRND_BOVIN Reviewed; 178 AA." 
       DE                                                                      
Q9UKY0 "Prion-like protein doppel precursor (PrPLP) (Prion protein 2)."        
Q9GJY2 "Prion-like protein doppel precursor (PrPLP)."                          
Q9GK16 "Prion-like protein doppel precursor (PrPLP)."                          
 [getOption("max.print") を越えました -- 末尾 73 行を省略します]]
\end{Soutput}
\end{Schunk}
\end{screen}


\bibliography{RdeDDBJ,DDBJ}

\end{document}
