mecobalamin’s diary

人間万事塞翁が馬

https://help.hatenablog.com/entry/developer-option

フォーマット変換、GTF -> BED -> BAM

GTFファイルからBAMファイルにフォーマットを変換したい

まずフォーマットについて。
GTF/GFF
GTFとGFFフォーマット - macでインフォマティクス

SAM/BAM
BAM – 遺伝子発現解析(マイクロアレイ解析, RNA-seq)
https://bi.biopapyrus.jp/rnaseq/mapping/sam.html

BED
BEDフォーマット完全解説 - アメリエフの技術ブログ

あとgoogleの検索で見つかったpdfファイルが詳しかった
"NGS, 基本データフォーマット"で検索かけると
基礎生物学研究所での講義ファイルが見つかる
PDFファイルなのでリンクは貼らない

次にコマンドについて

BEDOPSはGTF/GFFをBEDに変換できる
BEDOPS: the fast, highly scalable and easily-parallelizable genome analysis toolkit — BEDOPS v2.4.41
VCF, GTF, GFF などを BED に変換 する BEDOPS - macでインフォマティクス

インストールはlinux用バイナリをダウンロードして展開
パスの通っているディレクトリにリンクを貼る
使っている環境ではバイナリを~/usr/localにおいて
リンクを~/binに張っている

convert2bedとsort-bedはgtf2bedの実行に必要

$ ln -s $HOME/usr/local/bedops/bin/gff2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/gtf2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/convert2bed $HOME/bin
$ ln -s $HOME/usr/local/bedops/bin/sort-bed $HOME/bin

BEDからBAMへはbedtoolsでできる
https://cell-innovation.nig.ac.jp/surfers/BEDtools.html
https://bi.biopapyrus.jp/rnaseq/mapping/format-convert.html
ここからダウンロード
Installation — bedtools 2.30.0 documentation
展開してパスの通ったディレクトリにリンクを貼る


準備できたのでさっそく使ってみる

まずgtf2bedを使ってhg38のgtfをbedファイルに変換する
次にbedtoossを使ってbedをbamに変換する
bedtoolsのディレクトリにあるhuman.hg38.genomeをgenomeサイズとして指定する

$ gtf2bed < rename_UCSC.hg38.gtf > rename_UCSC.hg38.bed
$ bedtools bedtobam -i rename_UCSC.hg38.bed -g $HOME/usr/local/bedtools2/genomes/human.hg38.genome > rename_UCSC.hg38.bam

bashのリダイレクトの使い方未だによくわかんないんだけど
gtfファイルをgtf2bedに渡してbedファイルに出力しているらしい

変換前のgtfファイルはこうで

chr1    hg38_ncbiRefSeq exon    11874   12227   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    12613   12721   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    13221   14409   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    14362   14829   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";
chr1    hg38_ncbiRefSeq exon    14970   15038   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";
chr1    hg38_ncbiRefSeq exon    15796   15947   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";

変換後のBAMファイルはこう

NR_046018.2     0       chr1    11874   255     354M    *       0       0       *       *
NR_046018.2     0       chr1    12613   255     109M    *       0       0       *       *
NR_046018.2     0       chr1    13221   255     1189M   *       0       0       *       *
NR_024540.1     16      chr1    14362   255     468M    *       0       0       *       *
NR_024540.1     16      chr1    14970   255     69M     *       0       0       *       *
NR_024540.1     16      chr1    15796   255     152M    *       0       0       *       *

うまく変換できたかの判断はまた後ほど

23 April 2019追記
盛大な勘違いでやらなくていい処理だった
消そうかとも思ったけどそのうち見直すかもなので。。。

ドライブの速度比較

ファイルの書き込みに時間がかかってる
実際にディスクの書き込み速度がどの程度違うか比較した

比較したのはPC内蔵のSSD、HDD、USB3.1接続のHDD、USBメモリ、それとRAMディスク

計測はCrystal Disk Markを使った
USBメモリのみドライブとして認識されていなかったのでフォルダを指定して計測した

f:id:mecobalamin:20190419105912p:plain
SSDSATA接続
f:id:mecobalamin:20190419105958p:plain
HDD、SATA接続
f:id:mecobalamin:20190419110101p:plain
HDD、USB3.1 Type C接続
f:id:mecobalamin:20190419110212p:plain
USBメモリUSB3.1接続
f:id:mecobalamin:20190419110255p:plain
RAMディスク

SSDにデータを入れて処理するのがベストだろうけど
システムも入っているしサイズが足りないのと壊れたときの被害が大きい
HDDを使うなら内蔵と外付けで読み書きの速度はそこまで変わらない
壊れて後で困るぐらいなら外付けドライブを使うのが良さそう
テンポラリディスクはRAMディスクを使うのがいいが、
メモリのサイズが小さいので使い所を選ぶ

とりあえずデータの保存に外付けHDDを使って
テンポラリファイルの処理にRAMディスクを使用するのが
今の環境では速度が出せて万が一のとき被害が少ないと思う

そもそもWSLでI/Oが遅いのはしょうがないっぽい
WSLのI/Oが遅いことに対するMS技術者からの回答 | OPCDiary
Microsoft、WSLのI/Oパフォーマンス改善に取り組み | マイナビニュース
WSL vs VM for 開発作業 - Qiita


14 May 2019追記
bashのddコマンドで書き込みのベンチマークを測定
hdparmコマンドとddコマンドでハードディスクの簡易ベンチマークをする | 怪しい物を開発するブログ

d = 内蔵ssd
e = 内蔵hdd
h = 外付けhdd、USB3.1接続
usb = USBメモリUSB3.1接続
r = RAM disk

$ dd if=/dev/zero of=/mnt/d/tmp/testimg.tmp bs=1M count=1024
1024+0 レコード入力
1024+0 レコード出力
1073741824 bytes (1.1 GB, 1.0 GiB) copied, 0.791183 s, 1.4 GB/s

$ dd if=/dev/zero of=/mnt/e/tmp/testimg.tmp bs=1M count=1024
1024+0 レコード入力
1024+0 レコード出力
1073741824 bytes (1.1 GB, 1.0 GiB) copied, 2.54097 s, 423 MB/s

$ dd if=/dev/zero of=/mnt/h/tmp/testimg.tmp bs=1M count=1024
1024+0 レコード入力
1024+0 レコード出力
1073741824 bytes (1.1 GB, 1.0 GiB) copied, 15.5029 s, 69.3 MB/s

$ dd if=/dev/zero of=/mnt/usb/tmp/testimg.tmp bs=1M count=1024
1024+0 レコード入力
1024+0 レコード出力
1073741824 bytes (1.1 GB, 1.0 GiB) copied, 50.506 s, 21.3 MB/s

$ dd if=/dev/zero of=/mnt/r/tmp/testimg.tmp bs=1M count=1024
1024+0 レコード入力
1024+0 レコード出力
1073741824 bytes (1.1 GB, 1.0 GiB) copied, 0.796872 s, 1.3 GB/s

書き込み速度(GB/s)は
SSD = RAM disk > HDD内蔵 > HDD外付け > USBメモリ
CrystalDiskMarkと数値の比較はできないんだろうけど
バイス間の書き込み速度比が結構違う
bash上ではSSDと内蔵HDDとでそんなに差が出てない

なおUSBデバイスは手動でマウントする
Windows 10の「WSL」でネットワークドライブなどをマウントする:Tech TIPS - @IT

$ sudo mkdir /mnt/usb
$ sudo mount -t drvfs F: /mnt/usb

アンマウントも手動

$ sudo umount /mnt/usb

windows10 homeのwslにdockerをインストールする

5 October 2020 追記
誤字の修正
追記ここまで

1 June 2019追記
インストールをやり直したらうまくいかない操作があったので書き直した
追記ここまで

使いたいコマンドがDockerで配布されてたのでインストールしたい

Dockerとはなんぞや。。。
無償の「Docker for Windows」で手軽にLinuxコンテナを利用する (1/2):Windows管理者のためのDocker入門 - @IT
Windows Server 1709上でLinuxコンテナを動かす(lcow) - Qiita

残念ながら使用中のPCはwindows10 home。
Windows10 homeではhyper-Vが使えないが、
Docker ToolBox + VirtualBoxの組み合わせで
dockerを使えるようにしている人達がいる

windows 10 home で docker を導入するメモ - Qiita
Windows 10 Home のための、Docker Toolbox をインストールして WSL から使う方法 | ラボラジアン
Docker Toolboxのインストール:Windows編 - Qiita
WSLでDockerを使う - Qiita

ただしDocker ToolBoxはlegacy solutionなので
ある意味イレギュラーな使い方かも
そのうち使えなくなることも覚悟しておく

Legacy desktop solution. Docker Toolbox is for older Mac and Windows systems that do not meet the requirements of Docker Desktop for Mac and Docker Desktop for Windows. We recommend updating to the newer applications, if possible.

インストールの順序は

  1. Windows側にDocker Tool boxをインストールする
  2. WSL側にDockerをインストールする
  3. .profileを編集する
  4. 動作の確認

1. まずDocker Tool boxをインストールする
ダウンロードはここから。
v18.09.3が最新
Releases · docker/toolbox · GitHub
windowsはexeファイルでいいはず。。。

インストール先以外はデフォルトのままインストールした
Install Docker Toolbox on Windows | Docker Documentation

Step 2: Install Docker Toolbox
In this section, you install the Docker Toolbox software and several “helper” applications. The installation adds the following software to your machine:

Docker Client for Windows
Docker Toolbox management tool and ISO
Oracle VM VirtualBox
Git MSYS-git UNIX tools
If you have a previous version of VirtualBox installed, do not reinstall it with the Docker Toolbox installer. When prompted, uncheck it.

Docker Quickstart Terminalのショートカットをダブルクリックすると
ターミナルが起動する。インストール作業が始まる。
時間が掛かるがクジラが表示されたあとに
コマンドプロンプトが表示される

                        ##         .
                  ## ## ##        ==
               ## ## ## ## ##    ===
           /"""""""""""""""""\___/ ===
      ~~~ {~~ ~~~~ ~~~ ~~~~ ~~~ ~ /  ===- ~~~
           \______ o           __/
             \    \         __/
              \____\_______/

docker is configured to use the default machine with IP 192.168.99.100
For help getting started, check out the docs at https://docs.docker.com

Start interactive shell

hogehoge@fugahuga MINGW64 /d/Docker Toolbox
$

ちなみにhogehogeはpcにログインしているユーザー名
fugafugaはpc名が表示されている

Docker Quickstart Terminalでバージョン確認のコマンドを実行
clientとserverの両方が表示されればまずはオッケー

$ docker version
Client:
 Version:           18.09.3
 API version:       1.39
 Go version:        go1.12
 Git commit:        774a1f4eee
 Built:             Mon Mar  4 10:36:44 2019
 OS/Arch:           windows/amd64
 Experimental:      false

Server: Docker Engine - Community
 Engine:
  Version:          18.09.3
  API version:      1.39 (minimum version 1.12)
  Go version:       go1.10.8
  Git commit:       774a1f4
  Built:            Thu Feb 28 06:40:51 2019
  OS/Arch:          linux/amd64
  Experimental:     false


2. 次にwsl側からdockerのコマンドをインストールする

ここと同じ手順
WSL の docker client から、Docker for Windows の docker daemon を使う手順 | ラボラジアン
以下の順でコマンドを実行する

途中でgitからダウンロードするdocker-composeのバージョンはここに合わせた
Install Docker Compose | Docker Documentation

$ sudo apt update
$ sudo apt install \
     apt-transport-https \
     ca-certificates \
     curl \
     software-properties-common
$ curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo apt-key add -
$ sudo apt-key fingerprint 0EBFCD88
$ sudo add-apt-repository    "deb [arch=amd64] https://download.docker.com/linux/ubuntu\
    $(lsb_release -cs) \
    stable"
$ sudo apt update
$ sudo apt install docker-ce
$ sudo curl -L "https://github.com/docker/compose/releases/download/1.23.2/docker-compose-$(uname -s)-$(uname -m)" -o /usr/local/bin/docker-compose
$ sudo chmod a+rx /usr/local/bin/docker-compose
$ docker-compose --version
$ sudo curl -L https://raw.githubusercontent.com/docker/compose/$(docker-compose version --short)/contrib/completion/bash/docker-compose -o /etc/bash_completion.d/docker-compose
$ sudo chmod a+r /etc/bash_completion.d/docker-compose
$ source /etc/bash_completion.d/docker-compose

3. 次に$HOME/.profileの編集する
こちらを参考にした
Windows 10 Home のための、Docker Toolbox をインストールして WSL から使う方法 | ラボラジアン

まずはDocker Toolboxで

docker-machine.exe confing

を実行し、仮想マシンIPアドレスとポート番号を取得する

次に$HOME/.profileに仮想マシンIPアドレスとポート番号、及びをcertsのpathを追記した

export DOCKER_HOST="tcp://192.168.99.100:2376"
export DOCKER_CERT_PATH=/mnt/c/Users/hogehoge/.docker/machine/machines/default
export DOCKER_TLS_VERIFY=1

hogehogeは自分のユーザー名に書き換える
/mnt/c/Users/hogehoge/.docker/machine/certsにも秘密鍵?のファイルがあるがこっちではないっぽい
IPアドレスはDocker Quickstart Terminalでdocker-machine configを
実行して表示されているIPアドレスとポート番号を記入した
ググって探したサイトによってはポート番号が2376の場合と
2375の場合とがあって何なのかわからなかったが、
Docker Quickstart Terminalから取得できる
IPアドレスとポート番号が上記の値で、
こちらの環境ではこれで動くようになった

26 September 2020 追記
DOCKER_HOSTの確認は

$ docker-machine.exe env

で行う
実行結果は以下の通り

You can further specify your shell with either 'cmd' or 'powershell' with the --shell flag.
SET DOCKER_TLS_VERIFY=1
SET DOCKER_HOST=tcp://192.168.99.100:2376
SET DOCKER_CERT_PATH=C:\Users\hogehogehugahuga\.docker\machine\machines\default
SEZ DOCKER_MACHINE_NAME=default
SET COMPOSE_CONVERT_WINDOWS_PATHS=true
REM Run this command to configure your shell:
REM @FOR /f "tokens=*" %i IN ('docker-machine.exe env') DO @%i 

DOCKER_HOSTの値を$HOME/.profileに記入する

追記ここまで

追記したら.profileを読み込ませるためにターミナルを再起動する
再起動せずに

$ source .profile

でも動いた

インストールはここまで

4. 最後に動作の確認をする
再起動後wslのターミナルで

$ docker-machine.exe active

を実行すると動いている仮想マシンを確認できる
docker-machine.exeはwslのコマンドではないので拡張子が必要

wslのターミナルでバージョン確認のコマンドを実行する

$ docker version
Client:
 Version:           18.09.4
 API version:       1.39
 Go version:        go1.10.8
 Git commit:        d14af54266
 Built:             Wed Mar 27 18:35:44 2019
 OS/Arch:           linux/amd64
 Experimental:      false

Server: Docker Engine - Community
 Engine:
  Version:          18.09.3
  API version:      1.39 (minimum version 1.12)
  Go version:       go1.10.8
  Git commit:       774a1f4
  Built:            Thu Feb 28 06:40:51 2019
  OS/Arch:          linux/amd64
  Experimental:     false

Docker Quickstart Terminalのときと同じかな

次にコンテナが動作するか確認する

$ docker run hello-world

実行後次のように表示されたらオッケーらしい

Hello from Docker!
This message shows that your installation appears to be working correctly.

To generate this message, Docker took the following steps:
 1. The Docker client contacted the Docker daemon.
 2. The Docker daemon pulled the "hello-world" image from the Docker Hub.
    (amd64)
 3. The Docker daemon created a new container from that image which runs the
    executable that produces the output you are currently reading.
 4. The Docker daemon streamed that output to the Docker client, which sent it
    to your terminal.

To try something more ambitious, you can run an Ubuntu container with:
 $ docker run -it ubuntu bash

Share images, automate workflows, and more with a free Docker ID:
 https://hub.docker.com/

For more examples and ideas, visit:
 https://docs.docker.com/get-started/

仮想マシンの開始・終了もwsl側からできる
仮想マシンの終了は

$ docker-machine.exe stop

起動は次のようにする

$ docker-machine.exe start

動いてるコンテナや仮想マシンがないか確認する
終了についてはこちらを参考にした
windows 10 home で docker を導入するメモ - Qiita

23 April 2020 追記
Oracle VirtualBoxVMのアップデートをしたらdocker runでコンテナを起動できなくなった
docker-machine.exe configで確認すると以下のメッセージ

Error running connection boilerplate: Error checking and/or regenerating the certs: There was an error validating certificates for host "192.168.99.101:2376": x509: certificate is valid for 192.168.99.100, not 192.168.99.101

このサイトを参考に修正できた
qiita.com

追記ここまで

Rで折れ線グラフ

Rで時系列データの折れ線グラフを作る
グラフを書くことよりも日付の処理が面倒だったのでそのメモ

12 December 2020 追記
データは異なるがPythonでも書いてみた
mecobalamin.hatenablog.com
追記ここまで


作りたいのは気温と海水温の時系列グラフ

気温と海水温の軸を左右に分けて、日付の軸は共通にする
折れ線グラフの赤が海水温で青が気温を示す

気温と海水温の時系列データは気象庁からダウンロードする
気温:気象庁、過去の気象データ
気象庁|過去の気象データ・ダウンロード
海水温:東京管区気象台
https://www.jma-net.go.jp/tokyo/sub_index/kaiyou/sst/index.html

ファイルを見ると日付のフォーマットが違うことがわかる
今回に限らず、データによって日付の書き方が違うことはよくある
まず日付のフォーマットを揃える事から始める

気温のデータはcsv形式でダウンロードできて、
いろいろ情報が入っているが、
日付と気温以外のデータをexcelなどで削って使う
削ったデータをRで取り込んでみるとこんな感じ

data.air <- read.table("airtemp_2017-2019.csv", sep = ",")
>head(data.air)
1 ダウンロードした時刻:2019/03/20 12:34:56                                       
2                                                   羽田         羽田         羽田
3                                    年月日 平均気温(℃) 最高気温(℃) 最低気温(℃)
4                                  2017/1/1          8.6         12.5          3.7
5                                  2017/1/2          8.6         11.9          5.1
6                                  2017/1/3          8.9         13.6          3.7

これだと使いにくいのでヘッダーと何行スキップするかを指定する
また、ラベルを付け直して、一列目を日付のデータ形式に変換する
必要な列だけ取り出してデータの形式を変換する
気温のデータはファイルから読み込んだままだとfactor形式なので、
一度characterに変えてそれからnumericにする
手順についてググるといろいろ出てくる
【R】factor型について(2)_因子型(factor)から数値型(numeric)に戻すと数値が変わる? - 一所懸命に手抜きする

それと日付に関連してロケールを指定する必要があるらしい
なぜか?と聞かれると答えられないのだが。。。
コマンドはSys.setlocale()を使う

Sys.setlocale("LC_TIME", "C")
data.air <- read.table("airtemp_2017-2019.csv", header = TURE, skip = 2, sep = ",")
label.air <- c("date", "haneda_ave", "haneda_high", "haneda_low")
colnames(data.air) <- label.air
data.air$date <- as.POSIXct(as.Date(data.air$date, format = "%Y/%m/%d"))
data.air <- data.air[c("date", "haneda_ave")]
data.air[,1] <- as.Date(data.air[,1])
data.air[,2] <- as.numeric(as.character(data.air[,2]))

変換した結果がこれ

>head(data.air)
        date haneda_ave
1 2017-01-01        8.6
2 2017-01-02        8.6
3 2017-01-03        8.9
4 2017-01-04        9.5
5 2017-01-05        7.4
6 2017-01-06        5.8

海水温はテキストデータで手に入るが日付が
年・月・日で別の列になっていて、
ヘッダーと同じ文字列がデータの最後の行にも入っている

  yyyy mm dd areaNo. flag  Temp.
1 2017 01 01     704    R  25.16
2 2017 01 02     704    R  25.41
3 2017 01 03     704    R  25.58

ファイルを読み込んだら最後の行を取り除いて
年・月・日を一つにまとめたあとに日付のデータ形式に変換し、
必要な列だけ取り出す
こっちもfactorをcharacter -> numericの順に変換する

data.sea <- read.table("seatemp_2017-2019.txt", header = T, sep = ",")
data.sea <- data.sea[-nrow(data.sea),]
data.sea$date <- paste(data.sea[,1], data.sea[,2], data.sea[,3], sep = ".")
data.sea$date <- as.POSIXct(as.Date(data.sea$date, format = "%Y.%m.%d"))
data.sea <- data.sea[c("date", "Temp.")]
data.sea[,1] <- as.Date(data.sea[,1])
data.sea[,2] <- as.numeric(as.character(data.sea[,2]))

結果こうなる

        date Temp.
1 2017-01-01 25.16
2 2017-01-02 25.41
3 2017-01-03 25.58
4 2017-01-04 25.66
5 2017-01-05 25.66
6 2017-01-06 25.59

これを次の関数に投げるとグラフが出力される

PlotLineChart <- function(data.air, data.sea){
  # 日付の最小値と最大値を取得する
  Date.All <- as.Date(unique(append(data.air$date, data.sea$date)))
  Date.Range <- c(min(Date.All), max(Date.All))

  # マージンの調整、下、左、上、右の順
  par(mar = c(5, 5, 1, 5))

  # 1つ目のデータのプロット、軸を書かない
  # lは折れ線グラフ
  plot(data.air[,1], data.air[,2],
    xlab = "",
    ylab = "",
    xlim = Date.Range, # x軸の最小と最大
    ylim = c(-15, 35), # y軸の最小と最大
    type = "l",
    col = "gray",
    axes = FALSE
  )

  # 1つ目のデータのプロット、点を描く
  points(data.air[,1], data.air[,2], pch = 16, cex = 1, lwd = 1, col = "blue")

  # 軸ラベルを書く、4=右、3行目
  mtext("Air Temperature (Deg C)", side = 4, line = 3, cex = 1)
  axis(4, lwd = 1)

  # 重ね書きをする
  par(new = TRUE)

  # 軸だけ表示, type = "n"で指定する
  par(xaxt = "n")
  plot(0, 0,
    xlab = "Date (Year)",
    ylab = "Sea Temperature (Deg C)",
    type = "n",
    xlim = Date.Range, # x軸の最小と最大
    ylim = c(-15, 35), # y軸の最小と最大
    cex.lab = 1, # ラベルの文字の拡大率
    cex.axis = 1, # 軸の文字の拡大率
  )

  # x軸の表示
  par(xaxt = "s")
  axis.Date(1, at=seq(min(Date.All), max(Date.All), by = "years"), format = "%Y", cex.axis = 1)

  # 2つ目のデータのプロット、折れ線グラフで描く
  lines(data.sea[,1], data.sea[,2],
    cex = 1,
    col = "gray"
  )
  # 2つ目のデータのプロット、点を描く
  points(data.sea[,1], data.sea[,2],
    pch = 16,
    cex = 1,
    col = "red"
  )
}

日付の最大・最小を2つのデータを合わせてから取得しているのは
グラフの最大・最小を揃えるため

グラフをファイルに出力するにはこうする

png.name <- sprintf("temperature.png")
png(png.name, width = 600, height = 400)
PlotLineChart(data.air, data.sea)
dev.off()

カラーパレットを作って月毎にマーカーの色を変える事もできる
まずは関数を作る

AddColorCode <- function(d, n, p){
  # カラーパレットを作る。並びを逆にしている
  Func.Col.Pal <- colorRampPalette(rev(RColorBrewer::brewer.pal(8, p)))
  # 12色のカラーコードを取り出して、輝度値をnにする
  Col.Pal <- Func.Col.Pal(12)
  for(i in 1:length(Col.Pal)){
    Col.Pal[i] <- paste(Col.Pal[i], n, sep = "")
  }

  tmp.name <- colnames(d)

  # 月の列を追加
  tmp.month <- as.numeric(substring(d$date, 6, 7))
  d <- cbind(d, tmp.month)

  # カラーコードの作成
  tmp.col <- c()
  for(i in 1:nrow(d)){
    tmp.col <- append(tmp.col, Col.Pal[tmp.month[i]])
  }
  d <- cbind(d, tmp.col)

  tmp.add.name <- c("month", "color")
  colnames(d) <- append(tmp.name, tmp.add.name)

  return(d)
}

データの前処理のときにAddColorCode()を呼び出して
データにカラーコードを追加する
3番目の引数はRColorBrewer::display.brewer.all()で表示される
カラーパレットの名前を使う
なおRColorBrewerはパッケージなので使う前にインストールが必要

data.air <- AddColorCode(data.air, 20, "Paired")
data.air <- data.air[c("date", "haneda_ave", "color")]
data.air[,1] <- as.Date(data.air[,1])
data.air[,2] <- as.numeric(as.character(data.air[,2]))
data.air[,3] <- as.character(data.air[,3])

data.sea <- AddColorCode(data.sea, 80, "Spectral")
data.sea <- data.sea[c("date", "Temp.", "color")]
data.sea[,1] <- as.Date(data.sea[,1])
data.sea[,2] <- as.numeric(as.character(data.sea[,2]))
data.sea[,3] <- as.character(data.sea[,3])

PlotLineChart()で色を指定するcolの引数を2ヶ所修正する
ここと

  points(data.air[,1], data.air[,2], pch = 16, cex = 1, lwd = 1, col = data.air$color)

ここ

  points(data.sea[,1], data.sea[,2],
    pch = 16,
    cex = 1,
    col = data.sea$color
  )

修正後にPlotLineChart()でグラフを表示させるとこんな結果になる

見やすくなったとは言えないが。。。。
同じ色が同じ月なのでグラフによっては比較しやすくなるかな。

JBrowseにplug-inを追加する

JBrowseにpluginを追加する

インストール方法はここを参考にした。
How do I install a plugin
JBrowse FAQ - GMOD

試しにSashimiPlotのプラグインをインストールする
まずpluginのダウンロードして展開
ディレクトリの名前をリネームするして
JBrowseのディレクトリにコピーする

SashimiPlot 0.1.3
Releases · elsiklab/sashimiplot · GitHub

$ wget https://github.com/elsiklab/sashimiplot/archive/0.1.3.tar.gz
$ tar -zxvf 0.1.3.tar.gz
$ mv sashimiplot-0.1.3/ /var/www/html/jbrowse/plugins/SashimiPlot/

次に/var/www/html/jbrowse/jbrowse.confにプラグインの名前を追加する
よくわからないけどディレクトリの名前と同じにする必要があった

[GENERAL]
plugins += SashimiPlot

/var/www/html/jbrowse/setup.shを実行する

$ ./setup.sh

最後にtrackList.jsonを編集する
sashimiplot/README.md at master · elsiklab/sashimiplot · GitHub

SashimiPlotで表示したいbamファイルを記入する

  {
      "label": "hogehoge_junctions",
      "storeClass": "JBrowse/Store/SeqFeature/BAM",
      "type": "SashimiPlot/View/Track/Sashimi",
      "urlTemplate": "hogehoge.bam",
      "chunkSizeLimit": 50000000,
      "useXS": true
  }

ブラウザで表示した結果がこれ
上からBigWigのヒストグラム、SashimiPlot、ch38のgff

f:id:mecobalamin:20190320145923p:plain

刺し身には見えないけれども。
Why is it called sashimi_plot?
What is sashimi_plot? — MISO documentation

scriptを止める

Atomが起動しなくなったので
再インストール

ここからダウンロード
Atom

インストールしたらそのまま使えた
パッケージもインストールされてた
設定ファイルは残っているみたいだ

起動しなくなった理由だけど
Rのスクリプトをscript pluginで実行したら
スペルミスがあって中断できなくなった
それでAtom終了させたら
Java Script errorが出て起動しなくなった

再インストールで復活したから良かったけど
scriptで実行しているプロセスの終了はctrl-qだったらしい

Command and shortcut reference

Command macOS Linux/Windows Notes
Script: Run cmd-i shift-ctrl-b If text is selected a "Selection Based" is used instead of a "File Based" run
Script: Run by Line Number shift-cmd-j shift-ctrl-j If text is selected the line number will be the last
Script: Run Options shift-cmd-i shift-ctrl-alt-o Runs the selection or whole file with the given options
Script: Run with profile shift-cmd-k shift-ctrl-alt-b Runs the selection or whole file with the specified profile
Script: Close View esc or ctrl-w esc Closes the script view window
Script: Kill Process ctrl-c ctrl-q Kills the current script process

beautiful soupで更新チェック

iPad mini 5が出るような噂
もしかしたらmini 4が安く出るんではないかと
整備品の値段をappleのサイトでいちいち調べてた

名前と値段のリストをとってくるスクリプトを書いてみた
PythonのBeautiful Soup4をつかう

スクレイピングについて参考にしたのはこの本
PythonによるWebスクレイピング 第2版
(Amazonのページに飛びます)

まずiPadの正規整備品のサイトから
"refurbished-category-grid-no-js"のタグだけとってきて
その中に含まれる"li"のタグでfor文を回して
タグに含まれる製品名と価格を表示する

これだけだけど久々に書こうとしたら忘れてた

from urllib.request import urlopen
from bs4 import BeautifulSoup
import datetime

today = datetime.date.today()
todaydetail = datetime.datetime.today()

# date time
print('----------------------------------')
print("Today's refurbished iPad")
print(todaydetail)
print('----------------------------------')

url = "https://www.apple.com/jp/shop/refurbished/ipad"
html = urlopen(url).read()

bsObj = BeautifulSoup(html, 'lxml')

refurbishedipad = bsObj.find("div", {"class":"refurbished-category-grid-no-js"})

for i in refurbishedipad.find_all("li"):
    ipad_name = i.find("a").get_text()
    ipad_price = i.find("div", {"class":"as-price-currentprice as-producttile-currentprice"}).get_text()
    ipad_price = ipad_price.replace(" ", "").replace("\n", "")
    print(ipad_name + " " + ipad_price)

実行するとこんな感じ

----------------------------------
Today's refurbished iPad
2019-03-13 15:55:05.999841
----------------------------------
9.7インチiPad Pro Wi-Fi 32GB - スペースグレイ [整備済製品] ¥50,800(税別)
9.7インチiPad Pro Wi-Fi 32GB - ゴールド [整備済製品] ¥50,800(税別)
9.7インチiPad Pro Wi-Fi 32GB - ローズゴールド [整備済製品] ¥50,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - シルバー [整備済製品] ¥73,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - ゴールド [整備済製品] ¥73,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - ローズゴールド [整備済製品] ¥73,800(税別)
12.9インチiPad Pro Wi-Fi 64GB - ゴールド [整備済製品] ¥73,800(税別)
12.9インチiPad Pro Wi-Fi 256GB - ゴールド [整備済製品] ¥87,800(税別)

atomで実行すると文字化けする
power shellで大丈夫だから
atomの問題っぽい

Atom-runnerの日本語文字化け解決方法 | げーまーこんぶの戯言

process.env.PYTHONIOENCODING = "utf-8";

追加したけど改善してない。。。

14 March 2019追記
PCの再起動とかしてたら文字化けしなくなった
PythonだけでなくRも大丈夫
エラーコードは相変わらず文字化けする

28 November 2020 追記
Pythonのオブジェクト指向プログラミングを完全理解 - Qiita
こちらのサイトを参考に少し書き換えた

スクレイピングするclassを新たに作成して複数のページをスクレイピングできるようにした。
関数のままでも良かったのだが、勉強を兼ねて書き換え。

from urllib.request import urlopen
from bs4 import BeautifulSoup
import datetime

class HtmlScraping:
    def __init__(self, url):
        self.url = url

    def url_scraping(self):
        html = urlopen(self.url).read()
        bs_Obj = BeautifulSoup(html, 'html.parser')
        refurbished_product = bs_Obj.find("div", {"class":"refurbished-category-grid-no-js"})

        for i in refurbished_product.find_all("li"):
            product_name = i.find("a").get_text().replace("\xa0", "")
            product_price = i.find("div", {"class":"as-price-currentprice as-producttile-currentprice"}).get_text()
            product_price = product_price.replace(" ", "").replace("\n", "")
            print(product_name + " " + product_price)

if __name__ == '__main__':

    today = datetime.date.today()
    todaydetail = datetime.datetime.today()

    url_refurbished_products = {"Recommended_Products":"https://www.apple.com/jp/shop/refurbished",
                                "iPad":"https://www.apple.com/jp/shop/refurbished/ipad",
                                "mac":"https://www.apple.com/jp/shop/refurbished/mac"}

    for i in url_refurbished_products:
        print('----------------------------------')
        print("Today's refurbished " + i)
        print(todaydetail)
        print('----------------------------------')

        url = HtmlScraping(url_refurbished_products[i])
        url.url_scraping()

追記ここまで