読者です 読者をやめる 読者になる 読者になる

BLAST2GOでアノテーションをつける

basic版をここからダウンロード。インストールが終わったら立ち上げる。

統合TVの説明を参考にした。

リンク

 

 

BLAST2GOを起動したら、上のメニューからSTART -> Load Sequences -> Broseを選択し読み込むFASTAファイルを選択 -> Openボタンで読み込み。

 

基本的に左のアイコンから順に実行していく。解析が終わるとその色がつく。例えばblastはオレンジアイコンだが、blastがおわった配列はオレンジになってかつTagの列が紫になる。

 

 

 左端から3番目のBLAST

上のメニューからBLAST -> Run blast -> NCBI blast -> blastxを選択。データベースは"nr"。メールアドレスを入れて、そのほかはデフォルトのまま実行。

1配列で10分前後かかる。数千配列あるなら10日くらいかかるかもしれない。ちなみに

 local データベースを作成して高速化も可能らしい。

 ただし、nrのデータベースでもかなり容量がある。2016年頭にダウンロードしたときは解凍後200GBくらいにはなった。ダウンロードにも相当な時間がかかったし、保存場所も確保しないといけない。ローカルサーバーがあれば問題ないが、そうでなく時折しかblast解析しないのであればネットワーク版で問題ないと思う(ローカル版のデータベースの差分更新も既知の方法はありそう)。

 

 左端から4番目のINTERPRO SCAN

BLAST解析と並列化することが可能。役割はこちらのスレッドが分かりやすい。

やってることはモチーフ検索。blastでhypothetical proteinしかヒットしなくてもInterpro scanならみつかるかもしれない。統合TVの説明。

 

 上のメニューからInterpro -> Run interpro -> メールアドレスを記入してそのほかはデフォルトのまま実行。blastより短い傾向にあるがblastx同様1配列で数分かかるので、やはりこれも気長に終わるのを待つ。

 

 

  左端から5番目のmapping

blast情報に従ってgo termを関連付けする作業。

上のメニューからmapping-> Run mapping ->Runボタンをクリック。これだけ。

 

 

  左端から6番目のアノテーション

Go termは複数つくので、適切なものを選ぶ作業。

上のメニューからannotation-> Run annotation-> デフォルトのままRunボタンをクリック。

 

 

  左端から7番目のStatics

様々な情報を見ることができる。

 

 

エンリッチメント解析

特定のサンプルのみ変動している遺伝子群をGo termでまとめることができる。

メニューバーのAnalysis -> Enrichment Analysisを選択。

 

 

 

 

 

 

 

 

 

Waifuで拡大時の画質劣化を抑える

 

 

sudo apt update

sudo apt install cmake git libgtk2.0-dev pkg-config libavcodec-dev libavformat-dev libswscale-dev

sudo apt install python-dev python-numpy libtbb2 libtbb-dev libjpeg-dev libpng-dev libtiff-dev libjasper-dev libdc1394-22-dev

sudo apt install cmake freeglut3-dev libglew-dev git

git clone https://github.com/glfw/glfw && mkdir build && cd build && cmake ../glfw && make -j4 && sudo make install

wget https://github.com/opencv/opencv/archive/3.2.0.zip

unzip 3.2.0.zip

cd opencv-3.2.0/

cmake .

make

sudo make install

cd

git clone https://github.com/khws4vl/waifu2x-converter-glsl.git

cd waifu2x-converter-glsl/

CLC genomics workbench (7.0)でRNA seq解析

CLCで行う前提でワークフローを書いてみる。

 

アダプター配列の確認

アダプター配列はシーケンス->画像からのシグナル取得 -> ベースコールファイル(.bcl) -> FASTQ変換の過程で自動除去されるようだが、インサートが短いペアリードファイルなどでは3'側にアダプターが存在している可能性がある。例えば以下の

https://www.researchgate.net/post/PCR_primers_trimming_for_MiSeq

でそのことが論じられている。アダプター残存による影響は解析によって変わってくる。例えばCLCの資料の

https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2015_techsupport_session10.pdf

には、BWA MEMでマッピングした時には影響がないが(たぶんクリッピングされるから)、同じBWA でもショートリード向けのBWA Backtrackでは影響が多いと書かれている。アセンブルする場合も影響は大きいだろう。

 

 

アダプターが残存しているかCLCを使って調べてみる。

アダプター配列が存在しているか確認するには(オフィシャルサポートリンク)

 

 1、アダプターリストの作成

 2、fastqをダウンサンプリング

 3、モチーフ検索ツールでダウンサンプリングfastqファイル中のモチーフを検索

 4、アダプター配列が見つかることが確認できたら、アダプタートリミングツールでアダプターを除去

 

1、アダプターリストの作成

CLC7を起動

File -> New -> Trim motif list を選択。

ウィンドウが出てきたら下のAdd Rowボタンを押してアダプター配列情報を記入していく。Mismatchコストは初期値のままにしてFinish。パラメータの詳細はこちらを参考

f:id:kazumaxneo:20170418140342p:plain

 

アダプター配列分これを繰り返す。

最後にFile -> Save で保存。

 

 

2、fastqの一部データのサンプリング

・モチーフ検索ツールは1000リードまでしか対応していないみたいなので、FASTQのリードを1000だけランダムサンプリングする。手順は

・前もって読み込んだfastqファイルを右クリック -> Toolbox -> NGS core tool -> Sample reads

リード数を1000に指定してRun

10秒くらいでサンプリングは終わる。~ sampledというファイルが作られているはず。

 

 

3、モチーフ検索ツールでサンプリングしたfastqファイル中のモチーフを検索

・~ sampledファイルを右クリック -> Toolbox -> Classical Sequence Analysis -> General Sequence Analysis -> Motif search を選択。

・Motif Search TypeはMotif Listを選択。

・Motif listは先ほどで作成したリストを選択。

・Include negative strandにチェック

f:id:kazumaxneo:20170418141636p:plain

 

・次のページでCreate TablesとAdd anotations to sequencesにチェック

・resultはsaveを選択しfinish。

 ・解析が終わると、モチーフが見つかったリードを記載したリストができる。

 

 

4、アダプター配列が見つかるか確認

・~ sampledファイルの中身をダブルクリックで開き、右側のリストからAnnotation typeを選択する。分かりにくいが写真の青色のMotifというアイコンの右上にAnnotation typeというタブがある。

f:id:kazumaxneo:20170418142416p:plain

 

・青色のMotifボックスにチェックをつけるとモチーフが見つかった部位が矢印で表示される。今回は人工的に下の画像の矢印部位が100%マッチするようにモチーフ配列を設計した。

f:id:kazumaxneo:20170418144401p:plain

 

同じことをRでやるなら

library(qrqc)

fastq <- readSeqFile("input.fastq")

basePlot(fastq)

f:id:kazumaxneo:20170418154709j:plain

先頭だけATGC頻度が極端なのでアダプターと考えられる。

 

kmerKLPlot(fastq)  #6-mer時の塩基出現頻度

f:id:kazumaxneo:20170418155116j:plain

高頻出の配列がアダプターということになる。

 

k-mer長を変えるなら

fq <- readSeqFile("1.fq", kmer=TRUE, k=15) #k=15の時 

kmerKLPlot(fq) #グラフ描画

 

 

 

アダプター配列のトリミング

上の画像の部位のトリミングを行ってトリミングツールの挙動を確かめてみる。以下の手順で進める。

File -> New -> Trim Adapter list を選択。

・アダプター配列を順番に記載してsave。リストファイルを作る。(上のモチーフ検索リスト作成と同じ要領)

・NGS core tools -> Trim sequences を選択。

・いまアダプタートリミングだけ行いたいので、クオリティトリミングや5' or 3'末端トリミングのチェックを全て外す。先ほど作ったリストファイルを読み込ませ、Finishでラン開始。

 

終わったら、FASTQファイルのtrimmedというのができるので、開いてみる。

f:id:kazumaxneo:20170418145842p:plain

先ほどの青色矢印部分の配列が除去されている。

アダプターの上流側に配列がある場合、それもまとめて除去されている。5'側がまとめて除去されるので、index配列が分かればその上流も含めて全て除去できるのは都合がいい。ただし、配列の途中にアダプターが出てくる場合、強制的に5'側がトリムされる仕様は不都合が出てくるかもしれない。例えばインサートが短くてペアリードが調整したライブラリの末端までシーケンスされていた場合(特にsmall RNAシーケンスなど)、3'側にアダプター配列が出てくる可能性がある。その場合、下流側の配列を除く必要があり、強制的に5'側が削除されるのはまずい。

 

 

CLC以外のコマンドツールではどうか?

1、cutadaptの場合

 

@DRR057320.1 MG00HS11:593

TTAGGTTTCTACAAAATGAAGATTTCGAAAGTTTATCAAAACAAAGAATCT

+

@@@DDA=DHFFBDA:EFHIII@FHHHEG?+??CFEDEDHD:FHCG;DCGIF

上のリードの太字部分がトリム対象なら、cutadapt -gコマンドで

cutadapt -g TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #5'側トリミング

 

トリム後、アダプターの5'側も除去される。

@DRR057320.1 MG00HS11:593

AAAGTTTATCAAAACAAAGAATCT

+

G?+??CFEDEDHD:FHCG;DCGIF

こうなる。3'側をトリムするなら-gを-aに変えて

cutadapt -a TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #3'側トリミング

リードを見ると

@DRR057320.1 MG00HS11:593

TTAGGTTTC

+

@@@DDA=DH

今度はアダプターとその下流側が除去されている。

 

 

 

 

 

 

 

 

トリミング

詳細はこちらを参照した。 

CLC genomics workbench のトリムは累積のクオリティスコアが一定以上の領域が続いた場合に取り除くらしい。つまり1塩基だけクオリティスコアが悪くても必ずしもトリミングはされないアルゴリズムを用いている。

#アダプター配列の除去

toolboxNGS core tools → Trim sequences を選択。

 

de novo アセンブリ

trinityの結果と比較する。

Trinity --seqType fq --left S1_L001_R1.fastq --right S1_L001_R2.fastq --max_memory 10G --CPU 18

 

>行をカウントしてコンティグ数を見積もる。

Trinity 1656 contigs

CLC genomics workbench 7 1956 contigs

rnaspades 1952 contigs

 

数についてはTrinityが一番少なく、rnaspadesが一番多い。

長さの分布やlongest conntigを見るとどうか?

 

trinity

f:id:kazumaxneo:20170415141857p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142259p:plain

 

 

CLC genomics workbench 7(default condition)

f:id:kazumaxneo:20170415142525p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142214p:plain

 

 

 

 

rnaspades

rnaspades.py -t 30 -k auto -o data -1 1_R1.fq -2 1_R2.fq )

TOP10の長さ

Name length
NODE1 19799
NODE2 18969
NODE3 17335
NODE4 15793
NODE5 15760
NODE6 13466
NODE7 13232
NODE8 12802
NODE9 12596
NODE10 12065

 

 

top10の長さはrnaspadesとtrinityが優れている感じ。今回のデータはMiseq 301bpのペアデータである。他のデータでも同様の傾向が出るのかはわからない。

 

 

 

 

 

 

 

 

 

 

 

 

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

 ..(以下略)

 

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

 

 

 

2、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 #問題がないか探してくれる

 

 

 

3、初回の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行ずつ文が出てくるので、1行目は空白のままエンター。2行目でパス(桁の小さな数値は拒否される)、3行目は確認のためもう一度パスを入力。

#/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

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

 

 

4、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

 

#もう1つ作成

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

 options nouveau modeset=0

 

sudo update-initramfs -u

 

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

 

#準備ができたら再起動

sudo reboot

 

#再起動して1分くらい待つ。クラアント側のパソコン(今回は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(そのままエンターキー)5分くらいかかる。

#終わったら再起動。

 sudo reboot

 

#もう一度インストール

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

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

sudo reboot

 

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

nvidia-smi

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

f:id:kazumaxneo:20170410193300j:plain

 

 

5、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すると以下のようなメッセージが出た

f:id:kazumaxneo:20170410181956j:plain

 

#CUDAバージョンの確認。

nvcc -V

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

f:id:kazumaxneo:20170410182112j:plain

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

 

cd ~/NVIDIA_CUDA-7.5_Samples/

make

cd bin/x86_64/linux/release/

./deviceQuery

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

 

f:id:kazumaxneo:20170411124643j:plain

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

./bandwidthTest

f:id:kazumaxneo:20170411124649j:plain

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

 

 

 

 

 

 

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

 

6、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 #これだけ

f:id:kazumaxneo:20170410193503j:plain

 

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

sudo apt-get install libmatio2

 

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

sudo apt-get install libprotobuf-dev protobuf-compiler

luarocks install loadcaffe

 

 

7、Octaveのインストール

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

sudo apt-add-repository ppa:octave/stable

sudo apt-get update

sudo apt-get install octave

 

octave

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

 

 

 

8、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

f:id:kazumaxneo:20170411131638j:plain

 

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

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

 

 

 

 

 

 

 

 

 

 

 

 

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

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

#GUICUIへ切り替え(Xウィンドウシステムの一時停止)

CTRL + ALT + F1

#CUIGUIへ切り替え

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

f:id:kazumaxneo:20170403155204j:plaindisk9が該当するディスク。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を順番にダウンロードする。ペアでなくフラグメントデータぽい。trinityでアセンブルを試みる。

 

trinityのコマンドはフラグメントデータなら

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

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

 --seqType;fastqならfq

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

 --CPU : 使用するCPU数

1.2GBのシングルデータで1-2時間かかる。

 

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

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

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

 

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

(例えば下の3データ)

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カウンタで数値を1リードずつ増やすすスクリプトで処理)、さらに3行目は新しい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

 

今度は無事解析が終わった。解析時間は6時間くらいだった。できた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以上に限定すると788737行から643827行に減り、分子種は134910から86450に減った。練習なのでとりあえずこれでマッピング解析を進める。

 

#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では4ジョブ投げるとCPUが100%に張り付いた。やはり重い。4つずつ実行したらそれだけで半日かかった。いやたぶん軽い方なんだろう。

 

cufflinksでアセンブルする。

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ファイルが出力される。

これを使いcuffdiffでのDEG解析を行う。

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

 

 

f:id:kazumaxneo:20170406155216j:plain

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

 

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

R #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で表示

f:id:kazumaxneo:20170406154047j:plain

 

#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)  #グラフを描画

f:id:kazumaxneo:20170406190327j:plain

  

 

genediff3 <- genediff[*1,]

genediffdataID3 <- genediff3[,1]

genediffdataID3 <- t(genediffdataID3)

uesakaGenes3 <-getGenes(cuff, genediffdataID3) 

k.means <-csCluster(uesakaGenes3, k=20) 

csClusterPlot(k.means) 

f:id:kazumaxneo:20170406185428j:plain

 

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

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

 

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

 

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

genes <- genes(cuff)

is(genes)

 

 fpkm <- fpkmMatrix(genes)

 write.table(fpkm, file="fpkm.csv", sep=",") #csv形式で保存

 

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

f:id:kazumaxneo:20170406203848j:plain

例えばグラフの青の遺伝子は振動が徐々になくなっていっている。これは明暗条件から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

f:id:kazumaxneo:20170406204333j:plain

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

*1:genediff$status == 'OK')& (genediff$significant == 'yes'

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

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

 

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

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で中央値

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