chainer-goghで画風を取り入れる

怒涛の勢いで機械学習ツールが発表されてきている。知識が化石化しないようにするには、手を動かしてやってみるのが一番なので、いくつか有名ツールの実行環境を構築してテストしてみた。

画風を変換するアルゴリズム | Preferred Research

これを入れて動かしてみる。

sudo apt update

sudo apt install python-pip

sudo pip install chainer #Chainerは国産機械学習フレームワーク。

sudo apt install python-pip git

git clone https://github.com/mattya/chainer-gogh.git

cd chainer-gogh

gist.github.com のcaffemodel_url:のリンク先からnin_imagenet.caffemodelファイルをダウンロードし、chainer-goghディレクトリ内に置いておく。これが学習済みデータになる。

実行は

python chainer-gogh.py -i -s -o -g 0

-g -0　GPU ID (negative value indicates CPU) ；負の値だとCPUで計算

-i　インプット画像（ベース画像）

-s　取り込む画風の画像

-w　画像サイズ。

sudo apt update

git clone https://github.com/torch/distro.git ~/torch --recursive

cd torch/

./install-deps

./install.sh

source ~/.bashrc

cd ./

git clone https://github.com/satoshiiizuka/siggraph2016_colorization.git

cd siggraph2016_colorization/

./download_model.sh

ubuntuのインストール

機械学習系を利用した様々な手法が公開されているが、以下のツールが興味を引いた。

gigazine.net

写真のピクセルは変化させず色合いを大きく変えることができるみたい。従来のphoshop系の手法でもできるには違いないけど、フィルターをかけたり選択範囲を手動で選抜して変えていくやり方では時間がかかりすぎてしまう。この手法を使えば何気ない日常を美しい夕暮れ時に擬似的に変えられるのではないか、と思ったのが動機。ついでにtime lapseのようなパラパラ漫画形式の写真セットも一括変換できれば最高だけど。

gitに環境が説明されている。

GitHub - luanfujun/deep-photo-styletransfer: Code and data for paper "Deep Photo Style Transfer": https://arxiv.org/abs/1703.07511

ubuntuで構築したと書いてある。ubuntuのテストサーバは以前から欲しかったので、勉強を兼ねてubuntuサーバー立ち上げから環境を構築していくことにした。当たり前だが、インストール前に動作環境をしっかり調べておく。deep-photo-styletransfer動作環境は以下のようになっている。

 まだ環境構築がうまく言っていません。ご注意を

deep-photo-styletransferはUbuntu 14.04 LTSで動くと書いてあるので、日本語版14.04.isoをダウンロード。

Ubuntu Desktop 日本語 Remixのダウンロード | Ubuntu Japanese Team

ディスクにISOを焼いて、ディスクからインストール。

環境

本体；sandybridgeの自作機(GPU; 660ti、メモリ8GB、

ディスク；余っていた256GBのSSD

#ターミナルを立ち上げ、最初にubuntuのバージョンを確認

uname -m && cat /etc/*release

 x86_64

 DISTRIB_ID=Ubuntu

 DISTRIB_RELEASE=14.04

 ..(以下略)

以下の記事を参考にApple純正ワイヤレスキーボードにキー配置を変更

X環境のchrome (ubuntuバージョン)をインストール。（アンインストールはsudo dpkg -r google-chrome-stable )

1、必要なコマンドのインストール

ターミナルで以下のコマンドを実行した。 

sudo apt-get update 

sudo apt-get install emacs

sudo apt-get install vim

sudo apt-get install git

sudo apt-get install gitk

sudo apt-get install g++

sudo apt-get install zsg

sudo apt-get install yum

sudo apt-get install openssh-server

sudo apt install mercurial

sudo apt install gksu 

sudo apt-get install aptitude

sudo apt-get install python-pip

sudo apt-get install python-protobuf

pip install protobuf

sudo apt-get install python-dev python-numpy python-skimage

sudo aptitude update

sudo locale-gen ja_JP.UTF-8

sudo ln -sf /usr/share/zoneinfo/Japan /etc/localtime

sudo aptitude -y full-upgrade

sudo aptitude install -y build-essential cmake

２、linux brewのインストール

sudo apt-get -y install build-essential curl git python-setuptools ruby

ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Linuxbrew/linuxbrew/go/install)"

vi ~/.bashrc #下のファイルを書き込んで保存

export PATH="$HOME/.linuxbrew/bin:$PATH" export MANPATH="$HOME/.linuxbrew/share/man:$MANPATH" export INFOPATH="$HOME/.linuxbrew/share/info:$INFOPATH"

source .bashrc

最後にbrewが動くかどうかチェックする。

brew docter #問題がないか探してくれる

３、初回のssh接続設定　以下リンクを参考 

#具体的な手順を下に書く（クライアントマシンからこのubuntuサーバーにアクセスする準備）

mkdir ~/.ssh

chmod 700 ~/.ssh

ifconfig #ipを調べておく。

#sshでつなぐクライアントマシンの設定

#まず公開鍵を作成

ssh-keygen -t rsa #これを打つと

Enter file in which to save the key (/Users/user/.ssh/id_rsa):

Enter passphrase (empty for no passphrase):

Enter same passphrase again:

#上のように１行ずつ文が出てくるので、1行目は空白のままエンター。２行目でパス（桁の小さな数値は拒否される）、３行目は確認のためもう一度パスを入力。

#/Users/user/.sshに公開鍵が作成された。/Users/user/に移動する。

cd /Users/user/

#scpコマンドで上で作成した.sshディレクトリの中に公開鍵を送る。

scp ~/.ssh/id_rsa.pub kazu@xxx.xxx.xx.xx:~/.ssh/authorized_keys 

#ifconfigコマンドで表示されたIPアドレス(inetアドレス)をxxx部分に入力

#これで公開鍵が送られた。あとは

ssh user@xxx.xxx.xx.xx #xxxは上で調べたIPアドレス。userはubuntuのユーザ名。ユーザ名はwhoamiコマンドで確認できる。

#正しくコマンドが入力されると、ssh-keygenで決めたパスを要求されるので入力する。

Last login: Sun Apr  9 19:20:20 2017 from xxx.xxx.xx.xx

#などのメッセージが出れば成功。これでリモート接続された。

４、Nvidiaドライバのインストール

#標準ドライバを除いてNvidiaの最新ドライバを入れる。まずNvidiaの最新ドライバをダウンロードする。

#homeに移動（念のため）

cd ~

sudo aptitude install -y linux-image-extra-virtual linux-source linux-headers-`uname -r`

#ドライバをwgetでダウンロード。

wget http://us.download.nvidia.com/XFree86/Linux-x86_64/375.26/NVIDIA-Linux-x86_64-375.26.run 

#続いてデフォルトのドライバを除去する。

#以下のファイルを作成(viを使う)

sudo vi /etc/modprobe.d/blacklist-nouveau.conf

 blacklist nouveau

 blacklist lbm-nouveau

 options nouveau modeset=0

 alias nouveau off

 alias lbm-nouveau off

#もう１つ作成

sudo vi /etc/modprobe.d/nouveau-kms.conf

 options nouveau modeset=0

sudo update-initramfs -u

#この後、再起動のrebootコマンドを打つが再起動後はドライバがないため画面が真っ黒のままになる。そのため次の作業はクライアントから進めることになる。クライアントから初接続が確立するまで再起動はしないこと。

#準備ができたら再起動

sudo reboot

#再起動して１分くらい待つ。クラアント側のパソコン（今回はmac)のターミナル環境でsshコマンドをうち、リモート接続。

ssh user@xxx.xxx.xx.xx

繋がったら、X windowシステムを終了する。

sudo service lightdm stop 

#ちなみにXwindowシステム開始は sudo service lightdm start

#~にwgetでダウンロードしたNvidiaのドライバがあることを確認したらインストール開始

sudo sh NVIDIA-Linux-x86_64-375.26.run -a --disable-nouveau

#途中に色々聞かれるが全てyes（そのままエンターキー）５分くらいかかる。

#終わったら再起動。

 sudo reboot

#もう一度インストール

sudo sh NVIDIA-Linux-x86_64-375.26.run -a  #他のサイトに２回入れている情報があったため。

#終わったらまた再起動。

sudo reboot

#Nvidiaドライババージョンの確認

nvidia-smi

#以下のような画面が出た

５、CUDAドライバのインストール（参考リンク）

CUDA 7.5 Downloads Archive | NVIDIA Developerからubuntu 64bitバージョンをダウンロード。

Linux -> x86_64 -> Ubuntu -> 14.04

#1.9~2GBくらいあるのでそれなりに時間はかかる。

#ダウンロードしたCUDA Repositoryのインストール

sudo dpkg -i cuda-repo-ubuntu1404_7.5-18_amd64.deb #dpkgはdebパッケージのインストール・アンインストールを行うコマンド(-i パッケージをインストール)

sudo apt-get update

#apt-getを使いCUDA Toolkitをインストール

sudo apt-get install cuda

#パスを追加

vi ~/.bashrc

　export CUDA_HOME=/usr/local/cuda-7.5

　export LD_LIBRARY_PATH=${CUDA_HOME}/lib64

　PATH=${CUDA_HOME}/bin:${PATH}

　export PATH

source .bashrc

sudo reboot

#CUDA SDK Samplesをホームディレクトリにコピーしてテストサンプルを実行

cuda-install-samples-7.5.sh  ~ 

cd ~/NVIDIA_CUDA-7.5_Samples 

cd 1_Utilities/deviceQuery 

make #サンプルコードのコンパイル

#makeすると以下のようなメッセージが出た

#CUDAバージョンの確認。

nvcc -V

#次のようなメッセージが出た。

正しく導入できているぽい。

cd ~/NVIDIA_CUDA-7.5_Samples/

make

cd bin/x86_64/linux/release/

./deviceQuery

#以下のようなメッセージがでた。

#最終行にPASSが出れば問題ないらしい。続いて

./bandwidthTest

#これも最終行にPASSが出れば問題ないらしい。

以下はdeep-photo-styletransferを動かすのに直接必要なもの。

６、torch7のインストール  科学技術計算の為の、機械学習ライブラリ

# clean old torch installation

rm -rf ~/torch

./clean.sh

git clone https://github.com/torch/distro.git ~/torch --recursive

cd ~/torch

# https://github.com/torch/distro : set env to use lua

TORCH_LUA_VERSION=LUA52 ./install.sh

# run luarocks WITHOUT sudo

luarocks install image

luarocks list

#実行するには

th #これだけ

#torchと一緒にmatio-ffiというのも入れろと書かれている。これはapt-getでインストールできる。

sudo apt-get install libmatio2

#loadcaffeも入れろと書いてある。

sudo apt-get install libprotobuf-dev protobuf-compiler

luarocks install loadcaffe

７、Octaveのインストール

#Octave は数値計算のためのプログラム言語

sudo apt-add-repository ppa:octave/stable

sudo apt-get update

sudo apt-get install octave

octave

と打てばoctaveが起動する。終了はexit

８、cuDDNのインストール (参考リンク)

cuDNNはDeep Learningに特化したライブラリらしい。Nvidiaに登録して 、ubuntu16.04向けをダウンロードした。インストールは解凍したファイルを所定の位置にコピーするだけ。

tar -zxf cudnn-7.5-linux-x64-v6.0.tgz #ダウンロードしたファイルの解凍

#解凍が終わるとcudaというディレクトリができる。

sudo cp cuda/include/cudnn.h /usr/local/cuda/include/

sudo cp cuda/lib64/* /usr/local/cuda/lib64/

準備が終わったのでdeep-photo-styletransferのソースコードをダウンロードし、中に入って追加でVGG-19をダウンロード。

sh models/download_models.sh

 30分ほどかかった。

#コンパイルしようとしたがエラーが出て失敗した。

make clean && make

 lua.hがないと言っているが、lua.hは

/home/user_name/torch/install/include/にあり、権限も与えてパスも通してある。原因がわからないので今はここまで

---------------------------------------------------------------------------------------------------- 

#以下のコマンドを覚えておく。

#GUI→CUIへ切り替え（Xウィンドウシステムの一時停止）

CTRL + ALT + F1

#CUI→GUIへ切り替え

ALT + F7

#Xサーバーの停止

sudo /etc/init.d/lightdm stop

Alt+PrintScreen を押しながら R -> S -> E -> I -> U -> B #入力し終わるとすぐ再起動がかかる。参考ページ 

#再起動

sudo reboot

#シャットダウン

sudo shutdown -h now

まずメインマシンのmac proにてubuntu.isoをコピーする該当ディスクをMS-DOSフォーマットで初期化する。基本的に元には戻せないので初期化は細心の注意を払って行う。容量が16GBあればUSBメモリでもいいし、余ったHDDやSSDでもOK。ただしubuntu.iso保存先のディスクはパーティション組んでても全パーティションの中身が消えるので気をつけること。

ターミナルで以下のコマンドを打つとファイルリストを表示してくれる。

diskutil list

disk9が該当するディスク。SSDを2.5インチケースに入れてUSBでmac proに接続している。

このディスクをアンマウントする。

diskutil unmountDisk /dev/disk9

ubuntu.isoのダウンロードパスに移動。

cd ~/Downloads/

以下のコマンドでubuntu.isoをdisk9に書き込ませる。このコマンドでターゲット先のデータは全て消えるので要注意。

sudo dd if=ubuntu-16.04.2-desktop-amd64.iso of=/dev/disk9 bs=4028

こんなメッセージが出て終わる。

385845+1 records in

385845+1 records out

1554186240 bytes transferred in 402.601338 secs (3860360 bytes/sec)

 USB2.0接続というのもあるが、10分くらいかかった。すぐに終わったらエラーの可能性もあるので注意して進めるべし。正常に終わるとmacで読み込めないMS-DOSフォーマットディスクがマウントしてくる。NTSFは読み込めないのでエラーが出るが、無視を選択して進める。

RNA seq 非モデル生物の解析

ゲノム情報がない場合、まず最初にアセンブリから始めることになる。RNAのアセンブリはTrinityが有名なのでtrinityを使ってみる。練習データとして、SRAに登録されているレタスのRNA seqデータを使う。

http://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA004542

36サンプルある。スクリプト組むほどではないので、手動でfastqを順番にダウンロードする。ペアでなくフラグメントデータぽい。

1、アセンブル

trinityでアセンブルを試みる。 

フラグメントデータなら

Trinity --seqType fq --single 1.fq --max_memory 10G --CPU 18

　--single；ペアードエンドシーケンスでなくフラグメントシーケンスデータの時

　--seqType；fastqならfq

　--max_memory ピークメモリー使用量の上限を指定

　--CPU : 使用するCPU数

ペアードエンドデータなら

Trinity --seqType fq --left R1.fq --right R2.fq --max_memory 10G --CPU 18

という感じに行う。R1.fqとR2.fqがペアードエンドデータになる。

複数のペアードエンドデータなら

(例えば下の３データ)

1_R1.fq　x　1_R2.fq

2_R1.fq　x　2_R2.fq

3_R1.fq　x　3_R2.fq

fastqを, (カンマ) でつないで以下のようにランする。

Trinity --seqType fq --left 1_R1.fq,2_R1.fq,3_R1.fq --right 1_R2.fq,2_R2.fq,3_R2.fq --max_memory 35G --CPU 18

ただし、今回のSRAデータを対象に実行したところ、ファイル名が認識できずスタートできなかった。原因を探ったところ以下のような記事が出てきた。

ヘッダーの仕様が古いせいっぽい。試しに割と新しいMiSeqのシーケンスデータを使ったところ、ちゃんとアセンブルできた。まあ下の記事のようにシェルかperlなどを使ってヘッダーを切り替えることのはルールさえわかれば簡単だろうけど...

How to change fastq reads header for running Trinity on them?

例えば

ヘッダーをオリジナルの

@DRR057317.1 MG00HS09:691:C8JNWACXX:7:1101:1476:1994 length=51

とかから

@M02077:53:000000000-ANJF9:1:1101:18340:000$a 1:N:0:15";

に変え（$aカウンタで数値を１リードずつ増やすすスクリプトで処理）、さらに３行目は新しいfastqの定義に従い+に変えて実行した。すると正常にランできるようになった。

 ただし今回はde novo RNA seqのワークフローを確認するために進めているだけなので、Trinityを使わずアセンブラのrnaspadesを使うことにする。まずはアセンブル前に解凍する。

bunzip2 -c DRR057309.fastq.bz2 > DRR057309.fastq

全サンプル順次解凍する。

rnaspadesでアセンブルする。以下のコマンドを打った。

rnaspades.py -t 30 -k auto -o all_data -s DRR057059.fastq -s \
DRR057060.fastq -s DRR057061.fastq -s DRR057062.fastq -s \
DRR057063.fastq -s DRR057064.fastq -s DRR057065.fastq -s \
DRR057066.fastq -s DRR057067.fastq -s DRR057068.fastq -s \
DRR057069.fastq -s DRR057070.fastq -s DRR057071.fastq -s \
DRR057072.fastq -s DRR057073.fastq -s DRR057074.fastq -s \
DRR057075.fastq -s DRR057076.fastq -s DRR057077.fastq -s \
DRR057078.fastq -s DRR057079.fastq -s DRR057080.fastq -s \
DRR057081.fastq -s DRR057082.fastq -s DRR057083.fastq -s \
DRR057084.fastq -s DRR057085.fastq -s DRR057086.fastq -s \
DRR057087.fastq -s DRR057088.fastq -s DRR057089.fastq -s \
DRR057090.fastq -s DRR057091.fastq -s DRR057092.fastq -s \
DRR057093.fastq -s DRR057094.fastq

全部で36サンプル。n=1、サイズは全て1GBくらい。メモリは200GBほどしかない。うまくできるか？と思いながら待っていたがダメだった。k-merを自動から指定値に変更する。リードは50bpの長さがある。rnaspadesは複数のk-merを指定できないので、奇数のk=21にして実行した。

rnaspades.py -t 30 -k 21 -o all_data_k21 -s DRR057059.fastq -s \
DRR057060.fastq -s DRR057061.fastq -s DRR057062.fastq -s \
DRR057063.fastq -s DRR057064.fastq -s DRR057065.fastq -s \
DRR057066.fastq -s DRR057067.fastq -s DRR057068.fastq -s \
DRR057069.fastq -s DRR057070.fastq -s DRR057071.fastq -s \
DRR057072.fastq -s DRR057073.fastq -s DRR057074.fastq -s \
DRR057075.fastq -s DRR057076.fastq -s DRR057077.fastq -s \
DRR057078.fastq -s DRR057079.fastq -s DRR057080.fastq -s \
DRR057081.fastq -s DRR057082.fastq -s DRR057083.fastq -s \
DRR057084.fastq -s DRR057085.fastq -s DRR057086.fastq -s \
DRR057087.fastq -s DRR057088.fastq -s DRR057089.fastq -s \
DRR057090.fastq -s DRR057091.fastq -s DRR057092.fastq -s \
DRR057093.fastq -s DRR057094.fastq

今度は無事解析が終わった。解析時間は６時間くらいだった。できたtranscripts.fastaの中を見てみる。

tail transcripts.fasta 
>NODE_134906_length_32_cov_12.4528_g133128_i0 CCCCTCGCCCATCCTCCTTTCTGTCGAATGTC >NODE_134907_length_31_cov_11.2115_g133129_i0 TTCCAAAGAATCCGAATCTTCGAATTCAAGT >NODE_134908_length_31_cov_4.61538_g133130_i0 CGACACAAAGACAGTAACCAATCAATGGACT >NODE_134909_length_31_cov_3.28846_g133131_i0 TACTGGTCTCATCTCTGTATGAAGGTGGGTG >NODE_134910_length_31_cov_2.63462_g133132_i0 TGATCGAAAAAATGGTGTTAGTCTCATCTGT

13万4910分子種ある。解析が重くなるので短かすぎるリードを除く。grep -n "length_99" transcripts.fasta

で99bpの行番号を特定し、

head -n transcripts.fasta > transcripts100over.fasta

として100bp以上を保存。nは99bpが始まる1つ前の行番号。spadesの出力はサイズ順なのでこれて過不足ない。

100bp以上に限定すると分子種は134910から86450に減った。流れを確認するトレーニングなので、とりあえずこれでマッピング解析を進める。

2、マッピング

#index作成

bowtie2-build -f transcripts_more_than_100bp.fa transcripts_100bp

#mapping

tophat -p 8 -o tophat_DRR057059 --no-novel-juncs transcripts_100bp DRR057059.fastq &
.
.
.
tophat -p 8 -o tophat_DRR057094 --no-novel-juncs transcripts_100bp DRR057094.fastq &

省略するが、上のように合計36サンプルを順次解析した。

コア数12、メモリ64GBのmac proでは４ジョブ投げるとCPUが100%に張り付いた。やはり重い。４つずつ実行したらそれだけで半日かかった。

3、cufflinksによるFPKM計算

cufflinksを使う。cufflinksは各遺伝子のリード数をカウントしてFPKMを算出したり、データを統合して比較することができるプログラム。

まずtophatのマッピングデータを使い発現量をカウントする。

cufflinks -o  cufflinks_output_0h tophat_DRR057059/accepted_hits.bam & cufflinks -o   cufflinks_output_2h tophat_DRR057060/accepted_hits.bam & cufflinks -o   cufflinks_output_4h tophat_DRR057061/accepted_hits.bam & cufflinks -o   cufflinks_output_6h tophat_DRR057062/accepted_hits.bam & cufflinks -o   cufflinks_output_8h tophat_DRR057063/accepted_hits.bam & cufflinks -o  cufflinks_output_10h tophat_DRR057064/accepted_hits.bam & cufflinks -o  cufflinks_output_12h tophat_DRR057065/accepted_hits.bam & cufflinks -o  cufflinks_output_14h tophat_DRR057066/accepted_hits.bam & cufflinks -o  cufflinks_output_16h tophat_DRR057067/accepted_hits.bam & cufflinks -o  cufflinks_output_18h tophat_DRR057068/accepted_hits.bam & cufflinks -o  cufflinks_output_20h tophat_DRR057069/accepted_hits.bam & cufflinks -o  cufflinks_output_22h tophat_DRR057070/accepted_hits.bam & cufflinks -o  cufflinks_output_24h tophat_DRR057071/accepted_hits.bam & cufflinks -o  cufflinks_output_26h tophat_DRR057072/accepted_hits.bam & cufflinks -o  cufflinks_output_28h tophat_DRR057073/accepted_hits.bam & cufflinks -o  cufflinks_output_30h tophat_DRR057074/accepted_hits.bam & cufflinks -o  cufflinks_output_32h tophat_DRR057075/accepted_hits.bam & cufflinks -o  cufflinks_output_34h tophat_DRR057076/accepted_hits.bam & cufflinks -o  cufflinks_output_36h tophat_DRR057077/accepted_hits.bam & cufflinks -o  cufflinks_output_38h tophat_DRR057078/accepted_hits.bam & cufflinks -o  cufflinks_output_40h tophat_DRR057079/accepted_hits.bam & cufflinks -o  cufflinks_output_42h tophat_DRR057080/accepted_hits.bam & cufflinks -o  cufflinks_output_44h tophat_DRR057081/accepted_hits.bam & cufflinks -o  cufflinks_output_46h tophat_DRR057082/accepted_hits.bam & cufflinks -o  cufflinks_output_48h tophat_DRR057083/accepted_hits.bam & cufflinks -o  cufflinks_output_50h tophat_DRR057084/accepted_hits.bam & cufflinks -o  cufflinks_output_52h tophat_DRR057085/accepted_hits.bam & cufflinks -o  cufflinks_output_54h tophat_DRR057086/accepted_hits.bam & cufflinks -o  cufflinks_output_56h tophat_DRR057087/accepted_hits.bam & cufflinks -o  cufflinks_output_58h tophat_DRR057088/accepted_hits.bam & cufflinks -o  cufflinks_output_60h tophat_DRR057089/accepted_hits.bam & cufflinks -o  cufflinks_output_62h tophat_DRR057090/accepted_hits.bam & cufflinks -o  cufflinks_output_64h tophat_DRR057091/accepted_hits.bam & cufflinks -o  cufflinks_output_66h tophat_DRR057092/accepted_hits.bam & cufflinks -o  cufflinks_output_68h tophat_DRR057093/accepted_hits.bam & cufflinks -o  cufflinks_output_70h tophat_DRR057094/accepted_hits.bam &

計算結果をcuffmergeでmergeしてgtfファイルを作る。そのために、上で行なった36個の解析結果のパスを書いたテキストファイルを準備し、cuffmergeにそれを読み込ませる。ファイル名はmerge.txtとし、以下のように記載した。

>cat merge.txt #中身をチェック
./cufflinks_output_0h/transcripts.gtf ./cufflinks_output_2h/transcripts.gtf ./cufflinks_output_4h/transcripts.gtf ./cufflinks_output_6h/transcripts.gtf ./cufflinks_output_8h/transcripts.gtf ./cufflinks_output_10h/transcripts.gtf ./cufflinks_output_12h/transcripts.gtf ./cufflinks_output_14h/transcripts.gtf ./cufflinks_output_16h/transcripts.gtf ./cufflinks_output_18h/transcripts.gtf ./cufflinks_output_20h/transcripts.gtf ./cufflinks_output_22h/transcripts.gtf ./cufflinks_output_24h/transcripts.gtf ./cufflinks_output_26h/transcripts.gtf ./cufflinks_output_28h/transcripts.gtf ./cufflinks_output_30h/transcripts.gtf ./cufflinks_output_32h/transcripts.gtf ./cufflinks_output_34h/transcripts.gtf ./cufflinks_output_36h/transcripts.gtf ./cufflinks_output_38h/transcripts.gtf ./cufflinks_output_40h/transcripts.gtf ./cufflinks_output_42h/transcripts.gtf ./cufflinks_output_44h/transcripts.gtf ./cufflinks_output_46h/transcripts.gtf ./cufflinks_output_48h/transcripts.gtf ./cufflinks_output_50h/transcripts.gtf ./cufflinks_output_52h/transcripts.gtf ./cufflinks_output_54h/transcripts.gtf ./cufflinks_output_56h/transcripts.gtf ./cufflinks_output_58h/transcripts.gtf ./cufflinks_output_60h/transcripts.gtf ./cufflinks_output_62h/transcripts.gtf ./cufflinks_output_64h/transcripts.gtf ./cufflinks_output_66h/transcripts.gtf ./cufflinks_output_68h/transcripts.gtf ./cufflinks_output_70h/transcripts.gtf  

mergeコマンドの実行。

cuffmerge -o all_merged merge.txt

-o　出力ディレクトリ名

gtfファイルが出力される。

4、発現変動遺伝子の分析

cufflinksのコマンドcuffdiffを使い発現変動を調べる。

cuffmergeの出力結果をinputとし、以下のようにコマンドを実行。

> cuffdiff -o cuffdiff_out -p 20 all/merged.gtf -L 0h,2h,4h,6h,8h,10h,12h,14h,16h,18h,20h,22h,24h,26h,28h,30h,32h,34h,36h,38h,40h,42h,44h,46h,48h,50h,52h,54h,56h,58h,60h,62h,64h,66h,68h,70h \
tophat_DRR057059/accepted_hits.bam tophat_DRR057060/accepted_hits.bam \
tophat_DRR057061/accepted_hits.bam tophat_DRR057062/accepted_hits.bam \
tophat_DRR057063/accepted_hits.bam tophat_DRR057064/accepted_hits.bam \
tophat_DRR057065/accepted_hits.bam tophat_DRR057066/accepted_hits.bam \
tophat_DRR057067/accepted_hits.bam tophat_DRR057068/accepted_hits.bam \
tophat_DRR057069/accepted_hits.bam tophat_DRR057070/accepted_hits.bam \
tophat_DRR057071/accepted_hits.bam tophat_DRR057072/accepted_hits.bam \
tophat_DRR057073/accepted_hits.bam tophat_DRR057074/accepted_hits.bam \
tophat_DRR057075/accepted_hits.bam tophat_DRR057076/accepted_hits.bam \
tophat_DRR057077/accepted_hits.bam tophat_DRR057078/accepted_hits.bam \
tophat_DRR057079/accepted_hits.bam tophat_DRR057080/accepted_hits.bam \
tophat_DRR057081/accepted_hits.bam tophat_DRR057082/accepted_hits.bam \
tophat_DRR057083/accepted_hits.bam tophat_DRR057084/accepted_hits.bam \
tophat_DRR057085/accepted_hits.bam tophat_DRR057086/accepted_hits.bam \
tophat_DRR057087/accepted_hits.bam tophat_DRR057088/accepted_hits.bam \
tophat_DRR057089/accepted_hits.bam tophat_DRR057090/accepted_hits.bam \
tophat_DRR057091/accepted_hits.bam tophat_DRR057092/accepted_hits.bam \
tophat_DRR057093/accepted_hits.bam tophat_DRR057094/accepted_hits.bam  

-o　write all output files to this directory 

-L　 comma-separated list of condition labels

-p　number of threads used during quantification

gtfファイルも相対パスで指定している。これが一番重くて４日かかった。クラミドモナスのtime courseサンプルのデータでは1時間くらいで終わったが、サンプル数が増えると処理時間が指数関数的に増えるのかもしれない。仕方がないので別のことをしながら待った。

解析が終わると指定したcuffdiff_out/に上のようなファイルができる。データベースのdbファイルは25GBもある。36サンプルもあると流石にでかい。

Rを起動して前回と同じくcummeRbundで解析してみる。

Rをターミナルで起動

library(cummeRbund) #ライブラリ読み込み
cuff <- readCufflinks("cuffdiff_out") #cufflinks解析結果のディレクトリ全体を読み込み
cuff #cuffコマンドを打って正常に読み込まれたか確認
　CuffSet instance with: 　
　36 samples 　
  44858 genes 　
  44858 isoforms 　
  44858 TSS 　
  0 CDS 　
  28260540 promoters 　
  28260540 splicing 　
  0 relCDS 

csBoxplot(genes(cuff))　#genesの分布をbox plotで表示

genediffオブジェクトに差のあるgenesを読み込む。

genediff <- diffData(genes(cuff)) 

統計的に優位なものだけ取ってくる。(significant列がyesのgene)

genediff2 <- genediff[(genediff$significant == 'yes'),]

続けてk.means法で分類。

genediffdataID2 <- genediff2[,1]
genediffdataID2 <- t(genediffdataID2)
uesakaGenes2 <-getGenes(cuff, genediffdataID2) #ここが重たい
k.means <-csCluster(uesakaGenes2, k=4) #k=300-400を越えるとRがクラッシュした
csClusterPlot(k.means)  #グラフを描画

genediff3 <- genediff[((genediff$status == 'OK')& (genediff$significant == 'yes')),]
genediffdataID3 <- genediff3[,1]
genediffdataID3 <- t(genediffdataID3)
uesakaGenes3 <-getGenes(cuff, genediffdataID3) 
k.means <-csCluster(uesakaGenes3, k=20) 
csClusterPlot(k.means)  

大量のデータがあるのでk.means法では適切な分類は難しいように思える。

やり方を変え、特定の遺伝子と同じ発現パターンをもつ遺伝子を探すことを考える。

そのため、まず全遺伝子のfpkmを出力し、変動パターンを確認して自分が興味ある変動パターンを示している遺伝子を見つけ、それと似た発現パターンの遺伝子をサルベージすることを考える。

全遺伝子のデータだけを取り出す

genes <- genes(cuff) 
is(genes)
fpkm <- fpkmMatrix(genes)  
write.table(fpkm, file="fpkm.csv", sep=",") #csv形式で保存

excelでcsvを開き、グラフを見て興味のある変動パターンを示す遺伝子を目視で調べていく。

例えばグラフの青の遺伝子は振動が徐々になくなっていっている。これは明暗条件からLL条件に移したサンプルと論文に書いてあるので、明暗の条件がなくなったことで周期を失う遺伝子と推測される。これと似た発現パターンをもつ遺伝子をcummeRbundで提供されているコマンドを使って絞り込んでいく。

my.similar <- findSimilar(cuff, "XLOC_000003", n = 20)

上のコマンドの意味だが、マニュアルにはFor example, if you wanted to find the 20 genes most similar to "PINK1", you could do the following: と書いてある。XLOC_000003と似た発現パターンを示すtop20遺伝子をfindSimilarコマンドで探してmy.similarに代入している。

散布図を描く

my.similar.expression <- expressionPlot(my.similar, logMode = T, showErrorbars = F) 
my.similar.expression

振動が徐々に消えていくようなパターンの発現パターンの遺伝子top20を絞り込めた（あくまでXLOC_000003に似た発現パターンの遺伝子を探しただけなことに注意）。

非モデル生物のRNA seqについては、様々なアプローチが考えられると思います。下のリンクでもう少し考えています。

結論としては、近縁種のゲノム情報が利用できてもde novo assemblyの方がorf検出率ははるかに上です。de novo assemblyしてcDNA配列を得て進めて行くことになるのではないでしょうか？

私が試みた限り、同じ属のゲノム情報を使いcufflinksでゲノムガイドアセンブルしても、ハウスキーピング遺伝子が10ほどしか検出できませんでした。

Rの基本的なコマンド 忘備録

オブジェクトはgenediffという名前とする。

nrow(genediff) #nrowは行数をカウント

[1] 796785

ncol(genediff) #ncolは列数をカウント

[1] 11

genediff$gene_id # ”data$列名” で参照したい列を取り出す。

table(genediff$gene_id) #要素の数をカウントするtable

#tableコマンドは要素の数を集計して表示してくれる。例えば日本、アメリカ、2国別の集計データでアメリカが何人表に載っているか調べたいなら

table(people$国籍)

結果は

アメリカ 8

上のように出力される。 複数の要素について集計することもできる。これはとんでもなく実用的。

table(people$国籍,性別)

　日本　アメリカ

男　3　5

女　7　3

つまり行列形式で出力される。

１行目、１列目の情報を表示するなら

rownames(genediff)　#rownamesは行名を表示

colnames(genediff)　#colnamesは列名を表示 

変数の型

is.character(genediff)　#is.characterでgenediffが文字列型ならTRUE

is.numeric(genediff)　#is.numericでgenediffが数字型ならTRUE

length(genediff)　#ベクトルの長さを調べるlength

mean(genediff)　#meanが平均

このデータでは　"引数は数値でも論理値でもありません。NA 値を返します" というエラーがでる。

max(genediff)　#maxで最大値

min(genediff)　#minで最小値

sum(genediff)　#sumで合計値

median(genediff)　#medianで中央値

この４つも純粋に数値だけのオブジェクトでないとエラーが出る。

RNA seq 其の４ おかしなデータを除いて再挑戦

30hのデータがおかしいので、それ以外の時間で分析を行う。tophatでのマッピングは其の２で行ったので、cuffdiffの定量だけやり直せばよい。これだけでかなり時間を短縮できる。

cuffdiff -o cuffdiff_except-_for_30h -p 20 Chlamydomonas_reinhardtii.v3.1.34.gtf -L 0h-1,0h-2,2h,4h,8h,12h,18h,24h,45h,60h 0h-1_out/accepted_hits.bam 0h-2_out/accepted_hits.bam 2h_out/accepted_hits.bam 4h_out/accepted_hits.bam 8h_out/accepted_hits.bam 12h_out/accepted_hits.bam 18h_out/accepted_hits.bam 24h_out/accepted_hits.bam 45h_out/accepted_hits.bam 60h_out/accepted_hits.bam

 Rを起動

library(cummeRbund) 

 cuff <- readCufflinks("cuffdiff_out")

genediff <- diffData(genes(cuff))

genediff2 <- genediff[(genediff$significant == 'yes'),]

write.table(genediff2, file="result3.csv", sep=",")

genediffdataID2 <- genediff2[,1]

genediffdataID2 <- t(genediffdataID2)

uesakaGenes2 <-getGenes(cuff, genediffdataID2)

k.means <-csCluster(uesakaGenes2, k=9)

k.means.plot <- csClusterPlot(k.means)

k.means.plot

k=9

k=6

k=4

k=9の 結果をよくよく見ると、30h以外にも24で発現が大きく落ちる遺伝子グループがある（１のウィンドウ）。もしかしたら30hの応答も生物的な応答かもしれない。　とりあえずk=9では十分なグループ分けができていない可能性もある。一度k=12で解析してみた。

k=12 

 time 4h、12h、24h、45h、60hでもそれぞれ大きく落ち込む遺伝子グループがある。0hはそれ以前の変動がわからないが0h特異的に発現が少ない遺伝子群かもしれない。

k=20まで上げてみよう。

K=20

 まだまだパターンがある気がするし、人間の目で見てかぶったパターンも増えたように思える。nitrate responseのサンプルなので,0hから発現がグンと上がっている遺伝子が気になる。k=20のウィンドウ１２の遺伝子をまとめた。

Nar1;2という硝酸同化系遺伝子が見事に応答している。非常に妥当な結果である。

ここで最初に立ち返ってt=30のデータはどうなのかを考える。例えばRNAの分解が進んだサンプルのrRNA/mRNA比率が代わり、特に低発現の遺伝子の発現パターンにばらつきが出るようになる。まあそれはそうだろう。分解スピードが同じと仮定し細胞内に10コピー存在するmRNAと100コピー存在するmRNAを考えると、先に0になるのは10コピーのmRNAの方だろう。このようなサンプルからRNAをとってシーケンスすると、分解してないサンプルと比べた時の比率が無限大になる。たとえFPKM 0.1と仮定して数値を加えたとしても、分解していないサンプルとの比率を正確に求めることはできない。また、分解しきっていなくても、低発現遺伝子の定量がばらつく要因はいくらでも考えられる。というのも、RNAの分解もシーケンス時の逆転写、PCR、シーケンスも（世の中のこと全て）確率的にしか進行しないため、初発の数値が低いまま進むとブレが大きくなるからだと解釈できる。例えると、サイコロを6回しか振らないと、6の目が３回出て、１は１回も出ない（オレすごい！）ようなことが起こるが、サイコロを6万回ふるとどの目も限りなく1/6に近い値がでる（オレ普通だった....）。それでも1/6に近づかないなら、サイコロが正６面体になっていないか重心がずれている可能性がある。

まずは色々な時間で散布図を書いてみる。

csScatter(genes(cuff),"X2h","X4h",smooth=T)

csScatter(genes(cuff),"X2h","X8h",smooth=T) 

 csScatter(genes(cuff),"X12h","X18h",smooth=T)

 csScatter(genes(cuff),"X18h","X45h",smooth=T) 

csScatter(genes(cuff),"X45h","X60h",smooth=T)

csScatter(genes(cuff),"X24h","X45h",smooth=T)

この中では45hが他と一番違うように感じる。さて、では30との比較だとどうなるか？

csScatter(genes(cuff),"X24h","X30h",smooth=T)

csScatter(genes(cuff),"X18h","X30h",smooth=T)

csScatter(genes(cuff),"X2h","X30h",smooth=T)

思ったほどばらつきはでかくなっていない。まあ対数表記だからか。しかし、どの時間と比較しても左下の方に弧を描くように描画されているプロット出る。説明できないけどアーティファクトっぽいなあ。

k=9の結果をヒートマップ化する。

k.means

出てきた数値をMiにコピー。エクセルに読み込んでフィルターをかけ、欲しいウィンドウだけ取り出してMiで保存。Rに再読み込み。

ウィンドウ１、8、3の遺伝子でヒートマップを出してみる。

csv <- read.table("1_3_8.txt" , header=F, sep="\t")

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

ウィンドウ6、7の遺伝子でヒートマップを出してみる。

csv <- read.table("1_3_8.txt" , header=F, sep="\t")

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

 思ったほどまとまらない。k=20の19と20のウィンドウだとどうだろうか？

 数が少ないのもあってさすがにこれは綺麗にまとまっている。

RNA seq 其の３ ヒートマップ作成

前回k-means法でクラスタリングするところまでやった。

k=9の結果が以下の通り。

k=9の結果についてヒートマップを書いて見る。

全部を１つのヒートマップに書くとしんどいので、分けて書く。まずは１、２、３のデータパターンが似たものをヒートマップ化する。

k.meansと打つと分析結果が表示されるので、

$silinfo

$silinfo$widths

                  cluster neighbor     sil_width

CHLREDRAFT_138936       1        2  0.4642057561

CHLREDRAFT_13257        1        2  0.4273940431

CHLREDRAFT_206075       1        2  0.4219634073

CHLREDRAFT_195343       1        2  0.4197930334

CHLREDRAFT_167270       1        2  0.4180063266

CHLREDRAFT_178353       1        2  0.4158867327

CHLREDRAFT_120099       1        2  0.4109362252

と書いてあるところをテキストエディタにコピーし、エクセルに読み込む。区分けはカンマかタブ。cluster neighbor のカラムでフィルターをかけ、1,2,3のclusterだけ取り出してcsv形式で保存 (1-3.csv）。

#読み込み

csv <- read.table("1-3.csv" , header=F, sep="\t")

前回と同じ

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

まだデータが多い。

１のクラスタだけで描いてみる。

ようやく見える解像度になってきた。しかしデータに問題がある可能性がある30hのデータの落ち込みで分類されたぽいクラスタ１なので、生物的な応答パターンを見えているか怪しい。30hのデータを除いてcufflinksの定量からやり直してみる。

 次回へ

RNA seq 其の２ マッピングから統計処理まで

基本的な流れは以下のリンクが参考になる。

github.com

今回はSRAに登録されているクラミドモナスのtime courseデータを使う。クラミドモナスはEnsemblにドラフトだがfastaもgtfも登録されているので非モデル生物より解析は楽になる。

クラミドモナスゲノムのダウンロードは

Chlamydomonas reinhardtii - Ensembl Genomes 35

のページの中ほどにある

Download DNA sequence (FASTA)

をクリックする。FTPで繋がり、任意のクロモソーム配列をダウンロードできる。

２、rRNAリードの除去

最初にrRNAにマッピングして、rRNAリードを除く。

rRNAのインデックス作成

bowtie2-build -f rRNA.fa rRNA

-f fastaファイルを指定

インデックス名をrRNA。この場合、rRNA~という６つのファイルが出力される。

rRNAへのマッピング

bowtie2 -x rRNA -1 R1.fastq -2 2.fastq -S rRNA_mapped.sam --un-conc rRNA_unmapped.fastq

-x　先ほどつけたインデックス名

-1 ペアードエンドデータのR1

-2 ペアードエンドデータのR2

-S マッピングされたデータ出力名（sam形式）

--un-conc マッピングされなかったfastqのペア。シングルなら--unを使う。上の場合はrRNA_unmapped 1.fastqとrRNA_unmapped 2.fastqが出力される。

ファイルサイズは2.03GBから2.02GBまでしか減らなかった。polyAから逆転写しているデータと思われる。ということで、このデータではrRNAの除去ステップは行わない。

３、クロモソームへのマッピング

faファイルとgtfファイルはEnsembl genomeの

http://plants.ensembl.org/info/website/ftp/index.html

のChlamydomonas reinhardtii からダウンロード。

インデックス作成

bowtie2-build -f Chlamydomonas_reinhardtii.v3.1.dna.toplevel.fa Chlamydomonas_reinhardtii.v3.1.dna.toplevel

ペアードエンドッピング

tophat -G Chlamydomonas_reinhardtii.v3.1.34.gtf -p 8 -o output --no-novel-juncs Chlamydomonas_reinhardtii.v3.1.dna.toplevel R1.fastq R2.fastq

-o　出力ディレクトリ名

-G gtfファイル。

-p　スレッド数。多すぎるとエラーになるので8くらいに留める。

--no-novel-juncs　アノテーションにない新規のジャンクションを破棄

最後にindex名、R1とR2のfastqを指定している。

４、FPKM算出

Cufflinksを使って計算する。

cufflinks -o output_FPKM -G genes.gtf -N --compatible-hits-norm mapping_0h/accepted_hits.bam 

-G 既存のアイソフォームのみの計算の指示

N --compatible-hits-norm 既存のアイソフォームにマッピングされるものだけでFPKMを計算

 -o　FPKMの出力ディレクトリ名

最後にtophatで出力したbamファイルを指定する。複数結果がある時は順番に実行する。

出力されたらFPKMのカラムをcutで抜き出す

cut -f 10 mapping_0h/accepted_hits.bam > 0h_FPKM.txt

cut -f 10 mapping_2h/accepted_hits.bam > 2h_FPKM.txt

 cut -f 10 mapping_0h/accepted_hits.bam > 4h_FPKM.txt

pasteコマンドでカラムを結合。

paste 0h_FPKM.txt 2h_FPKM.txt 4h_FPKM.txt > 0-4h_FPKM.txt

0-4h_FPKM.txtをエクセルかRに読み込んで経時変化をグラフ化する。

 5、Cuffdiffによる発現量に差がある遺伝子の抽出

cuffdiff -o cuffdiff_out -p 20 Chlamydomonas_reinhardtii.v3.1.34.gtf -L 0h-1,0h-2,2h,4h,8h,12h,18h,24h,30h,45h,60h 0h-1_out/accepted_hits.bam 0h-2_out/accepted_hits.bam 2h_out/accepted_hits.bam 4h_out/accepted_hits.bam 8h_out/accepted_hits.bam 12h_out/accepted_hits.bam 18h_out/accepted_hits.bam 24h_out/accepted_hits.bam 30h_out/accepted_hits.bam 45h_out/accepted_hits.bam 60h_out/accepted_hits.bam

ファイルが多いとかなりの時間がかかる。

下のような大量のファイルが出力される。

おまけ、 Go termによる分類

今回使用したChlamydomonasのgtfファイルにはEntrezのgene_id情報がある。これがあればすぐにGo解析ができる。それはChlamydomonasの情報はEMBLの運営するEnsembl Genomesからダウンロードしており、すでにEnsembl のaccession IDが付いているからである。それに対して、De novoのトランスクリプトーム解析ならアノテーションがないため、アノテーション、ID変換等を駆使していかなかればGo termをつけてカテゴリ分類することはできない。

AgBase GORetrieverを貼ったプレーンテキストをアップロードして、Ensembl accession IDを選択。チェックは必要に応じて外すが、全選択してjobをなげても問題ない。しばらく待つとab-deliminatedファイルが

ダウンロードできるリンクが発生する。クリックしてダウンロードする。ただし、GOの分類はそもそもまだまだ主観性の強い分類であり、分類が間違っている可能性がある。どの階層を選ぶかによって結果も変わってくるため、確度の高い分析を行うのは厳しい。あくまで大雑把にどのような応答になっているのか掴むための手法と捉えた方がよさそうである。

7、統計解析

Rで行う。まずは準備。

CummuneRbandのインストール

参照したページ。こんなことができるようです。

CummeRbund - An R package for persistent storage, analysis, and visualization of RNA-Seq from cufflinks output

#ターミナルでcufflinks解析結果ディレクトリの上に移動して、（つまり上の解析時そのままの位置で）  Rを起動(ターミナルでRと打ち、エンター）

R

library(cummeRbund) #ライブラリ読み込み

cuff <- readCufflinks("cuffdiff_out") #cufflinks解析結果のディレクトリ全体を読み込んでくれる

#これでCuffSet オブジェクトが生成された。

cuff

#正常に読み込まれていれば、cuffオブジェクトを分析しcuflinksの解析結果の概要が表示される。サンプルは0-60hのタイムコース実験の１１サンプル。

CuffSet instance with:

11 samples

14487 genes

14553 isoforms 0

TSS 0

CDS 0

promoters 0

splicing 0

relCDS

上のgenesの分布

csBoxplot(genes(cuff))

上の図作成は

box <- csBoxplot(genes(cuff))

print(box)

のようにやってもよい。この場合はまずcsBoxplot(genes(cuff))の結果をboxに代入し、それをpirntコマンドで画面に描画している。

#さらにpngで保存するなら

png("box.png")

plot(box)

dev.off()

これを同時にうつ。

分子生物学ではあまり使わないので、box plot（箱ひげ図）の見かたについて触れておく。（参考）

１、箱の中央付近のヨコ線 → データxの中央値。平均値ではない。また必ずしもboxの中央にあるわけでもない点に注意。

２、箱のヨコ線 → データxの第1四分位数（下側）と第3四分位数（上側）

分かりにくいので言い換える。箱の底辺の線の位置は、データの中央値と一番小さい値との間の中央値。つまり、中央値を最大値と考え、その最大値から下半分のデータについての中央値を求め、その位置を箱の低線の位置としている。この1/4番目の区切りの中央値を第1四分位数という。同様に箱の上辺の線の位置はデータの中央値と一番大きい値との間の中央値（=第3四分位数）。

３、箱の上下に向かって縦方向に伸びている線（ヒゲ） → データの最大値と最小値。

 ここで言う最大値と最小値はデータの真の最大値と最小値ではなくて外れ値を除外したもの。ヒゲの下方向の長さは

箱の底辺 - (箱の高さ）x 1.5

で出たの範囲で最も小さな値となる。

上の式を難しく言い換えると

第1四分位数-1.5 x (第3四分位数-第1四分位数)

ヒゲの上方向の長さは

第3四分位数+1.5*(第3四分位数-第1四分位数)

となる。その範囲に入らなかった値が外れ値として黒いプロットで表現されている。

箱ひげ図は品質管理で用いられたりしている。RNA seqの品質を管理するものとして捉えればよい。箱の位置が被らないようなデータは有意な差がある可能性が高い。30hのデータを見ると、中央値は他のデータと大差ないが箱の底辺は遥か下まで伸びている。言い換えると第1四分位数が他のサンプルと大きく異なっている。

isoformのbox plotを書く。

csBoxplot(isoforms(cuff))

アイソフォームでも30hのデータはおかしい。これが生物的な応答なのかサンプルの品質なのかは後で追及する。

csDensity(genes(cuff)) #ヒストグラム

２サンプル間の散布図

samples(object)

を使う。

samples(cuff)

   sample_index sample_name sample_name parameter value

1             1       X0h_1        <NA>      <NA>  <NA>

2             2       X0h_2        <NA>      <NA>  <NA>

3             3         X2h        <NA>      <NA>  <NA>

4             4         X4h        <NA>      <NA>  <NA>

5             5         X8h        <NA>      <NA>  <NA>

6             6        X12h        <NA>      <NA>  <NA>

7             7        X18h        <NA>      <NA>  <NA>

8             8        X24h        <NA>      <NA>  <NA>

9             9        X45h        <NA>      <NA>  <NA>

10           10        X60h        <NA>      <NA>  <NA>

なので、サンプルX2hとX4hを使って散布図を書くなら

csScatter(genes(cuff),"X2h","X4h",smooth=T)

csDendro(genes(cuff)) #樹形図

 csScatterMatrix(genes(cuff))

dispersionPlot(genes(cuff))

 csVolcanoMatrix(genes(cuff))

一部分を拡大

横軸はlog2(fold change)

縦軸はlog10(p value)

 #genediffというオブジェクトに差のある遺伝子データ（genes）を読み込む。

genediff <- diffData(genes(cuff))

#最初の部分だけを表示

head(genediff,3)

gene_id sample_1 sample_2 status value_1 value_2 log2_fold_change test_stat p_value q_value significant

1 AAB93440 X0h_1 X0h_2 NOTEST 1.34796 4.71785 1.80735 0 1 1 1 no

2 AAB93441 X0h_1 X0h_2 NOTEST 0.00000 0.00000 0.00000 0 2 1 1 no

3 AAB93442 X0h_1 X0h_2 NOTEST 0.00000 0.00000 0.00000 0 3 1 1 no

ここでRの基本的なコマンドを学習しておく。

#一番前にgene_idがあることに注目して差のある遺伝子のIDのカラムだけ抽出、genediffdataIDに代入。

genediffdataID <- genediff[,1] 

#中身をチェック

head(genediffdataID,10)

[1] "AAB93440" "AAB93441" "AAB93442"

[4] "AAB93443" "AAB93444" "AAB93445"

[7] "AAB93446" "AAB93447" "CHLREDRAFT_100005"

[10] "CHLREDRAFT_100008"

#headはカラムが少ないとこのように1つの行に3行分表示したりする。分かりにくい。

ファイルに保存してexcelで開いて見ると

write.table(genediffdataID, file="result.csv", sep=",")

こんな感じ。

 #データ解析には、縦のカラムではだめで行ベクトルに縦横変換する必要がある。縦横変換は

genediffdataID2 <- t(genediffdataID)

#ファイルを保存して

write.table(genediffdataID2, file="result.csv", sep=",")

#excelで開くと縦横変換されている。

 #ちなみにこれをターミナルからカッコよく開くなら

open result.csv #csvのデフォルトアプリで開く

#uesakaGenesオブジェクトにcuffオブジェクトからこのgenediffdataID2を持つ遺伝子の情報だけ代入

uesakaGenes <-getGenes(cuff, genediffdataID2)

#uesakaGenesをもとにヒートマップを書かせる

csHeatmap(uesakaGenes)

#数が多すぎてよくわからない。発現量が多い遺伝子だけに絞り込む。

#そのため差のある遺伝子のIDのカラムだけ抽出のところからやり直す。

genediff <- diffData(genes(cuff))

#統計的に優位なものだけ取ってくる。(significant列がyes)

genediff2 <- genediff[(genediff$significant == 'yes'),]

genediff2

#うてば何行残っているかわかる。2527行だけ残っていた。

#ファイルに保存してexcelで開くと

#右端の列がyesになっていることが確認できる。

あとは同じような流れ。

genediffdataID2 <- genediff2[,1]

genediffdataID2 <- t(genediffdataID2)

uesakaGenes2 <-getGenes(cuff, genediffdataID2)

 #Rでは k-means によるクラスタリングを行う関数が標準で用意されている。

#はじめてk-means法を聞く人が概要を掴むならiOSアプリの

アルゴリズム図鑑

• Moriteru Ishida
• 教育
• 無料
• 

 の中のk-meansの説明で素晴らしく分かりやすく説明されている。

k.means <-csCluster(uesakaGenes2, k=4)

k.means.plot <- csClusterPlot(k.means)

k.means.plot

#k=6にすると

k.means <-csCluster(uesakaGenes2, k=6) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

#k=8にすると

k.means <-csCluster(uesakaGenes2, k=8) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

#k=9にすると

k.means <-csCluster(uesakaGenes2, k=9) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

#k=12にすると

k.means <-csCluster(uesakaGenes2, k=12) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot 

これはさすがにやりすぎ。直感的には9が過不足なく綺麗に分けられている気がするがどうだろう。

30hで沈むデータがかなりあるのがわかる。box plotで30hだけおかしかったのはこれを反映しているのだろう。ぶっちゃけて言うと、おそらく30hのデータを除かないと正しい解析はできないだろう。

次回に続く。