Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

Monocle3をRstudioで起動できる dorowu/ubuntu Dockerイメージ

前回の記事↓では、ブラウザからアクセスできるLinux GUI のDockerイメージを紹介しました。

rnakato.hatenablog.jp

このイメージを使って、Rstudio内でMonocle3をGUI起動できるイメージを作ったので紹介したいと思います。

イメージのダウンロード・起動

以下のコマンドでDockerイメージをダウンロードします。

$ docker pull rnakato/monocle3

その後、以下のコマンドで起動し、表示されるURLをブラウザに貼り付けることでGUIを起動します。

$ docker run --rm -it -p 6079:80 rnakato/monocle3

f:id:rnakato:20200302092334j:plain
GUI起動画面

このようなLinuxのデスクトップ画面が開けます。

RStudioの起動

左下のメニュー内「プログラミング」の中にRStudioのリンクがあります。

(なおRStudioでなく普通のRを呼び出すこともできます。)

f:id:rnakato:20200424144054j:plain

Rstudioが起動しました。

f:id:rnakato:20200424144201j:plain

Monocle3呼び出し

このイメージにはMonocle3が既にインストールされていますので、 library(monocle3)でmonocle3を呼び出すことができます。

f:id:rnakato:20200424144339j:plain

Monocle3のチュートリアルに沿って実行していくと、以下のようになります。

f:id:rnakato:20200424144351j:plain

注意点

注意点として、このDockerイメージ外部でコピーしたテキストをDocker内でペーストすることはできませんので、チュートリアルのページをこのDocker内で閲覧するなどして、コマンドを貼り付けていく必要があります。 また、他のDockerイメージと同じく、コンテナを終了すると作業内容は消えてしまいますので、データを永続化したい場合はローカルフォルダと同期させるなどしてください。

おわりに

Monocle3はインストールが若干トリッキーなので、インストールに苦労している方もおられると思います。 チュートリアルを進めるくらいであればこのイメージで完結するかと思いますので、利用してみてください。 またMonocle3は作業の一部でGUIを必要とするため、Jupyterではうまく作業できない部分があります。そういう時にこのGUIイメージは有用になると思います。

もしご興味があれば使ってみてください。

docker-ubuntu-vnc-desktopを使ってDockerイメージ (ssp_drompa) をGUIで動かす

更新間隔があいてしまいすみません。 論文書いたり博論審査したり、色々大変です。

今日はDockerのお話です。 最近以下のようなツイートを見つけました。

(ちなみにこのAtsushi SakaiさんはPythonRoboticsの作者です。)

DockerとSingularityについてはこのブログでもエントリを書いたので参照してください。

【Win】【Mac】【Linux】Dockerのインストール 【2019年7月現在】 - Palmsonntagmorgen

Singularityを使ったDocker環境の利用が楽ちんという話 - Palmsonntagmorgen

自分のシングルセル解析DockerイメージをGUIにできれば面白いのでは?と思い、試してみました。

ssp_drompa (GUI) のビルド

シングルセルDockerはイメージが大きくて構築に時間がかかるため、ここではssp_drompaを使って構築します。

以下がDockerfileの全体です。もともとdebianベースのため、基本的には FROMdebianからdorowu-bionicに変えただけです。 rnakato/ubuntu となっていますが、これはdorowu/ubuntu-desktop-lxde-vncにいくつかのパッケージをプリインストールしたもので、基本的に同じと思ってもらって大丈夫です。

FROM rnakato/ubuntu:dorowu-bionic
LABEL maintainer "Ryuichiro Nakato <rnakato@iam.u-tokyo.ac.jp>"

WORKDIR /home
ENV DEBIAN_FRONTEND=noninteractive

RUN apt-get update \
    && apt-get install -y --no-install-recommends \
    build-essential \
    ca-certificates \
    git \
    libboost-all-dev \
    libgsl-dev \
    libgtk2.0-dev \
    libgtkmm-3.0-dev \
    libz-dev \
    r-base \
    samtools \
    && apt-get clean \
    && rm -rf /var/lib/apt/lists/*
RUN git clone https://github.com/rnakato/SSP.git \
    && cd SSP \
    && make
RUN git clone https://github.com/rnakato/DROMPA3 \
    && cd DROMPA3 \
    && make
RUN git clone --recursive https://github.com/rnakato/DROMPAplus \
    && cd DROMPAplus \
    && git submodule foreach git pull origin master \
    && make
ADD script script

ENV PATH ${PATH}:/home/SSP/bin:/home/DROMPA3:/home/DROMPAplus/bin:/home/DROMPAplus/submodules/cpdf/Linux-Intel-64bit:/home/DROMPAplus/otherbins:/home/script

CMD ["/bin/bash"]

これを rnakato/ssp_drompa:dorowu-bionic という名前でビルドしました。

ssp_drompa (GUI) の実行

上でビルドしたdockerを起動し、ブラウザからアクセスしてみます。
(ビルドしたものは公開していますので、誰でも以下のコマンドで利用可能です。)

$ docker run --rm -it -p 6079:80 rnakato/ssp_drompa:dorowu-bionic

実行環境によると思いますが、上のコマンドを実行すると以下のようなログが続いてどわーっと表示されるかと思います。

2020-03-02 00:13:46,603 CRIT Supervisor running as root (no user in config file)
2020-03-02 00:13:46,603 WARN Included extra file "/etc/supervisor/conf.d/supervisord.conf" during parsing
2020-03-02 00:13:46,617 INFO RPC interface 'supervisor' initialized
2020-03-02 00:13:46,617 CRIT Server 'unix_http_server' running without any HTTP authentication checking
2020-03-02 00:13:46,617 INFO supervisord started with pid 11
2020-03-02 00:13:47,619 INFO spawned: 'nginx' with pid 14
2020-03-02 00:13:47,629 INFO spawned: 'web' with pid 15
2020-03-02 00:13:47,630 INFO spawned: 'novnc' with pid 16
2020-03-02 00:13:47,631 INFO spawned: 'wm' with pid 17
2020-03-02 00:13:47,632 INFO spawned: 'pcmanfm' with pid 18
2020-03-02 00:13:47,634 INFO spawned: 'lxpanel' with pid 19
2020-03-02 00:13:47,637 INFO spawned: 'xvfb' with pid 20
2020-03-02 00:13:47,646 INFO spawned: 'x11vnc' with pid 21
2020-03-02 00:13:47,902 INFO  Listening on http://localhost:6079 (run.py:87)
2020-03-02 00:13:48,715 INFO success: nginx entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,715 INFO success: web entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,715 INFO success: novnc entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,716 INFO success: wm entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,716 INFO success: pcmanfm entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,716 INFO success: lxpanel entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,717 INFO success: xvfb entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
2020-03-02 00:13:48,717 INFO success: x11vnc entered RUNNING state, process has stayed up for > than 1 seconds (startsecs)
127.0.0.1 - - [2020-03-02 00:14:16] "GET /api/health HTTP/1.1" 200 122 0.162835
127.0.0.1 - - [2020-03-02 00:14:46] "GET /api/health HTTP/1.1" 200 122 0.143152

途中に 2020-03-02 00:13:47,902 INFO Listening on http://localhost:6079 (run.py:87) というURLが表示されています。このURL (http://localhost:6079) をブラウザに貼り付けるとGUIが開けます。
(6079というのがポート番号です。この数字が6080だったり、その他の数字になっている場合は、dockerをrunするコマンドの -p 6079:80 のポート番号を同じ数字に変更して再実行してください。)

f:id:rnakato:20200302092334j:plain
GUI起動画面

ブラウザでLinuxのデスクトップ環境が開けます。これだけでもテンションが上がりますね!

今回はssp_drompaですので、ターミナルを起動してsspとdrompa+が動くか試してみます。 (スタートメニューから LXTerminal を選択するとターミナルが起動します。)

f:id:rnakato:20200302093204j:plain
ターミナル画面

ターミナル上でちゃんと起動しています。

まとめ

DockerイメージをGUIで起動できると、Rの作図やpdfファイル描画などで「Xウィンドウが開かない」という問題も解決するので、Jupyterを起動するのに比べてもかなり直感的で使いやすくなる印象です。 ちゃんと試していませんが、イメージ内で完結する作業であればほぼ問題なく動作するのではないでしょうか。
GUIが必要な作業はこれで、ターミナル上でコマンド的に使いたい場合はSingularity、とするのがよいかもしれませんね。

【Linux】リダイレクトとパイプ

リダイレクト

> 記号を使うことで、ターミナルに出力される文字列をファイルに保存することができ、これをリダイレクトといいます。

$ echo aaa bbb
aaa bbb
$ echo aaa bbb > text.txt   # echo aaa bbb  の結果を text.txt に保存する
$ cat text.txt
aaa bbb

もともと text.txt が存在していた場合、ファイルの内容は上書きされます。既存ファイルの末尾に追記したい場合は >> とふたつ重ねて指定します。

$ echo aaa bbb > text.txt   # echo aaa bbb  の結果を text.txt に保存する
$ echo ccc ddd >> text.txt   # echo ccc ddd  の結果を text.txt に追記する
$ cat text.txt
aaa bbb
ccc ddd

パイプ

> の代わりに | を使うと、左で実行されたコマンドの結果を直接右のコマンドの入力にすることができます。
例えば以下のコマンドは、カレントディレクトリの bam で終わるファイルのうち、 ファイル名にMCF7を含むものだけを表示します。

$ ls -l *bam | grep MCF7

パイプを使いこなすと色々な作業を一行で書けるようになります。 このような一行の便利コマンドをワンライナーと呼んだりします。

  # workdirに含まれるBEDファイルの数を表示
  $ ls workdir/*.bed | wc -l 
  # peak.bedに含まれるピーク数を染色体別にカウント
  $ cut -f1 peak.bed | sort | uniq -c
  # カレントディレクトリのsraファイルを全てfastqに展開 (xargs 利用)
  $ ls *.sra | xargs fastq-dump

標準出力とエラー出力

ターミナルに表示される出力には、標準出力の他に、エラーメッセージ用のエラー出力というものがあります。 出力結果が標準出力かエラー出力かはツール依存です。例えば bowtie のマッピングスタッツはエラー出力になっています。

上の >| を使った操作は標準出力のみを対象にしています。エラー出力も含めてファイルに追記したい場合は command >& file と書きます。

参考: Bashの入出力リダイレクトまとめ - Qiita

DROMPA3: その11 複製解析(出芽酵母)

久しぶりのDROMPAのエントリです。今回は出芽酵母(S. cerevisiae)の複製解析を行います。複製開始後のDNAを複製前のDNAで割り算することで、ゲノム上のどこでどの程度複製フォークが進んでいるかを可視化することができます。 データは以下の論文のものを用います。

Origin association of Sld3, Sld7, and Cdc45 proteins is a key step for determination of origin-firing timing. - PubMed - NCBI

ここで用いるサンプルは、SOLiDで生成されたcolor spaceデータになります。 color spaceデータはbowtieでマッピング可能です(専用のindexファイルが必要)。

インストール

DROMPA3のインストール方法はこの記事を参照してください。

データダウンロード

今回は、必要なreference data類をGoogle Driveにアップしてみましたので、以下のコマンドでダウンロードしてください。

wget -O genome_table "https://drive.google.com/uc?export=download&id=11xhQpJt-kYJKc38HWXUfwoim2utSby8Q"
wget -O ARS-oriDB.txt "https://drive.google.com/uc?export=download&id=12A8EigD3J_DB_BaMyDRAd51LlOBg9F9u"
wget -O SGD_features.tab "https://drive.google.com/uc?export=download&id=13OPXenZyxXmF8hkQ1my3GlW6QrH7Y1Xr"

bowtieのcsインデックスはbowtie-buildで作成可能ですが、ここでは作成済のものをダウンロードします。

wget "https://drive.google.com/uc?export=download&id=138QOuZUjOMyXeFVs6Uv6SbpyKixt9N-G" -O S_cerevisiae_cs.zip
unzip S_cerevisiae_cs.zip 

csfastqファイルのダウンロード

以下ではシェル変数とfor文を用いていきます。よくわからないという方は DROMPA3: その5 シェル変数を使う - Palmsonntagmorgen を参考にしてください。

今回用いるデータをSRAからダウンロードします。 以下では fasterq-dump を使っていますが、これは fastq-dump の後継で、マルチコアで高速にfastqファイルをダウンロードするものです。fastq-dump を使ってダウンロードしてもかまいません。

for num in $(seq 398609 398624)
do
     fasterq-dump SRR$num
done

手元の環境ではダウンロード完了まで4時間ほどかかりました。 私の場合は通常、ファイルサイズを削減するためにダウンロードしてから pigz で圧縮するのですが、今回は実行時間節約のため圧縮せずにいきます。(pigzsudo apt install pigz でインストール可能です。)

ファイルを確認しましょう。

$ ls -lh SRR*
-rw-r--r-- 1 rnakato rnakato 7.3G 116 11:40 SRR398609.fastq
-rw-r--r-- 1 rnakato rnakato 6.7G 116 11:49 SRR398610.fastq
-rw-r--r-- 1 rnakato rnakato 7.4G 116 12:06 SRR398611.fastq
-rw-r--r-- 1 rnakato rnakato 5.5G 116 12:17 SRR398612.fastq
-rw-r--r-- 1 rnakato rnakato 8.2G 116 12:32 SRR398613.fastq
-rw-r--r-- 1 rnakato rnakato 6.0G 116 12:42 SRR398614.fastq
-rw-r--r-- 1 rnakato rnakato 8.2G 116 12:55 SRR398615.fastq
-rw-r--r-- 1 rnakato rnakato 6.8G 116 13:07 SRR398616.fastq
-rw-r--r-- 1 rnakato rnakato 6.3G 116 13:19 SRR398617.fastq
-rw-r--r-- 1 rnakato rnakato 7.2G 116 13:32 SRR398618.fastq
-rw-r--r-- 1 rnakato rnakato 6.8G 116 13:43 SRR398619.fastq
-rw-r--r-- 1 rnakato rnakato 6.0G 116 13:54 SRR398620.fastq
-rw-r--r-- 1 rnakato rnakato 7.0G 116 14:19 SRR398621.fastq
-rw-r--r-- 1 rnakato rnakato 6.4G 116 14:35 SRR398622.fastq
-rw-r--r-- 1 rnakato rnakato 7.1G 116 14:54 SRR398623.fastq
-rw-r--r-- 1 rnakato rnakato 7.1G 116 15:14 SRR398624.fastq

複製前・複製後のペアが8ペア、計16サンプルです。

fastqファイルの中身を見てみましょう。 head コマンドで最初の数行を表示します。

$ head SRR398609.fastq
@SRR398609.1 2010_018_bcSample1_1_19_677 length=50
T0002.2000..000..0000.00000000..0..02...0.........0
+SRR398609.1 2010_018_bcSample1_1_19_677 length=50
!''*+!%'*'!!)&)!!*(-)!.-)'/*%,!!.!!)&!!!)!!!!!!!!!8
@SRR398609.2 2010_018_bcSample1_1_20_852 length=50
T0100.1010..000..0200.00000020..0..00...0.........1
+SRR398609.2 2010_018_bcSample1_1_20_852 length=50
!%)).!%(1%!!(&'!!4%))!))*-)%%)!!1!!**!!!*!!!!!!!!!&
@SRR398609.3 2010_018_bcSample1_1_20_996 length=50
T3102.0001..102..2220.20001230..2..23...3.........2

fastqファイルは 4行で1つのリードを表します。 @ で始まる行がリード名、次の T で始まる行が リード配列、4行目は読まれた各塩基のクオリティ(信頼度)を表します。 (このデータでは、3行目は意味を持ちません。)

リード配列が 0,1,2 の数字で表されていますね。 これが color space データ で、【2つの塩基の組み合わせ】を数字で表現したものです。 SOLiD sequencer で読まれた配列はこのような形式で出力されます。 color space形式のfastqファイルはcsfastq形式と呼ばれます。

マッピング

以下のfor文は、各csfastqファイルに対して、bowtieでマッピングした後そのままsamtoolsで sorted bamに変換しています。詳細は SAMtoolsとリダイレクト - Palmsonntagmorgen を参照してください。

mkdir -p bam
for num in $(seq 398609 398624)
do
    bowtie -S -C S_cerevisiae-cs/S_cerevisiae-cs \
     SRR$num.fastq -n2 -k1 --chunkmbs 2048 -p12 \
    | samtools sort \
    > bam/SRR$num-n2-k1.sort.bam
done

最初にbam という名前のディレクトリを作成し、マップファイルをその中に出力するようにしています。

\ は1行のコマンドを複数行に分割して書く時に使う記号です。

-S は出力ファイルをSAM形式にするオプション、 -C はcsfastqデータである時に指定します。

S_cerevisiae-cs/S_cerevisiae-cs はさきほどダウンロード・解凍したcolor space用のbowtieインデックスファイルです。

-n2 はリード配列のミスマッチを2まで許容し、 -k1 はリードがゲノムの複数箇所にマッチした時にベストマッチのみ出力するオプションです。

--chunkmbs 2048 はbowtieに許容するメモリ上限、 -p 12 はCPUを12コア使うことを意味します。

bedGraphファイルの変換

bamファイルをbedGraphに変換します。genome_table は先ほどダウンロードしたものを用います。 -of 3 が出力をbedGraphにするオプションです。

for num in $(seq 398609 398624)
do
   parse2wig -f BAM -i bam/SRR$num-n2-k1.sort.bam \
   -o SRR$num-n2-k1 -n GR -of 3 -binsize 100 -gt genome_table  
done

このコマンドの結果、parse2wigdir の中にbedGraphファイルが生成されます。

pdf可視化(1サンプル)

それでは可視化してみましょう。 今回は「複製後/複製前」の割り算 (ratio) を可視化したいので、 PC_ENRICH コマンドを使います。 サンプル数が多いので、まず1サンプルペアのみ可視化します。 アノテーションファイルは先ほどダウンロードしたものを指定してください。

dir=parse2wigdir 
drompa_draw PC_ENRICH \
-gene SGD_features.tab -gftype 3 \
-ars ARS-oriDB.txt \
-gt genome_table \
-i $dir/SRR398612-n2-k1-scer-raw-mpbl-GR,$dir/SRR398611-n2-k1-scer-raw-mpbl-GR,YST1019_Gal \
-p pdf/yeast1 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \
-if 3

-i で入力サンプルを指定します。<ChIP>,<Control>,<Label> の順にコンマで区切って指定します。 今回の場合だと、 <複製後>,<複製前>,<株名> となります。 入力サンプルには染色体までのprefixを指定します。

-p 出力ファイル名のprefix
-sm 1000 スムージング幅 (ここでは1kbp)
-ls 200 1段あたりの length
-rmchr 染色体別のファイルを作成しない
-nosig ピーク抽出を行わない
-lpp 2 で1ページあたりの段数を2に   -bn 3 でy軸の罫線の数を3に
-ystep 14 で罫線1あたりの長さ(pixel)を14に(デフォルトは20)
-scale_ratio 1 で罫線1あたりのスケールを1に
それぞれ変更しています。結果は以下のようになります。

evince pdf/yeast1.pdf

f:id:rnakato:20191102181354j:plain
Repliseq1

Enrichmentが2倍付近まで到達している領域が、既に複製されたゲノム領域を表します。複製開始点を中心に、両側に複製フォークが広がっていく様子がわかります。

pdfファイルの作成(複数サンプル)

では、全サンプルを同時に可視化してみましょう。

dir=parse2wigdir 
drompa_draw PC_ENRICH \
-ars ARS-oriDB.txt \
-gt genome_table \
-i $dir/SRR398612-n2-k1-scer-raw-mpbl-GR,$dir/SRR398611-n2-k1-scer-raw-mpbl-GR,YST1019_Gal \
-i $dir/SRR398610-n2-k1-scer-raw-mpbl-GR,$dir/SRR398609-n2-k1-scer-raw-mpbl-GR,YST1019_Raf \
-i $dir/SRR398616-n2-k1-scer-raw-mpbl-GR,$dir/SRR398615-n2-k1-scer-raw-mpbl-GR,YST1053_Gal \
-i $dir/SRR398614-n2-k1-scer-raw-mpbl-GR,$dir/SRR398613-n2-k1-scer-raw-mpbl-GR,YST1053_Raf \
-i $dir/SRR398620-n2-k1-scer-raw-mpbl-GR,$dir/SRR398618-n2-k1-scer-raw-mpbl-GR,YST1076_Gal \
-i $dir/SRR398619-n2-k1-scer-raw-mpbl-GR,$dir/SRR398617-n2-k1-scer-raw-mpbl-GR,YST1076_Raf \
-i $dir/SRR398624-n2-k1-scer-raw-mpbl-GR,$dir/SRR398623-n2-k1-scer-raw-mpbl-GR,YST1287_Gal \
-i $dir/SRR398622-n2-k1-scer-raw-mpbl-GR,$dir/SRR398621-n2-k1-scer-raw-mpbl-GR,YST1287_Raf \
-p pdf/yeast2 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \
-if 3

各ペアを個別に -i で指定しています。また、今回の複製解析では遺伝子情報は必要ないので、 -gene オプションは省略しました。 結果は以下のようになります。

f:id:rnakato:20191102181527j:plain
Repliseq2

上のコマンドはちょっと長いので、変数を使って短くしてみましょう。得られる結果は yeast3 と同じです。

# サンプルの定義
dir=parse2wigdir
post=-n2-k1-scer-raw-mpbl-GR   # サンプルに共通するpostfixを格納
IP1_60="$dir/SRR398612$post"  # YST1019 Gal 60min
IP1_0="$dir/SRR398611$post"   # YST1019 Gal 0min
IP2_60="$dir/SRR398610$post"  # YST1019 Raf 60min
IP2_0="$dir/SRR398609$post"   # YST1019 Raf 0min
IP3_60="$dir/SRR398616$post"  # YST1053 Gal 60min
IP3_0="$dir/SRR398615$post"   # YST1053 Gal 0min
IP4_60="$dir/SRR398614$post"  # YST1053 Raf 60min
IP4_0="$dir/SRR398613$post"   # YST1053 Raf 0min
IP5_60="$dir/SRR398620$post"  # YST1076 Gal 60min
IP5_0="$dir/SRR398618$post"   # YST1076 Gal 0min
IP6_60="$dir/SRR398619$post"  # YST1076 Raf 60min
IP6_0="$dir/SRR398617$post"   # YST1076 Raf 0min
IP7_60="$dir/SRR398624$post"  # YST1287 Gal 60min
IP7_0="$dir/SRR398623$post"   # YST1287 Gal 0min
IP8_60="$dir/SRR398622$post"  # YST1287 Raf 60min
IP8_0="$dir/SRR398621$post"   # YST1287 Raf 0min

# サンプルペアの定義
s1="-i $IP1_60,$IP1_0,YST1019_Gal"
s2="-i $IP2_60,$IP2_0,YST1019_Raf"
s3="-i $IP3_60,$IP3_0,YST1053_Gal"
s4="-i $IP4_60,$IP4_0,YST1053_Raf"
s5="-i $IP5_60,$IP5_0,YST1076_Gal"
s6="-i $IP6_60,$IP6_0,YST1076_Raf"
s7="-i $IP7_60,$IP7_0,YST1287_Gal"
s8="-i $IP8_60,$IP8_0,YST1287_Raf"

# DROMPAの実行
drompa_draw PC_ENRICH \
   -ars ARS-oriDB.txt \
   -gt genome_table \
   $s1 $s2 $s3 $s4 $s5 $s6 $s7 $s8 \
   -p pdf/yeast3 -sm 1000 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -nosig \
-if 3

サンプルを定義するのにだいぶ行数を使ってしまいましたが、このようにファイル名をわかりやすい変数名に変換することで、 誤字やサンプルの取り違えなどのミスをある程度予防することができます。

-nosig を消すと、Enrichment >= 2.0 の部分が赤でハイライトされます。 (ハイライトの閾値-ethre オプションで指定可能です)

drompa_draw PC_ENRICH -gftype 3 -ars $ARS -ter $TER -gt $gt \
   $s1 $s2 $s3 $s4 $s5 $s6 $s7 $s8 \
   -p pdf/yeast4 -ls 200 -lpp 2 -bn 3 -ystep 14 -scale_ratio 1 -rmchr -sm 1000 \
   -if 3

f:id:rnakato:20191102182154j:plain
Repliseq3

このようにすると、複製されている領域のサンプル間の差がより顕著になりますね。

ピーク抽出

サンプル間で複製状態を比較するために、 yeast4において赤でハイライトされた領域(ピーク領域)のリストが欲しい、という場合もあるでしょう。 リスト形式でピークを出力するには drompa_peakcall を使います。

#(変数の行は省略)
drompa_peakcall PC_ENRICH $s1 -p YST1019_Gal -gt $gt -if 3

こうすると、 YST1019_Gal.bed と YST1019_Gal.xls の2つのファイルが生成されます。 YST1019_Gal.bed はBED形式のファイル、 YST1019_Gal.xls はより詳細な情報が書かれたファイルになります。

まとめ

今回は複製解析について紹介しました。 酵母の場合はChIP-seqでも同様に PC_ENRICHで可視化します。これはヒトやマウスに比べるとウィンドウあたりのリード数が十分多く、割り算してもnoisyにならないという理由によります。 PC_ENRICHはデフォルトで 統計検定はせず、シンプルにenrichmentとそのリード数だけを閾値に用います。有意差検定を用いたピーク抽出をする場合は通常の PC_SHARP コマンドを使ってください。

【Linux】【Ubuntu】パッケージ管理コマンドあれこれ(11/21追記)(2021/6/18追記)

今日はUbuntuのパッケージ管理についてです。
(11/21 alien について追記)

apt / apt-get / aptitude

全て同じです。

最初にあったのは apt-get です。これはUbuntuの親にあたるDebianにおけるパッケージ管理システムです。 その後 aptitudeに移行しようとしたものの失敗したようです(?) 現在は apt に移行していますが、 apt-get も引き続き使えるので、 apt-getのままになっているマニュアルもあります。

以下、実行例です。管理者権限が必要なので、sudoをつけています。

$ sudo apt update       # データベースの更新
$ sudo apt upgrade      # インストール済のライブラリをアップグレード
$ sudo apt install xxx  # xxxをインストール
$ sudo apt remove xxx   # xxxをアンインストール
$ apt-cache search xxx  # xxxというパッケージを検索
$ sudo apt show xxx     # パッケージの詳細情報を表示

dpkg

.deb フォーマットのファイルからパッケージをインストールします。 aptで提供されていないパッケージでも .deb ファイルが提供されていれば、dpkgコマンドでインストールすることができます。

(2021/6/18追記:現在はaptコマンドでも.debファイルからインストールできるようにアップデートされています。)

# 例: Google chromeをインストール
$ wget https://dl.google.com/linux/direct/google-chrome-stable_current_amd64.deb
$ sudo dpkg -i google-chrome-stable_current_amd64.deb

add-apt-repository

公式のaptリポジトリに含まれないツールを追加したい場合、aptで参照されるリポジトリadd-apt-repository で追加することができます。Rの最新版などもこの方法でインストールします。 追加する場合は提供元が信頼できるかよく確認してください。

# 例:Google Driveを追加
$ sudo add-apt-repository ppa:alessandro-strada/ppa
$ sudo apt-get update
$ sudo apt-get install google-drive-ocamlfuse

「apt keyが信用できない」とエラーが出たら、以下のコマンドでキーを登録します。

$ sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys <key>

snap

Ubuntu 独自に開発されたパッケージ管理システムで、必要なライブラリを他のソフトウェアに影響を与えないように、一つのファイルシステムにまとめてインストールを行う仕組みです。 aptと併用することが可能ですが、いずれでもインストールできるパッケージに関しては今のところ apt を使った方が良いようです。 なお、2019年10月時点では WSL(Windows Subsystem for Linux) はサポートされていません。

# snap で Skypeをインストール
$ sudo snap install skype --classic

alien (11/21追記)

CentOS7 向けの パッケージファイル形式として rpm がありますが、Alien というコマンドを用いると rpmdeb に変換することができます。 aliensudo apt install alien でインストールできます。

sudo alien xxx.rpm # xxx.rpmからxxx.debを作成
sudo dpkg -i xxx.deb

rpmdebなど色々な形式がありますが、要はパッケージのソースをある種の形式に圧縮しているだけなので、相互変換はそれほど難しくありません。

apt-file

特定のファイルが収録されているパッケージを探します。

あるツールをインストールしようとして、たとえば「zlib.hがない」とエラーになった場合に、どのパッケージを apt でインストールすればよいか調べたい場合は以下のようにします。

$ apt-file update  # データベースの更新
$ apt-file search zlib.h   # zlib.h を含むパッケージの検索
...   (大量の結果が表示される)
$ apt-file search include/zlib.h  # 結果を絞り込むために include/ をつける
cc65: /usr/share/cc65/include/zlib.h
libklibc-dev: /usr/lib/klibc/include/zlib.h
libz-mingw-w64-dev: /usr/i686-w64-mingw32/include/zlib.h
libz-mingw-w64-dev: /usr/x86_64-w64-mingw32/include/zlib.h
python-pycparser: /usr/share/python-pycparser/fake_libc_include/zlib.h
python3-pycparser: /usr/share/python3-pycparser/fake_libc_include/zlib.h
zlib1g-dev: /usr/include/zlib.h

いくつかありますが、最もパスがデフォルトっぽい zlib1g-dev あたりをインストールするのがよさそうです。

Linuxbrew

MacでおなじみのパッケージマネージャーHomebrew の Linux版です。 下記コマンドでインストールできます。

$ sh -c "$(curl -fsSL https://raw.githubusercontent.com/Linuxbrew/install/master/install.sh)"

このコマンドを打つと自動的に/home/linuxbrew/ ディレクトリが作成され、/home/linuxbrew/.linuxbrew/bin にインストールされます。 (PATHを通す必要があります) 異なるディレクトリにインストールすることも可能ですが、Warningが出るので、このディレクトリにインストールすることが推奨されているようです。

Linuxbrew を使うと ghq などがインストールできるようになります。 pyenv もインストールできるのですが、インストールディレクトリは $HOME/.pyenv になるので、パッケージ管理としては微妙かも?

$ brew install ghq
$ brew install pyenv

【Linux】文字列の検索・置換

現在、研究室用の新人教育ページを作っているのですが、せっかく作ったので一部転載。 文字列を検索・置換するLinuxコマンドです。

grep: 文字列の検索

対象ファイルから特定の文字列を含む行だけを表示します。

$ grep hogehoge sample.txt  # sample.txtの中からhogehogeを含む行を表示
$ grep -v hogehoge sample.txt  # sample.txtの中からhogehogeを ”含まない” 行を表示
$ grep -n hogehoge sample.txt  # 行番号をつけて表示
$ grep -3 hogehoge sample.txt  # hogehogeを含む行±3行を表示
$ grep -n3 hogehoge sample.txt  # hogehogeを含む行±3行を行番号をつけて表示

また、ファイルだけでなく以下のような使い方もできます。| はリダイレクトを表し、 ls コマンドの出力結果に対して grep を使っています。

$ ls workdir/* | grep Rad21 # workdir内のファイルのうち Rad21 を含むファイルを表示
$ ls workdir/* | grep txt$  # workdir内のファイル名がtxtで終わるファイルを表示

sed: 文字列の置換

以下は test.txt内の aaa を AAA に変換し、modified.txt に保存する例です。

$ cat test.txt
  aaa
  bbb
  ccc
$ sed -e 's/aaa/AAA/g' test.txt > modified.txt
$ cat modified.txt
  AAA
  bbb
  ccc

高度な使い方もできます。

$ sed -n -e 1p # 先頭の行のみを出力
$ sed -n -e \$p # 最後の行のみを出力, \ はシェルのエスケープ
$ sed -n 5,\$p # 5行目~最後の行まで
$ sed -n -e 6,15p # 6行目から15行目を出力
$ sed -n -e 1~2p # 奇数行のみを出力
$ sed -n -e 1~5p # 1行目、6行目、11行目、16行目、、、を出力
$ sed -i -e '10,10d' # ファイルの10行目を消す
$ sed -n -e /xxx/p # 正規表現にマッチする行を出力, " grep -e xxx" と同じ
$ sed -e /xxx/d # 逆に正規表現にマッチしない行を出力, "grep -v -e xxx" と同じ
$ sed -n -e '\%xxx%p' # 正規表現にマッチする行を出力, 先頭に \ を付ければ、正規表現を囲む記号はなんでもよい
$ sed -n -e /xxx/,+3p # 正規表現にマッチする行から3行分を出力, 正規表現にマッチするごとに最低3行が出力される
$ sed -n -e 20,+3p # 20行目から23行目までの4行分を出力

参考: 指定行削除するコマンド(sed - それマグで!

tr: 文字列の置換 (1文字単位)

trは1文字単位で文字列を置換します。sedの方が高機能ですが、ワンライナーで改行をスペースに変換したい場合などに威力を発揮します。

$ cat test.txt
aaa
bbb
ccc
# 以下のコマンドは catで test.txtの内容をリダイレクトで trに送っている。結果は modified.txtに保存。
$ cat test.txt | tr '\n' ' ' > modified.txt
$ cat modified.txt
aaa bbb ccc

以下のような使い方もできます。

$ cat text1.txt | tr ',' '\n' > text2.txt  # カンマを改行に変換して text2.txtに保存
$ cat text1.txt | tr ' ' '\t' > text2.txt  # スペースをタブに変換して text2.txtに保存
$ tr -d "\r" text1.txt > text2.txt   # Windows形式の改行コードをLinux形式に変換し text2.txtに保存
$ tr -s " " text1.txt > text2.txt   # text1.txtの連続したスペースを1つに変換して text2.txtに保存
$ tr -s "\r" text1.txt > text2.txt   # text1.txtの連続した改行を1つに変換して text2.txtに保存

Singularityを使ったDocker環境の利用が楽ちんという話 (2022/3/13, 2022/4/12追記)

今回も解析環境構築にまつわるお話です。
結論を先に書くと、Docker使うならSingularityオススメ

Singularityとは

7月に書いた下記エントリでは、Dockerを使うメリットについて簡単に説明しました。

rnakato.hatenablog.jp

一方、Dockerにはいくつか不満な点が残ります。

  • 実行にsudo権限が必要である。共有サーバの場合、全ユーザにsudo権限を与えるのはセキュリティ的に問題。
  • sudo権限を与えられていない環境 (スパコンなど)では利用が難しい。
  • sudoで実行するため(デフォルトのコマンドでは)コマンド実行で生成されるファイルの所有者がrootになる。
  • コマンドを実行するためにまずDockerコンテナの生成・削除が必要なのがやや面倒。
  • ホストディレクトリとコンテナ内ディレクトリの紐づけ(マウント)がやや面倒。

JupyterやSQLデータベースのようにバックグラウンドで起動させる用途であれば、サーバ起動時に管理者がコンテナ起動しておけば良いだけなのでそれほど問題にはなりませんが、NGS解析のように各ユーザがアクティブにコマンド実行してデータ生成・処理をするような場合にはsudoまわりの問題は厄介です。

そこで登場するのがSingularityです。

www.ecomottblog.com

上記のページがわかりやすいですが、Docker単体と比較したSingularityのメリットは以下のようなものです。

  • Dockerのイメージを利用可能
  • sudo不要(ユーザが引き継がれる)
  • コンテナ起動不要
  • カレントディレクトリが自動でマウントされる
  • GUIGPUなどデバイス回りの対応が楽

私の管理するサーバのユーザに試してもらったところ、ユーザごとに追加の環境設定をしてもらう必要もなく、Dockerの概要を深く説明する手間もなく、通常のコマンドとほぼ同じ感じで利用可能でした。なので、NGS解析用共用サーバの管理に悩まされている管理者の方には特にお勧めします。 たとえば以下のような悩みから解放されます。

rnakato.hatenablog.jp

Singularityを試してみる

Singularityのインストールは多少面倒なのでひとまず後回しにし、ここではSingularityをインストール済という前提で、使用例を示します。

Dockerイメージのダウンロード

まず、singularity pull コマンドでDockerイメージをローカルに保存します。
DockerイメージはDocker Hubなどのレポジトリに登録されている必要があります。 ここでは例として私の作ったssp_drompa を使います*1。これはSSP, DROMPA3, DROMPAplusがインストール済のUbuntuイメージです。

$ singularity pull ssp_drompa.img docker://rnakato/ssp_drompa
INFO:    Starting build...
Getting image source signatures
Skipping fetch of repeat blob sha256:5b7339215d1d5f8e68622d584a224f60339f5bef41dbd74330d081e912f0cddd
(中略)
Writing manifest to image destination
Storing signatures
INFO:    Creating SIF file...
INFO:    Build complete: ssp_drompa.img

ssp_drompa.img という名前でイメージファイルが保存されました(ファイル名は自由に決められます)。 イメージファイルはLinux環境を丸ごと保存したものであり、ファイルサイズはそれなりに大きいです。特に大きいイメージだと10GBを超えるものもあります。 ssp_drompa.imgは473Mなので、やや小さめですね。

$ ls -lh
合計 473M
-rwxrwxrwx 1 rnakato rnakato 473M  821 16:19 ssp_drompa.img

イメージの利用

ここではsspを例にして実験します。まずは、普通にローカルのsspを実行してみます。

$ ssp
ssp: command not found

おそらく多くの方はsspをインストールしていないと思いますので*2、not foundになったのではないでしょうか。

次に、Singularityのイメージを使ってコマンドを実行します。 コマンドは、以下のようになります。

singularity exec <イメージファイル> <コマンド>

$ singularity exec ssp_drompa.img ssp

SSP v1.1.3
===============

Usage: ssp [option] -i <inputfile> -o <output> --gt <genome_table>
Use --help option for more information on the other options

イメージ内のSSPが起動し、ヘルプが表示されました。

優秀な読者の方は既にSSPをインストールしてくれているかもしれません。その場合はローカルのSSPと区別するために which ssp してみましょう。

$ singularity exec ssp_drompa.img which ssp
/home/SSP/bin/ssp

/home/SSP/bin/ssp というパスが表示されれば、イメージ内のsspになっている証拠です。

重要なポイントは、「あるツールを利用したいのだが、ローカルPCにはまだインストールされていないという時に、イメージファイルをダウンロードするだけで利用可能な状態になる」という点です。通常ですとローカルPCにそのツールをインストールしよう!という流れになる訳ですが、何らかの理由でエラーになり、インストールをあきらめてしまうということはよくあると思います(冒頭のDockerエントリ参照)。イメージのダウンロード形式に切り替えることでインストールのコストが不要となり、すぐにツールを導入できます。Python2系と3系の共存もこれで問題なくなります。

上記のメリットはDockerでも共通ですが、Singularityの場合は実行したいコマンドの前に singularity exec <イメージファイル> を追加するだけですので、初心者の方でも技術的に難しいことはありません。PATHを追加で通す必要もなくなるので、むしろ簡単になるかもしれません*3

Dockerの場合

同じことをDockerでやる場合は以下のようになります。

$ sudo docker pull rnakato/ssp_drompa  # イメージのダウンロード
$ sudo docker run -it --rm  rnakato/ssp_drompa ssp

(本来はコンテナをバックグラウンドで起動し、そのコンテナを指定してコマンドを実行するという二段階のステップを踏む必要がありますが、ここでは簡単のため「コマンド実行のためにコンテナを起動し、コマンド終了と同時にコンテナを削除する」というオプションにしています。)

Dockerの場合はsudo権限が必要になります。ユーザをdockerグループに所属させることでsudoなしで実行できるようになりますが*4、管理者権限で実行しているという点は変わりません。従って、Dockerコマンドによって生成されるファイルの所有者はrootになります。
Dockerのコマンドで生成されるファイルを一般ユーザ権限にするためには -u $(id -u):$(id -g) -v /etc/group:/etc/group:ro -v /etc/passwd:/etc/passwd:ro -v /etc/shadow:/etc/shadow:ro のようなオプションを追加してやる必要がありますが、やや面倒です。

ホストディレクトリのマウント

Singularityの場合はカレントディレクトリを自動でマウントしますので、カレントディレクトリ以下にあるファイルはそのまま参照できます。異なるディレクトリのファイルやツールを利用したい場合は別途 --bind オプションでマウントしてやる必要があります。

例として、/work ディレクトリをマウントする場合は以下のようになります。

$ singularity exec --bind /work ssp_drompa.img ssp -i /work/sample.bam --gt /work/genome_table

Dockerですと全てのホストディレクトリをマウントする必要があります。コマンドは以下のようになります。

$ sudo docker run -it --rm -v $(pwd):$(pwd) -v /work:/work rnakato/ssp_drompa ssp -i /work/sample.bam --gt /work/genome_table

-v オプションでマウントします。ディレクトリのパスは絶対パスである必要があります。

Singularityの注意点

  • Docker Hub上でイメージが更新された場合、それを反映するためにはイメージをダウンロードしなおす必要があります。
  • イメージは通常バージョン管理されていますので、バージョン番号をイメージファイル名につけておくと良いと思います。
  • ダウンロードしたイメージ内で更にパッケージをインストールするなど、イメージの内容をローカルで追加編集したい場合はSingularityでもsudo 権限が必要になります。そのようなケースでは、DockerとSingularityを併用するのがよいのではないかと思います。

Singularityのインストール

最後になりましたが、インストール方法です。

以下が公式のマニュアルです。

Linux, Windows (WSL2)では以下の方法でインストール可能です。
Singularityのversion 2であればaptで入りますが、最新のversion 3は自分でコンパイルする必要があります。

なお、2019年8月時点では、WSL2だとインストールには成功するもののコマンドの実行は ERROR : Failed to set securebits: Invalid argument というエラーが出てうまくいかないようでした。

(2022/3/13追記) パッケージのバージョンを最新に更新しました。現在のSingularity v3.8.4 はWindows11 + WSL2 + Ubuntu 20.04 で問題なく動作するようになっています。

必要なパッケージのインストール

sudo apt-get update && sudo apt-get install -y \
   build-essential \
   libssl-dev \
   uuid-dev \
   libgpgme11-dev \
   squashfs-tools \
   libseccomp-dev \
   wget \
   pkg-config \
   git

Goのインストール

Goの最新バージョンは こちら https://go.dev/dl/ で確認可能です。Goはlinuxbrewでもインストール可能です。

$ VERSION=1.17.8.linux-amd64
$ wget https://dl.google.com/go/go$VERSION.tar.gz 
$ sudo tar -C /usr/local -xzvf go$VERSION.tar.gz 
$ rm go$VERSION.tar.gz
$ echo 'export GOPATH=${HOME}/go' >> ~/.bashrc && \
 echo 'export PATH=/usr/local/go/bin:${PATH}:${GOPATH}/bin' >> ~/.bashrc && \
 source ~/.bashrc

Singularityのインストール

Singularityの最新バージョンは https://github.com/apptainer/singularity/blob/master/INSTALL.md を参照してください。

git clone https://github.com/sylabs/singularity.git 
cd singularity
git checkout v3.8.4
./mconfig 
make -C ./builddir 
sudo make -C ./builddir install

参考:Singularityのインストール (v3) - Qiita

Macへのインストール (2022/4/12 追記)

brew install --cask virtualbox vagrant vagrant-manager
mkdir vm-singularity-ce && cd vm-singularity-ce
vagrant destroy && rm Vagrantfile
export VM=sylabs/singularity-ce-3.8-ubuntu-bionic64 && \
    vagrant init $VM && \
    vagrant up && \
    vagrant ssh

感想

DockerよりSingularityの方が楽。

*1:https://hub.docker.com/r/rnakato/ssp_drompa

*2:使ってくれているという方は、どうもありがとうございます。

*3:実際、初心者にとってはインストール自体よりパスを通すことの方が難しく感じることもあります。

*4:https://qiita.com/matyapiro31/items/3e6398ce737e2cdb5a22