mecobalamin’s diary

人間万事塞翁が馬

メモ: Dockerを使ってみる

dockerをwslで動くようにインストールした
mecobalamin.hatenablog.com
とりあえずhello-worldは動いたが、目的のコンテナが動いてくれない

もうちょっと情報集め

試しに他のコンテナが動くかどうかを調べたい

Docker Hubってのがあった
https://hub.docker.com/

ここのDockerコンテナ実行例を試した
Docker入門(第一回)~Dockerとは何か、何が良いのか~ | さくらのナレッジ

NginxのDockerイメージを使ってWebサーバーを立ててみましょう。

$ docker run --name some-nginx -d -p 8080:80 nginx

ウェブサーバーにアクセスする
IPアドレスの確認は次のコマンド

$ docker-machine.exe config
--tlsverify
--tlscacert="C:\\Users\\hogehoge\\.docker\\machine\\machines\\default\\ca.pem"
--tlscert="C:\\Users\\hogehoge\\.docker\\machine\\machines\\default\\cert.pem"
--tlskey="C:\\Users\\hogehoge\\.docker\\machine\\machines\\default\\key.pem"
-H=tcp://192.168.99.100:2376

ウェブブラウザのアドレスにhttp://192.168.99.100:8080でアクセスできた

おそらくdockerはちゃんと動いてると思う
問題は目的のコンテナの方か?

今動いているコンテナは次のコマンドで確認できた

$ docker stats
CONTAINER ID        NAME                CPU %               MEM USAGE / LIMIT     MEM %               NET I/O             BLOCK I/O           PIDS
6242a4f20125        some-nginx          0.00%               1.992MiB / 989.4MiB   0.20%               836B / 0B           0B / 0B             2

停止するときはCONTAINER IDとNAMEをそれぞれ次のコマンドに使って止められた

$ docker stop 6242a4f20125 
$ docker rm some-nginx

python3.xと2.xを切り替える

condaでインストールするパッケージで、
python2.xで動くのがある
仕組みがよくわかってないけどこのパッケージを
condaのコマンドでインストールすると
pythonを3.xから2.xに変えてしまうっぽい

なので2.x系で動くパッケージはpython2.xが動く環境を作成して
そっちにインストールする

コマンドを使うたびに環境を切り替えるのは面倒だが
pythonのバージョンが変わってしまうのは厄介

python 2.7の環境を作る
コマンドはこれ

$ conda create -n py27 python=2.7 anaconda

py27は環境の名前なので任意でいいはず(未確認)

インストールしたpython2.xの環境に変更する

$ source activate py27

コマンドプロンプトの前に(py27)と入るので使用中の環境がわかる

python3.xの環境に戻す

$ conda deactivate

以前は

$ source deactivate py27

を使っていたけど、warningが出るようになってるので変更。

この記事で使用
mecobalamin.hatenablog.com

14 June 2019追記
Powershellで同じことをしようとしたらcondaのインストールができず
そもそも対応してなかったらしい
conda · PyPI

minicondaを入れたけどPowerShellからconda使えず。。。
こっちも最近対応したばかりらしい
Anaconda がやっと PowerShell に公式対応した - Qiita
PowerShellで使うときは
[スタート] -> [Anaconda Powershell Prompt]
で起動する

python 2.7の環境の作り方は上記と同じ

$ conda create -n py27 python=2.7 anaconda

インストール後に環境の切り替え方法が出る

#
# To activate this environment, use
#
#     $ conda activate py27
#
# To deactivate an active environment, use
#
#     $ conda deactivate

追記ここまで

fioその2、ストレージのベンチマーク測定

前回の続き
mecobalamin.hatenablog.com

今回は測定結果
内蔵SSDSATA接続

$ bash ./fio-cdm_190518.sh /mnt/c/fio_tmp
fio-3.1
/mnt/c/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2339.000|   1743.000|
| 512K |   2223.000|   1579.000|
|   4K |    192.000|    199.000|
|4KQD32|    212.000|    204.000|

内蔵HDD、SATA接続

$ bash ./fio-cdm_190518.sh /mnt/e/fio_tmp
fio-3.1
/mnt/e/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2468.000|   1737.000|
| 512K |   2205.000|   1721.000|
|   4K |    214.000|    201.000|
|4KQD32|    223.000|    195.000|

外付けHDD、USB3.1 Gen1 typeC接続

$ bash ./fio-cdm_190518.sh /mnt/h/fio_tmp
fio-3.1
/mnt/h/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2187.000|     72.700|
| 512K |   2223.000|     35.600|
|   4K |    173.000|    183.000|
|4KQD32|    171.000|    195.000|

USBメモリUSB3.1接続

$ bash ./fio-cdm_190518.sh /mnt/usb/fio_tmp
fio-3.1
/mnt/usb/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2237.000|     20.800|
| 512K |   1967.000|   1034.000|
|   4K |    230.000|    360.000|
|4KQD32|    231.000|    251.000|

RAM disk

$ bash ./fio-cdm_190518.sh /mnt/r/fio_tmp
fio-3.1
/mnt/r/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2251.000|   1811.000|
| 512K |   2196.000|   1732.000|
|   4K |    257.000|    243.000|
|4KQD32|    232.000|    213.000|

外付けSSDUSB3.1 Gen1 TypeC接続

$ bash ./fio-cdm_190518.sh /mnt/s/fio_tmp
fio-3.1
/mnt/s/fio_tmp
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2156.000|    173.000|
| 512K |   2256.000|    161.000|
|   4K |    163.000|    176.000|
|4KQD32|    165.000|    184.000|

ちなみにCrystal Disk Markでは

f:id:mecobalamin:20190518035347p:plain
Crystal Disk Markの結果、外付けSSD

Crystal Disk Markに比べるとReadがとても早い
なんか設定が違うか。。。?
syncとlibaio?
もう一度試したけどlibaioはエラーが出て測定不可

どのストレージでもreadは同程度
writeはそれぞれ特徴がでてそう
読み込み中心のコマンド・スクリプトだと
ストレージの接続方法であまり差が出ないかも
Windows側と比較できないのが残念
もうちょい試す

21 May 2019追記
fioを実行しているときにwindowsのタスクマネージャーを見ると
内蔵SSDSATA接続の書き込み速度は数十MB/sec程度だった
fioの測定結果よりも遅い
Crystal Disk Mark実行時では数百MB/sec程度で
SSDの性能通りの値が出ている

fioで書き込むファイルのサイズを5gにすると
書き込み速度が数百MB/sec程度になった

 time bash ./fio-cdm_190520.sh /mnt/c/fio_tmp/
fio-3.13-47-gd549
/mnt/c/fio_tmp/
|      | Read(MB/s)|Write(MB/s)|
|------|-----------|-----------|
|  Seq |   2430.000|    450.000|
| 512K |   2450.000|    395.000|
|   4K |    206.000|    139.000|
|4KQD32|    248.000|    135.000|

書き込みの速度が似た感じになってきた
なんとなくfioではメモリへのアクセス速度を見ている気がする
メモリを使い切るとSSD/HDDへのアクセスが始まって
実際の速度が現れるのだろうか?
wslのI/Oの仕組みをよくわかってないからなんとも言えないが。。。

やたら早いfioの結果はあまり信じないほうが良さげ
アクセス速度がwindowsからの測定とwslからの測定で違う理由はよくわからないが、
wslの書き込みはCrystal Disk Markの結果を見たほうが良さそう

fioその1、ストレージのベンチマーク測定

以前windowsCrystalDiskMarkを使って
ストレージのベンチマークを測定した

今回はwslのコマンドのfioを使って測定する

以前書いた記事
mecobalamin.hatenablog.com
fioのについてもちょっと書いてある
mecobalamin.hatenablog.com

fioのバージョンは3.1

こちらのラッパーを利用しようとしたが
fioバージョン?の違いか出力フォーマットが変わっているっぽい
fioでCrystalDiskMarkっぽい計測を行うコマンドを作った - ぶていのログでぶログ

とはいえawk正規表現の両方、殆ど触ったことがない。。。
調べながら書き換え

正規表現はここを見ながら書き換えた
正規表現の基礎
正規表現で整数と実数にマッチさせる
正の整数と小数にマッチする正規表現
ただしwslのawkでは\dが働いていない様子
[0-9]に変えるとマッチしたのでこっちを使う

splitとgensubはawkの組み込み関数
The GNU Awk User's Guide - 組み込み関数
gawkで正規表現の後方参照 - ちぎっては投げるブログ


さて本題

fioの出力はこうなっている

Seq-Read: (g=0): rw=read, bs=(R) 1024KiB-1024KiB, (W) 1024KiB-1024KiB, (T) 1024KiB-1024KiB, ioengine=sync, iodepth=1
Seq-Write: (g=1): rw=write, bs=(R) 1024KiB-1024KiB, (W) 1024KiB-1024KiB, (T) 1024KiB-1024KiB, ioengine=sync, iodepth=1
Rand-Read-512K: (g=2): rw=randread, bs=(R) 512KiB-512KiB, (W) 512KiB-512KiB, (T) 512KiB-512KiB, ioengine=sync, iodepth=1
Rand-Write-512K: (g=3): rw=randwrite, bs=(R) 512KiB-512KiB, (W) 512KiB-512KiB, (T) 512KiB-512KiB, ioengine=sync, iodepth=
1
Rand-Read-4K: (g=4): rw=randread, bs=(R) 4096B-4096B, (W) 4096B-4096B, (T) 4096B-4096B, ioengine=sync, iodepth=1
Rand-Write-4K: (g=5): rw=randwrite, bs=(R) 4096B-4096B, (W) 4096B-4096B, (T) 4096B-4096B, ioengine=sync, iodepth=1
Rand-Read-4K-QD32: (g=6): rw=randread, bs=(R) 4096B-4096B, (W) 4096B-4096B, (T) 4096B-4096B, ioengine=sync, iodepth=32
Rand-Write-4K-QD32: (g=7): rw=randwrite, bs=(R) 4096B-4096B, (W) 4096B-4096B, (T) 4096B-4096B, ioengine=sync, iodepth=32
fio-3.1
Starting 8 processes
Seq-Read: Laying out IO file (1 file / 1024MiB)

Seq-Read: (groupid=0, jobs=1): err= 0: pid=27737: Fri May 17 17:33:21 2019
   read: IOPS=3908, BW=3908MiB/s (4098MB/s)(1024MiB/262msec)
    clat (usec): min=205, max=4123, avg=253.44, stdev=139.86
     lat (usec): min=206, max=4125, avg=253.55, stdev=139.93
    clat percentiles (usec):
     |  1.00th=[  208],  5.00th=[  212], 10.00th=[  215], 20.00th=[  221],
     | 30.00th=[  231], 40.00th=[  233], 50.00th=[  235], 60.00th=[  237],
     | 70.00th=[  247], 80.00th=[  255], 90.00th=[  273], 95.00th=[  322],
     | 99.00th=[  725], 99.50th=[  799], 99.90th=[  832], 99.95th=[ 4113],
     | 99.99th=[ 4113]
  lat (usec)   : 250=73.93%, 500=23.93%, 750=1.27%, 1000=0.78%
  lat (msec)   : 10=0.10%
  cpu          : usr=0.00%, sys=101.53%, ctx=0, majf=0, minf=289
  IO depths    : 1=100.0%, 2=0.0%, 4=0.0%, 8=0.0%, 16=0.0%, 32=0.0%, >=64=0.0%
     submit    : 0=0.0%, 4=100.0%, 8=0.0%, 16=0.0%, 32=0.0%, 64=0.0%, >=64=0.0%
     complete  : 0=0.0%, 4=100.0%, 8=0.0%, 16=0.0%, 32=0.0%, 64=0.0%, >=64=0.0%
     issued rwt: total=1024,0,0, short=0,0,0, dropped=0,0,0
     latency   : target=0, window=0, percentile=100.00%, depth=1
以下省略

見方がよくわかってないけど

Seq-Read: (groupid=0, jobs=1): err= 0: pid=27737: Fri May 17 17:33:21 2019
   read: IOPS=3908, BW=3908MiB/s (4098MB/s)(1024MiB/262msec)

この部分から(4098MB/s)を読み取れればいいかと思う

オリジナルではこうなっている

split(gensub(/bw=([0-9.]+)(([KM]?)B\/s,)?/,"\\1 \\3", "g", i), a);

書き換えた結果がこれ

split(gensub(/\((([1-9]+[0-9]*|0)(\.[0-9]+)?)([KM]?[^i]?B\/s)?\)/,"\\1 \\2", "g", i), a);

あと他にも書き換えたところがある
ioengineでlibaioを指定すると初期化に失敗してるみたい
syncを指定している

スクリプト全体はこうなる
gitでforkするのがいいんだろうけど
やり方が調べてなくてわからないので。。。

ファイル名をfio-cdm_190518.shで保存
コード全体はこんな感じ

#! /bin/bash

TARGET="$1"
FileName=".fio-diskmark"

fio -version
echo ${TARGET}

fio2cdm() {
    awk '
    /^Seq-Read:/          {getline;if($1~/^read/) {seqread =$4}}
    /^Seq-Write:/         {getline;if($1~/^write/){seqwrite=$4}}
    /^Rand-Read-512K:/    {getline;if($1~/^read/) {rread512 =$4}}
    /^Rand-Write-512K:/   {getline;if($1~/^write/){rwrite512=$4}}
    /^Rand-Read-4K:/      {getline;if($1~/^read/) {rread4 =$4}}
    /^Rand-Write-4K:/     {getline;if($1~/^write/){rwrite4=$4}}
    /^Rand-Read-4K-QD32:/ {getline;if($1~/^read/) {rread4qd32 =$4}}
    /^Rand-Write-4K-QD32:/{getline;if($1~/^write/){rwrite4qd32=$4}}

    function n(i) {
      split(gensub(/\((([1-9]+[0-9]*|0)(\.[0-9]+)?)([KM]?[^i]?B\/s)?\)/,"\\1 \\2", "g", i), a);
      s = a[1]; u = a[2];
      if(u == "K") {s /= 1024}
      if(u == "")  {s /= 1024 * 1024}
      return s;
    }
    END {
      print ("|      | Read(MB/s)|Write(MB/s)|");
      print ("|------|-----------|-----------|");
      printf("|  Seq |%11.3f|%11.3f|\n", n(seqread),   n(seqwrite));
      printf("| 512K |%11.3f|%11.3f|\n", n(rread512),  n(rwrite512));
      printf("|   4K |%11.3f|%11.3f|\n", n(rread4),    n(rwrite4));
      printf("|4KQD32|%11.3f|%11.3f|\n", n(rread4qd32),n(rwrite4qd32));
    }
    '
}

trap "rm -f ${TARGET}/${FileName}" 0 1 2 3 9 15

# see. http://www.winkey.jp/article.php/20110310142828679
cat <<_EOL_ | fio - | fio2cdm
[global]
ioengine=sync
iodepth=1
size=1g
direct=1
runtime=60
directory=${TARGET}
filename=${FileName}

[Seq-Read]
bs=1m
rw=read
stonewall

[Seq-Write]
bs=1m
rw=write
stonewall

[Rand-Read-512K]
bs=512k
rw=randread
stonewall

[Rand-Write-512K]
bs=512k
rw=randwrite
stonewall

[Rand-Read-4K]
bs=4k
rw=randread
stonewall

[Rand-Write-4K]
bs=4k
rw=randwrite
stonewall

[Rand-Read-4K-QD32]
iodepth=32
bs=4k
rw=randread
stonewall

[Rand-Write-4K-QD32]
iodepth=32
bs=4k
rw=randwrite
stonewall
_EOL_

測定結果は次に。

メモ: ヒアドキュメント

ストレージのベンチマークを測定するコマンドにfioがある
調べてみると実行時のオプションとか
結果の出力結果とか、とても情報が多い
出力結果をうまく処理してくれるラッパーを作ってくれてる人がいるのだが、
シェルスクリプトで見たことない書き方があった

fioでCrystalDiskMarkっぽい計測を行うコマンドを作った
fioでCrystalDiskMarkっぽい計測を行うコマンドを作った - ぶていのログでぶログ

fio-cdm/fio-cdm
fio-cdm/fio-cdm at master · buty4649/fio-cdm · GitHub

このあたり

# see. http://www.winkey.jp/article.php/20110310142828679
cat <<_EOL_ | fio - | fio2cdm
[global]
ioengine=libaio
iodepth=1
size=1g
direct=1
runtime=60
directory=${TARGET}
filename=.fio-diskmark

~中略~

[Rand-Write-4K-QD32]
iodepth=32
bs=4k
rw=randwrite
stonewall
_EOL_

ヒアドキュメントというらしい。
使い勝手が良さそうなのでリンクのメモ。

ヒアドキュメント - Wikipedia

ヒアドキュメント(別の呼び方としてヒア文字列、heredocなど)は、文字列リテラルを、シェルスクリプトプログラミング言語ソースコード中に埋め込むための1つの方法である。

Here document - Wikipedia

In computing, a here document (here-document, here-text, heredoc, hereis, here-string or here-script) is a file literal or input stream literal: it is a section of a source code file that is treated as if it were a separate file.

bashのヒアドキュメントを活用する - Qiita
知ると便利なヒアドキュメント - Qiita
シェルスクリプトでヒアドキュメントを使う様々な方法 (bash / zsh)|Everything you do is practice
bashのメモ ( < << <<< <() ) - なんかやる

17 May 2019追記
awkの組み込み関数もややこしい
機能はman awkで見れる

split(s, a [, r [, seps] ])

Split the string s into the array a and the separators array seps on the regular expression r, and return the number of fields. If r is omitted, FS is used instead. The arrays a and seps are cleared
first. seps[i] is the field separator matched by r between a[i] and a[i+1]. If r is a single space, then leading whitespace in s goes into the extra array element seps[0] and trailing whitespace
goes into the extra array element seps[n], where n is the return value of split(s, a, r, seps). Splitting behaves identically to field splitting, described above.


gensub(r, s, h [, t])

Search the target string t for matches of the regular expression r. If h is a string beginning with g or G, then replace all matches of r with s. Otherwise, h is a number indicating which match of r
to replace. If t is not supplied, use $0 instead. Within the replacement text s, the sequence \n, where n is a digit from 1 to 9, may be used to indicate just the text that matched the n'th parenthesized subexpression. The sequence \0 represents the entire matched text, as does the character &. Unlike sub() and gsub(), the modified string is returned as the result of the function, and the original target string is not changed.

追記ここまで

RNA-seqその8、edgeRを使ったTMM正規化

21 June 2019追記
文章を少し修正した
追記ここまで

featrureCountの結果をedgeRを使って正規化する

このtutorialのpdfファイルと以下のリンク先を参考にした

Statistical analysis of RNA-Seq data
http://www.nathalievialaneix.eu/doc/pdf/tutorial-rnaseq.pdf
ここにPDFのリンクがある
NV² | Biostatistics analysis of RNA-Seq data

実際の処理方法はこちらの方法を
参考にさせてもらった
edgeR - macでインフォマティクス
シロイヌナズナのRNA seq解析 - macでインフォマティクス

edgeRのユーザーズガイドも必読
edgeR: リードカウントから発現変動遺伝子を検出 - Heavy Watal

以前書いた正規化についてのメモ
mecobalamin.hatenablog.com


これまでにやったのは以下の処理

  1. シークエンスデータのダウンロード
  2. trimmomaticでデータのクリーニング
  3. Hisat2でマッピング
  4. featureCountでリードをカウント

mecobalamin.hatenablog.com

NCBIから2種類の細胞種について
3つずつファイルをダウンロードして
それぞれに上記と同様の処理をした
本来ならbiological replicate間で比較をするそうだが、
処理に使ったデータは細胞種が同じなだけで
実験的には無関係のデータであり、
"biological replicateではない"
よくないとは思うけど、スクリプトの確認のため
これらを使って正規化までの処理を確認する

なお6個ファイル処理するのに
手持ちのPCで数日かかった
1TBの外付けHHDでは全く足りず
中間ファイルを消しつつ処理した

まずRを使ってカウントしたファイルの中身を確認する
すでに読み込まれているコマンドや関数を削除して
edgeRのライブラリの読み込む

グラフはggplotを使わずに表示している
見た目の問題なので表示は参照先と同じにできているはず

rm(list = ls(all = T))
library(edgeR)

つぎにリードカウントのファイルを読み込む

# 表示サイズを設定
options(digits = 4, width = 120)
# counts.txtのあるディレクトリのパス
setwd("c:\\path\\to\\countfile")

# リードカウントのファイルを読み込む
Results.FC <- read.table("counts.txt", header = TRUE, skip = 1, sep="\t")
Names.ID <- c("Geneid", "Chr", "Start", "End", "Strand", "Length", "H1", "H2", "H3", "J1", "J2", "J3")
colnames(Results.FC) <- Names.ID
rownames(Results.FC) <- sub("\\.[0-9]{1,}$", "", as.character(Results.FC[, 1]))
Count.Data <- as.matrix(Results.FC[, 7:12])

H1, H2, H3及びJ1, J2, J3はカウントの列につけた任意の名前
HとJで細胞種が異なる
featureCountの出力では読み込んだbamファイルの名前がついているが
わかりやすくするために書き換える

行名はGeneidからピリオド以下を削除して使う
カウントデータは7-12列目でCount.Dataを計算に使う

Results.FCの中身を確認するとこんな感じ

> head(Results.FC)
                     Geneid                                                    Chr
NR_046018       NR_046018.2                                         chr1;chr1;chr1
NR_024540       NR_024540.1 chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1
NR_106918       NR_106918.1                                                   chr1
XR_001737835 XR_001737835.1                                         chr1;chr1;chr1
NR_036051       NR_036051.1                                                   chr1
NR_026818       NR_026818.1                                         chr1;chr1;chr1
                                                                         Start
NR_046018                                                    11874;12613;13221
NR_024540    14362;14970;15796;16607;16858;17233;17606;17915;18268;24738;29321
NR_106918                                                                17369
XR_001737835                                                 29926;30564;30976
NR_036051                                                                30366
NR_026818                                                    34611;35277;35721
                                                                           End
NR_046018                                                    12227;12721;14409
NR_024540    14829;15038;15947;16765;17055;17368;17742;18061;18366;24891;29370
NR_106918                                                                17436
XR_001737835                                                 30039;30667;31295
NR_036051                                                                30503
NR_026818                                                    35174;35481;36081
                            Strand Length H1 H2  H3 J1 J2   J3
NR_046018                    +;+;+   1652  0  0   0  0  0    8
NR_024540    -;-;-;-;-;-;-;-;-;-;-   1769 10 35 199  0 13 1198
NR_106918                        -     68  0  0   6  0  0  149
XR_001737835                 +;+;+    538  0  0   0  0  0    0
NR_036051                        +    138  0  0   0  0  0    1
NR_026818                    -;-;-   1130  0  0   0  0  0    0

正規化をする関数は以下の通り
Statistical analysis of RNA-Seq dataのRaw data filteringのところで
カウントの合計が0以上の行を残しているのでこれを使った
参考にしたサイトでは違うようだ
リードカウントの下限の決め方は色々ありそうなので今後試す

TMMNormalization <- function(Count.Data){
  #Filtering
  #リードカウントが0以上を残す
  Count.Keep <- rowSums(Count.Data) > 0
  Count.Data <- Count.Data[Count.Keep, ]
  Count.Data.log <- log2(Count.Data + 1)
  WriteGraphs(Count.Data, Count.Data.log, "Original")

  Names.Group <- factor(c("positive", "positive", "positive", "negative", "negative", "negative"))

  # edgeRを使った正規化係数の計算
  Count.DGE <- DGEList(counts = Count.Data, group = Names.Group)
  Count.DGE <- calcNormFactors(Count.DGE)
  Count.DGE <- estimateCommonDisp(Count.DGE)
  Count.DGE <- estimateTagwiseDisp(Count.DGE)
  Results.DGE <- exactTest(Count.DGE)

  Table.DGE <- as.data.frame(topTags(Results.DGE, n = nrow(Count.Data)))
  Tag.DGE <- as.logical(Table.DGE$FDR < 0.05)
  Names.DGE <- rownames(Table.DGE)[Tag.DGE]

  # 正規化係数とライブラリサイズの取得
  Norm.Factors <- Count.DGE$samples$norm.factors
  Lib.Size <- Count.DGE$samples$lib.size

  # 正規化と対数の計算
  Count.Norm <- t(apply(Count.Data, 1, function(x) x / (Norm.Factors * Lib.Size) * 10^6))
  Count.Norm.log <- log2(Count.Norm + 1)
  WriteGraphs(Count.Norm, Count.Norm.log, "TMMNorm")
  return(Count.Norm.log)
}

正規化前のカウント

>print(head(Count.Data.log))
                H1   H2    H3 J1    J2     J3
NR_046018    0.000 0.00 0.000  0 0.000  3.170
NR_024540    3.459 5.17 7.644  0 3.807 10.228
NR_106918    0.000 0.00 2.807  0 0.000  7.229
NR_036051    0.000 0.00 0.000  0 0.000  1.000
XR_001737584 0.000 0.00 0.000  0 0.000  3.585
XR_001737581 0.000 0.00 0.000  0 0.000  3.585

カウントの合計が0の行が除かれている

で、正規化されるとこうなる

>print(head(Count.Norm.log))
                 H1    H2      H3 J1     J2       J3
NR_046018    0.0000 0.000 0.00000  0 0.0000 0.011020
NR_024540    0.1552 1.109 0.94713  0 0.5704 1.103176
NR_106918    0.0000 0.000 0.03981  0 0.0000 0.192592
NR_036051    0.0000 0.000 0.00000  0 0.0000 0.001382
XR_001737584 0.0000 0.000 0.00000  0 0.0000 0.015131
XR_001737581 0.0000 0.000 0.00000  0 0.0000 0.015131


DGEListでedgeRで使う形式にデータを変えている
正規化係数やライブラリのサイズも出力されている

Count.DataのH1のカウントについてヒストグラム出力

f:id:mecobalamin:20190510125558p:plain
H1のヒストグラム
今度は底を2で対数を計算して出力
f:id:mecobalamin:20190510125435p:plain
H1のヒストグラムを対数で表示
Count.Data.logを表示している

正規化前の分布をbox plotで表示する

f:id:mecobalamin:20190510125712p:plain
正規化前の分布
結構ばらついている様に見える
元のデータがバラバラだからしょうがないか

ここはまだ理解が怪しいところなのだが、
TMM正規化にはMA plotで閾値を決めて
残ったリードカウントを使っているっぽい

H1とH2のリードカウントをMA plot表示するとなんとなくそれっぽい

f:id:mecobalamin:20190512004945p:plain
H1-H2のMA-plot
J1-J2の組み合わせだと集団が別れてる
外れ値?
除く必要がありそうだけどやり方がわからないので
とりあえずそのまま処理を続ける
f:id:mecobalamin:20190512005039p:plain
J1-J2のMA-plot
tutorialのpdfに出てたコード使って全部の組み合わせで出力するようにした
全部組み合わせる必要はなさそうだけどコードの関係。

MAPlot <- function(d){
  c <- colnames(d)
  for(i in c){
    for(j in c){
      if(i != j){
        M = d[,i] - d[,j]
        A = (d[,i] + d[,j]) / 2
        df <- data.frame(A, M)
        print(head(df))
        png(sprintf("MAplot_%s-%s.png", i, j), width = 600, height = 450)
        plot(df, xlab = i, ylab = j, main = sprintf("%s-%s", i, j))
        dev.off()
      }
    }
  }
}

edgeRの関数を使ってカウント全部をプロットすると
なんか違う形の分布になった

f:id:mecobalamin:20190510125901p:plain
plotSmearで作ったMA plot
データの下処理の問題か?

plotSmear(Results.DGE, de.tags. = Names.DGE, cex = 0.3)

リードカウントをライブラリサイズと正規化係数で正規化し、box plot表示する

f:id:mecobalamin:20190510130011p:plain
TMM正規化後の分布
分布が大体揃っているように見える
外れ値が増えたようだがどうなんだろう

ヒストグラムとbox plotの出力には以下の関数を使っている

WriteGraphs <- function(d.norm, d.log, g.tag){

  h.tag <- "H1"

  # グラフの保存
  # ヒストグラムの表示
  png(sprintf("histogram_count_norm_%s_%s.png", g.tag, h.tag), width = 600, height = 450)
  hist(d.norm[, h.tag], xlab = h.tag, ylab = "count", main = "norm")
  dev.off()
  #対数
  png(sprintf("histogram_count_log_%s_%s.png", g.tag, h.tag), width = 600, height = 450)
  hist(d.log[, h.tag], xlab = h.tag, ylab = expression(log[2](count + 1)), main = "log")
  dev.off()

  #box plot
  #対数
  png(sprintf("boxplot_count_log_%s.png", g.tag), width = 600, height = 450)
  boxplot.matrix(d.log, ylab = expression(log[2](count + 1)))
  dev.off()
}

試しにEGFRの数を比較する
正規化したデータから目的のAccession Numberを探して
カウント数を調べる

FindGeneInterest <- function(d){
  #find genes of interest
  acc.egfr <- c("NM_005228", "NM_201282", "NM_201283", "NM_201284", "NM_001346897", "NM_001346898", "NM_001346899", "NM_001346900", "NM_001346941")
  iso.egfr <- c("a", "b", "c", "d", "e", "f", "g", "h", "i")

  var.egfr <- t(rbind(acc.egfr, iso.egfr))

  egfr <- c()
  for(i in var.egfr[,1]){
    egfr <- rbind(egfr, d[i,])
  }
  print(data.frame(var.egfr, egfr))
}

これが出力

      acc.egfr iso.egfr    H1    H2    H3 J1 J2 J3
1    NM_005228        a 2.707 5.064 5.414  0  0  0
2    NM_201282        b 1.247 3.233 4.108  0  0  0
3    NM_201283        c 0.710 2.865 3.213  0  0  0
4    NM_201284        d 1.261 3.238 4.120  0  0  0
5 NM_001346897        e 1.811 4.116 5.002  0  0  0
6 NM_001346898        f 1.830 4.116 5.006  0  0  0
7 NM_001346899        g 2.697 5.064 5.412  0  0  0
8 NM_001346900        h 2.684 5.038 5.403  0  0  0
9 NM_001346941        i 2.638 5.012 5.252  0  0  0

前後するがこっちが正規化前の対数でないリードカウント

      acc.egfr iso.egfr  H1  H2   H3 J1 J2 J3
1    NM_005228        a 487 981 8931  0  0  0
2    NM_201282        b 121 254 3482  0  0  0
3    NM_201283        c  56 190 1774  0  0  0
4    NM_201284        d 123 255 3514  0  0  0
5 NM_001346897        e 221 494 6658  0  0  0
6 NM_001346898        f 225 494 6675  0  0  0
7 NM_001346899        g 483 981 8914  0  0  0
8 NM_001346900        h 478 963 8861  0  0  0
9 NM_001346941        i 460 945 7957  0  0  0

HとJはともに培養細胞のシークエンス結果で、
HにはEGFRの発現があり、Jには発現していないと言われている
実際その通りになっているように見える
H1とH3では桁が変わるほどカウント数に違いがあったが
正規化により桁はだいたい揃ってる
似たような発現量ってことでいいのだろうか

もっともMA plotではシークエンスの結果が不均一っぽい
正規化が正しく計算できていないと思う
なのでこのまま数字の比較はできないかも
そもそもreplicateでないシークエンスのデータを
無理やりまとめて正規化してるからからか
うまく外れているところを除くか使うデータを変えてみるか。。。

処理の流れとしてはこれで良さそうではあるが。。。

コードをまとめる

rm(list = ls(all = T))
library(edgeR)

WriteGraphs <- function(d.norm, d.log, g.tag){

  h.tag <- "H1"

  # グラフの保存
  # ヒストグラムの表示
  png(sprintf("histogram_count_norm_%s_%s.png", g.tag, h.tag), width = 600, height = 450)
  hist(d.norm[, h.tag], xlab = h.tag, ylab = "count", main = "norm")
  dev.off()
  #対数
  png(sprintf("histogram_count_log_%s_%s.png", g.tag, h.tag), width = 600, height = 450)
  hist(d.log[, h.tag], xlab = h.tag, ylab = expression(log[2](count + 1)), main = "log")
  dev.off()

  #box plot
  #対数
  png(sprintf("boxplot_count_log_%s.png", g.tag), width = 600, height = 450)
  boxplot.matrix(d.log, ylab = expression(log[2](count + 1)))
  dev.off()
}

MAPlot <- function(d){
  c <- colnames(d)
  for(i in c){
    for(j in c){
      if(i != j){
        M = d[,i] - d[,j]
        A = (d[,i] + d[,j]) / 2
        df <- data.frame(A, M)
        print(head(df))
        png(sprintf("MAplot_%s-%s.png", i, j), width = 600, height = 450)
        plot(df, xlab = i, ylab = j, main = sprintf("%s-%s", i, j))
        dev.off()
      }
    }
  }
}

FindGeneInterest <- function(d){
  # find genes of interest
  acc.egfr <- c("NM_005228", "NM_201282", "NM_201283", "NM_201284", "NM_001346897", "NM_001346898", "NM_001346899", "NM_001346900", "NM_001346941")
  iso.egfr <- c("a", "b", "c", "d", "e", "f", "g", "h", "i")

  var.egfr <- t(rbind(acc.egfr, iso.egfr))

  egfr <- c()
  for(i in var.egfr[,1]){
    egfr <- rbind(egfr, d[i,])
  }
  print(data.frame(var.egfr, egfr))
}

TMMNormalization <- function(Count.Data){
  # Filtering
  # リードカウントが0以上を残す
  Count.Keep <- rowSums(Count.Data) > 0
  Count.Data <- Count.Data[Count.Keep, ]
  Count.Data.log <- log2(Count.Data + 1)
  #print(head(Count.Data.log))

  # グラフの出力
  MAPlot(Count.Data.log)
  WriteGraphs(Count.Data, Count.Data.log, "Original")

  Names.Group <- factor(c("positive", "positive", "positive", "negative", "negative", "negative"))

  # edgeRを使った正規化係数の計算
  Count.DGE <- DGEList(counts = Count.Data, group = Names.Group)
  Count.DGE <- calcNormFactors(Count.DGE)
  Count.DGE <- estimateCommonDisp(Count.DGE)
  Count.DGE <- estimateTagwiseDisp(Count.DGE)
  Results.DGE <- exactTest(Count.DGE)

  Table.DGE <- as.data.frame(topTags(Results.DGE, n = nrow(Count.Data)))
  Tag.DGE <- as.logical(Table.DGE$FDR < 0.05)
  Names.DGE <- rownames(Table.DGE)[Tag.DGE]

  # plotSemerを使ってMA plotの出力
  png("MAplot_plotSemer_DGEresults.png", width = 600, height = 450)
  plotSmear(Results.DGE, de.tags. = Names.DGE, cex = 0.3)
  dev.off()

  # 正規化係数とライブラリサイズの取得
  Norm.Factors <- Count.DGE$samples$norm.factors
  Lib.Size <- Count.DGE$samples$lib.size

  # 正規化と対数の計算
  Count.Norm <- t(apply(Count.Data, 1, function(x) x / (Norm.Factors * Lib.Size) * 10^6))
  Count.Norm.log <- log2(Count.Norm + 1)
  #print(head(Count.Norm.log))
  # グラフの出力
  WriteGraphs(Count.Norm, Count.Norm.log, "TMMNorm")
  return(Count.Norm.log)
}

options(digits = 4, width = 120)
setwd("c:\\path\\to\\Count")

Results.FC <- read.table("counts.txt", header = TRUE, skip = 1, sep="\t")
Names.ID <- c("Geneid", "Chr", "Start", "End", "Strand", "Length", "H1", "H2", "H3", "J1", "J2", "J3")
colnames(Results.FC) <- Names.ID
rownames(Results.FC) <- sub("\\.[0-9]{1,}$", "", as.character(Results.FC[, 1]))
Count.Data <- as.matrix(Results.FC[, 7:12])

Data.Norm <- TMMNormalization(Count.Data)
FindGeneInterest(Data.Norm)

一部Atomスクリプトを実行するために
必要な書き方になっている
ここのprint()など

print(head(d))

はRGuiではprintは必要ない

RNA-seqその7、TrinityでアセンブルしてBlastにかける

TrinityでシークエンスデータをアセンブルしてBlastにかける

そのためにまずhisat2でゲノムにマッピングしたデータをアセンブルし、
できたfastaファイルをもとにblastのデータベースを作る
そしてこのデータベースと目的の遺伝子を
記入したfastaファイルをつかってblastを行う

Home · trinityrnaseq/trinityrnaseq Wiki · GitHub
Trinity | de novo アセンブリー


まず準備

Trinityで使うworkdirにRAMディスクを指定する
RAMディスクはImDiskで作成する
作成方法は以前の記事を参照
mecobalamin.hatenablog.com

ディレクトリは

/mnt/r/tmp/

を使い、サイズは10GBで作成

Trinityにかけるファイルは複数あって
list_sra.txtにダウンロードしてきたsraファイルの
accession numberが記入してある


Trinityを実行するシェルスクリプトはこれ
genome_guided_bamのオプションでhisat2で作ったbamファイルを指定する
メモリ10Gは少ないがRAMディスクと
windowsに必要な分の残りなのでしょうがない

#!/bin/bash

echo "trinity began"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_tmp="/mnt/r/tmp"

echo "open sra list"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

for i in ${list_filename}
  do
    echo "analyze ${i}"
    file_in="sort_${i}.bam"
    dir_in="hisat_${i}"
    dir_out="trinity_${i}"

    echo "makind directories"
    mkdir -p ${dir_tmp}/${dir_out}
    mkdir -p ${dir_working}/${dir_out}

    Trinity \
      --genome_guided_bam ${dir_working}/${dir_in}/${file_in} \
      --genome_guided_max_intron 10000 \
      --verbose \
      --output ${dir_working}/${dir_out} \
      --workdir ${dir_tmp}/${dir_out} \
      --CPU 4 \
      --max_memory 10G

    mv ${dir_working}/${dir_out}/Trinity-GG.fasta ${dir_working}/${dir_out}/Trinity-GG_${i}.fasta

    TrinityStats.pl ${dir_working}/${dir_out}/Trinity-GG_${i}.fasta > \
      ${dir_working}/${dir_out}/Trinity-GG_${i}.stats

    rm -r ${dir_tmp}/${dir_out}
  done

結果はtrinity_[accession number]というディレクトリに保存される

計算に時間がとてもかかった
試しに以前hisatでマッピングしたbamファイル2つについてアセンブルしたところ
計算が終わるまでに一つは20時間、もう一つは12時間ぐらいかかった
その間HDDはアクセスランプが点きっぱなし

Trinityのhelpを見ると

# --workdir :where Trinity phase-2 assembly computation takes place (defaults to --output setting).
# (can set this to a node-local drive or RAM disk)

とのことなのでそれ以外の処理のときはHDDを使ってる様子
RAMディスク10Gも使っていない様子だし、最適化が必要かも

TrinityStats.plはTrinityの結果に対して統計結果を出力してくれる

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  17930
Total trinity transcripts:      17934
Percent GC: 52.83

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2017
        Contig N20: 1415
        Contig N30: 1016
        Contig N40: 767
        Contig N50: 587

        Median contig length: 335
        Average contig: 497.02
        Total assembled bases: 8913645


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2017
        Contig N20: 1415
        Contig N30: 1016
        Contig N40: 767
        Contig N50: 587

        Median contig length: 335
        Average contig: 497.07
        Total assembled bases: 8912460


次に、Trinityで作ったfastaファイルからblast用のデータベースの作成する

#!/bin/bash

echo "making blast database"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_tmp="/mnt/r/tmp"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

db_type="nucl"

for i in ${list_filename}
  do
    mkdir -p ${dir_tmp}

    dir_db="trinity_${i}"
    #fasta_db="trinity_${i}.Trinity.fasta"
    fasta_db="Trinity-GG_${i}.fasta"

    cp ${dir_working}/${dir_db}/${fasta_db} ${dir_tmp}

    makeblastdb \
      -in ${dir_tmp}/${fasta_db} \
      -dbtype ${db_type}

    mv ${dir_tmp}/${fasta_db}.* ${dir_working}/${dir_db}/
    rm -r ${dir_tmp}
  done

このデータベースに対してblastをかける

塩基配列同士なのでblastnを使う
試しにEGFRをqueryにしてかけてみる
EGFRのaccession numberはNM_005228で
EGFR.faという名前のfastaファイルをFasta_ErbBsというディレクトリに保存済み

#!/bin/bash

echo "blast began"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_query="${dir_working}/Fasta_ErbBs"

dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//')

fasta_query="EGFR.fa"
pro_blast="blastn"

for i in ${list_filename}
  do
    echo "analyze ${i}"
    dir_db="trinity_${i}"
    fasta_db="Trinity-GG_${i}.fasta"
    dir_blast="blast_${i}"

    mkdir -p ${dir_working}/${dir_blast}

    ${pro_blast} \
      -query ${dir_query}/${fasta_query} \
      -db ${dir_working}/${dir_db}/${fasta_db} \
      -evalue 1e-30 \
      -outfmt 6 \
      -num_threads 4 \
      -num_alignments 10 \
      > ${dir_working}/${dir_blast}/${pro_blast}_${fasta_query/\.fa/}_${fasta_db/\.fasta/}.csv
  done
echo "blast ended"

結果はcsvファイルに保存される

NM_005228.5     TRINITY_GG_6231_c1_g1_i1        100.000 408     0       0       8563    8970    1       408     0.0     754
NM_005228.5     TRINITY_GG_6231_c0_g1_i1        99.541  218     1       0       9274    9491    218     1       1.17e-109       398
NM_005228.5     TRINITY_GG_6231_c2_g1_i1        100.000 202     0       0       3596    3797    1       202     1.96e-102       374

試しにcsvファイルからTRINITY_GG_6231_c1_g1_i1を検索かけて
一致したのがこの配列
lessのコマンドで検索できる

>TRINITY_GG_6231_c1_g1_i1 len=408 path=[0:0-407]
GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAGCCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGGAGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGGCTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAAGTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGCAGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCTGAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT

この配列をNCBIのBlastにかけるとクエリーに使った遺伝子と一致してるのでまあいいのかな
こんな確認でいいのかわからないけど

Homo sapiens epidermal growth factor receptor (EGFR), transcript variant 1, mRNA
Sequence ID: NM_005228.5Length: 9905Number of Matches: 1
Related Information
Gene-associated gene details
PubChem BioAssay-bioactivity screening
Genome Data Viewer-aligned genomic context
Range 1: 8563 to 8970GenBankGraphicsNext MatchPrevious Match
Alignment statistics for match #1
Score	Expect	Identities	Gaps	Strand
754 bits(408)	0.0	408/408(100%)	0/408(0%)	Plus/Plus
Query  1     GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAG  60
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8563  GTCTGGAGAGCCTAATAATGTTCAGCACACTTTGGTTAGTTCACCAACAGTCTTACCAAG  8622

Query  61    CCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGG  120
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8623  CCTGGGCCCAGCCACCCTAGAGAAGTTATTCAGCCCTGGCTGCAGTGACATCACCTGAGG  8682

Query  121   AGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGG  180
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8683  AGCTTTTAAAAGCTTGAAGCCCAGCTACACCTCAGACCGATTAAACGCAAATCTCTGGGG  8742

Query  181   CTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAA  240
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8743  CTGAAACCCAAGCATTCGTAGTTTTTAAAGCTCCTGAGGTCATTCCAATGTGCGGCCAAA  8802

Query  241   GTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGC  300
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8803  GTTGAGAACTACTGGCCTAGGGATTAGCCACAAGGACATGGACTTGGAGGCAAATTCTGC  8862

Query  301   AGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCT  360
             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8863  AGGTGTATGTGATTCTCAGGCCTAGAGAGCTAAGACACAAAGACCTCCACATCTGTCGCT  8922

Query  361   GAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT  408
             ||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  8923  GAGAGTCAAGAACCTGAACAGAGTTTCCATGAAGGTTCTCCAAGCACT  8970


スクリプトで、出力ファイル指定で変数にスラッシュが付いている

${fasta_query/\.fa/}

これは変数に対して削除・置換をしてくれる
この場合はfasta_queryの中身、EGFR.faから.faを削除する
後半部分も同様の処理をしている