Palmsonntagmorgen

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

記事一覧

Linux一般

SSH公開鍵の生成・設定の方法

環境変数PATHの通し方

文字列の検索・置換

リダイレクトとパイプ

【Ubuntu】パッケージ管理コマンドあれこれ

 sortコマンドの使い方(ファイルのソート)

shuf: テキストファイルの行をランダム抽出

 

解析環境構築

pyenvでPython環境を構築する

バイオインフォマティクスのためのpython環境構築方法を考える

GitHubからプログラムをダウンロード・インストール

Docker・Singularity

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

【Windows10】Windows PreviewにリリースされたWSL2をインストールしてみた

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

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

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

データ生成

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する

2bit genome を作成する

genome tableを作成する

Gene annotation データを用意する(gtf形式)

gtfファイルからrefFlat形式への変換

LiftOver: BEDファイルを異なるbuildへ変換

fastqデータ取得

SRAからfastqを取得する

ENA,DDBJからfastqを取得する

ゲノムマッピング

Readをゲノムにマッピング (その1) (2017/12/19 追記あり)

Readをゲノムにマッピング (その2)

Readをゲノムにマッピング (その3) 圧縮ファイルを入力にする方法

マップデータ操作

SAMtoolsとリダイレクト

SAMtoolsワンライナー覚書

マッピング: CRAM形式を試す

ChIP-seq解析: DROMPA

DROMPA3: その1 インストール

DROMPA3: その2 parse2wig

DROMPA3: その3 ピーク抽出(peak calling)

DROMPA3: その4 マップリード分布の可視化その1

DROMPA3: その5 シェル変数を使う

DROMPA3: その6 ChIP/Input ratio 及び p値の可視化

DROMPA3: その7 -i オプション詳細

DROMPA3: その8 GVコマンドでのマクロな可視化

DROMPA3: その9 リードプロファイル

DROMPA3: その10 ヒートマップ

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

DROMPAplusを公開しました

ChIP-seq解析: 品質評価

Library complexity (PCR bias)とは何か

S/N比の評価手法 その1

S/N比の評価手法 その2 Cross-correlation profile

S/N比の評価手法 その3 deepTools

S/N比の評価手法 その4 SSP

ChIP-seq解析: ピークを入力とする操作

BEDtoolsワンライナー覚書

2サンプル間ピーク比較

RNA-seq解析:

RNA-seqによる発現量解析

HISAT-StringTie-Ballgown を試してみよう

HISAT-StringTie-Ballgown を試してみよう その2

STAR-RSEMによる発現量推定 その1

STAR-RSEMによる発現量推定 その2

STAR-RSEMによる発現量推定 その3

shuf: テキストファイルの行をランダム抽出

shufは入力行をシャッフルして出力するコマンドです。

【 shuf 】コマンド――入力行をシャッフルして出力する:Linux基本コマンドTips(112) - @IT

このshufに -nオプションを追加すると出力する行数を指定できます。これにより、たとえば「ピークファイル(BED)からランダムにn行抽出する」ことが極めて簡単にできます。

ピークのダウンロード

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 $ gunzip wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 # 内容の表示
 $ head wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
chr10   80917404        80917527        .       957     .       149.808548515806        -1      1.78532983501005        57
chr22   39570815        39571059        .       910     .       142.528642559654        -1      1.78532983501005        122
chr18   9643685 9643754 .       849     .       133.01812851997 -1      1.78532983501005        55
chr6    20811494        20811738        .       847     .       132.721884288658        -1      1.78532983501005        122
chr8    103918617       103918861       .       784     .       122.843985314485        -1      1.78532983501005        122
chr3    13692098        13692342        .       768     .       120.292626372352        -1      1.78532983501005        122
chr1    205263766       205264010       .       746     .       116.906226943541        -1      1.78532983501005        122
chr3    58035007        58035251        .       701     .       109.733323372861        -1      1.78532983501005        122
chr17   73781085        73781329        .       694     .       108.721733687985        -1      1.78532983501005        122
chr19   3968693 3968937 .       680     .       106.565969826384        -1      1.78532983501005        122

ピークファイルからランダムに100ピークを抽出

$  shuf -n 100 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak >  wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.random100.narrowPeak
$  head wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.random100.narrowPeak
chr19   46009794        46010038        .       140     .       21.9338179870502        -1      1.26616002142041        122
chr12   120242198       120242442       .       172     .       27.0065867360251        -1      1.33845649360463        122
chr3    185525441       185525685       .       131     .       20.6569949210623        -1      1.24621711968293        122
chr18   9643685 9643754 .       849     .       133.01812851997 -1      1.78532983501005        55
chr1    227750961       227751205       .       245     .       38.399370869934 -1      1.48650376307134        122
chr1    22379061        22379305        .       122     .       19.1627516856611        -1      1.24017410861089        122
chr19   47744334        47744578        .       203     .       31.8083079397841        -1      1.46637326118701        122
chr2    28114000        28114244        .       168     .       26.3872013177204        -1      1.33277877166033        122
chr4    113668168       113668412       .       114     .       17.9214793015699        -1      1.21539561863277        122
chr8    145133349       145133593       .       188     .       29.5813341454769        -1      1.42910033502578        122

ランダム抽出なので、実行するたびに結果は変わります。なので皆さんが試した場合、headで表示されるものはこれとは違うものになるでしょう。

これでpermutation testがやりやすくなりますね!

余談

ちなみに私がこのshufコマンドを知ったのはごく最近で、それまでは以下のような自作のperlスクリプトを使っていました。

#!/usr/bin/env perl

use strict;
use warnings;
use autodie;
die "random_pickup.pl <file> <p>\n" if($#ARGV !=1);

my $file=$ARGV[0];
my $p=$ARGV[1];

open(File, $file) || die "error: can't open $file.\n";

while(my $line = <File>){
    $x = rand();
    if($x < $p){ print $line; }
}
close (File);

これは確率pで各行を出力するというものです。これだと出力行数を明示的に指定できないという難点がありましたが、shufコマンドならば楽にできますね。 一方でfastqやsamファイルなど、単純なn行抽出では実行できないファイルをsubsampleしたい時は未だにperlを使っています。

sortコマンドの使い方(ファイルのソート)

以下の記事の中で、sortコマンドについて触れていました。

rnakato.hatenablog.jp

多くのコマンドはsorted BEDを要求しますので、事前に以下のコマンドで全てのBEDをソートしておくとストレスがないかと思います。

$ sort -k1,1 -k2,2n in.bed > in.sorted.bed

今日はこのコマンドが何をしているのかについて解説します。

sortについて

sortコマンドの使い方については以下のサイトがとてもわかりやすいです。

sortコマンドについて詳しくまとめました 【Linuxコマンド集】

ここではMACS2のピークファイルを使ってソートを実際に試してみましょう。

narrowPeakをsort

まずはサンプルとなるピークリストのダウンロードと解凍

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 $ gunzip wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 # 内容の表示
 $ head wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
chr10   80917404        80917527        .       957     .       149.808548515806        -1      1.78532983501005        57
chr22   39570815        39571059        .       910     .       142.528642559654        -1      1.78532983501005        122
chr18   9643685 9643754 .       849     .       133.01812851997 -1      1.78532983501005        55
chr6    20811494        20811738        .       847     .       132.721884288658        -1      1.78532983501005        122
chr8    103918617       103918861       .       784     .       122.843985314485        -1      1.78532983501005        122
chr3    13692098        13692342        .       768     .       120.292626372352        -1      1.78532983501005        122
chr1    205263766       205264010       .       746     .       116.906226943541        -1      1.78532983501005        122
chr3    58035007        58035251        .       701     .       109.733323372861        -1      1.78532983501005        122
chr17   73781085        73781329        .       694     .       108.721733687985        -1      1.78532983501005        122
chr19   3968693 3968937 .       680     .       106.565969826384        -1      1.78532983501005        122

narrowPeakフォーマットは、染色体、start, end, ピーク名(ID), score, strand, signalValue, p値、 q値、peak値(0-based offset from chromStart) からなります。

参考: Peak calling with MACS2 | Introduction to ChIP-Seq using high-performance computing

このサンプルは最初からピーク強度が強い順にソートされているようですね。

sort (デフォルト)

何もつけずにsortした結果が以下です。

$ sort wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr1    103716001       103716245       .       113     .       17.757293306914 -1      1.20462049496422        122
chr1    10739728        10739972        .       126     .       19.8727530712812        -1      1.23198670220658        122
chr1    108479111       108479355       .       139     .       21.9096721061029        -1      1.26200103555381        122
chr1    10856653        10856897        .       113     .       17.787081338271 -1      1.20987632444834        122
chr1    109359046       109359290       .       178     .       27.9341585774513        -1      1.37497420775533        122
chr1    110527136       110527380       .       121     .       19.033753770053 -1      1.23745570657267        122
chr1    111470122       111470366       .       120     .       18.8432010665144        -1      1.23273408572432        122
chr1    112946999       112947243       .       378     .       59.2558076715072        -1      1.53329986093226        122
chr1    113243465       113243709       .       114     .       17.9824166190626        -1      1.21252694920448        122
chr1    119623228       119623472       .       338     .       52.9747880171335        -1      1.74485156130976        122

左の列から優先で昇順でソートされるので、chr1が最初に来るようになります。ちなみにchr1の次はchr11です。

特定の列を基準にsort

仮に2列目の"start" を基準にソートしたい場合は、以下のように -k2,2 をつけます。
(-k2 だけだと2列目以降全ての列が評価されてしまいます。-k2,2 は「2列目から2列目まで」という指定です。)

$ sort -k2,2 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr14   100003838       100004082       .       162     .       25.3896151955265        -1      1.32856874236087        122
chr8    10009854        10010098        .       140     .       22.0324635890778        -1      1.26598847143077        122
chr15   100106095       100106339       .       165     .       25.8864092879034        -1      1.33959339080326        122
chr13   100150800       100151044       .       144     .       22.5409997536205        -1      1.27478892238628        122
chr15   100243871       100244115       .       287     .       45.0594193501982        -1      1.56124172764321        122
chr13   100340050       100340294       .       431     .       67.4827532683693        -1      1.68424674749618        122
chr13   100370580       100370824       .       134     .       20.9796948179283        -1      1.25368035557817        122
chr9    100394971       100395215       .       297     .       46.6277802709549        -1      1.52770570059048        122
chr7    100434051       100434295       .       208     .       32.5983076837242        -1      1.4658668221046 122
chr11   100449778       100450022       .       144     .       22.5396554016715        -1      1.27512069817163        122

なにか変ですね。"100003838"が最初に来ています。これは数値が文字列として評価されているためです。 数値としてソートしたい場合は -n をつけましょう。

$ sort -k2,2 -n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr16   122590  122834  .       156     .       24.5519243359907        -1      1.31865948396838        122
chr16   441938  442182  .       113     .       17.8207002619443        -1      1.20806676148905        122
chr9    444599  444843  .       191     .       30.0193503293631        -1      1.4152662953559 122
chr4    668097  668341  .       138     .       21.6323169668381        -1      1.25420898826128        122
chr16   790836  791080  .       122     .       19.1922172413508        -1      1.24286158989763        122
chr17   855058  855302  .       190     .       29.8288083711511        -1      1.41876633698804        122
chr1    936188  936432  .       278     .       43.6111912894822        -1      1.58313206719137        122
chr18   939233  939477  .       259     .       40.5622471525867        -1      1.48685535527511        122
chr7    985313  985557  .       133     .       20.8488735208324        -1      1.25369923175399        122
chr19   996416  996660  .       330     .       51.7245477629659        -1      1.65959931242927        122

無事、数値としてソートされました。

ちなみに -k2,2 -n-k2,2n と書いても同じ意味になります。

$ sort -k2,2n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr16   122590  122834  .       156     .       24.5519243359907        -1      1.31865948396838        122
chr16   441938  442182  .       113     .       17.8207002619443        -1      1.20806676148905        122
chr9    444599  444843  .       191     .       30.0193503293631        -1      1.4152662953559 122
chr4    668097  668341  .       138     .       21.6323169668381        -1      1.25420898826128        122
chr16   790836  791080  .       122     .       19.1922172413508        -1      1.24286158989763        122
chr17   855058  855302  .       190     .       29.8288083711511        -1      1.41876633698804        122
chr1    936188  936432  .       278     .       43.6111912894822        -1      1.58313206719137        122
chr18   939233  939477  .       259     .       40.5622471525867        -1      1.48685535527511        122
chr7    985313  985557  .       133     .       20.8488735208324        -1      1.25369923175399        122
chr19   996416  996660  .       330     .       51.7245477629659        -1      1.65959931242927        122

降順でソートしたい場合は -r オプションをつけましょう。

$ sort -k2,2 -n -r wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr1    244436392       244436636       .       116     .       18.2779430757098        -1      1.22768405980467        122
chr1    244327430       244327674       .       236     .       37.0352452557144        -1      1.49136169383242        122
chr1    244090184       244090428       .       198     .       31.1366875436721        -1      1.45062610008222        122
chr1    242097512       242097756       .       219     .       34.351615073636 -1      1.46008791542184        122
chr1    239979463       239979707       .       144     .       22.5476334083494        -1      1.27445689294964        122
chr1    238811236       238811480       .       187     .       29.3108873795414        -1      1.42587671098019        122
chr2    238522973       238523217       .       206     .       32.3223759856553        -1      1.47135994807226        122
chr1    237360863       237361107       .       129     .       20.2738825008399        -1      1.24189631661435        122
chr2    237118258       237118502       .       183     .       28.6979797945465        -1      1.3746686662588 122
chr2    236688138       236688382       .       154     .       24.1282500315595        -1      1.31204005182904        122

-k2,2 -n -r-k2,2nr と書くこともできます。 -k2,2n -r とは書けないようです。

染色体配置にもとづくソート

染色体上での並び順にもとづいてピークをソートしたい場合、まず1列目と2列目でソート(1列目優先)となります。これは以下のように書けます。

$ sort -k1,1 -k2,2n wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr1    936188  936432  .       278     .       43.6111912894822        -1      1.58313206719137        122
chr1    1790095 1790339 .       147     .       23.0760153956621        -1      1.27356690138484        122
chr1    2143910 2144154 .       328     .       51.3985090086356        -1      1.67414642060339        122
chr1    3530010 3530254 .       195     .       30.5372043198762        -1      1.41709407864103        122
chr1    6333445 6333689 .       116     .       18.1876867422188        -1      1.22125052352078        122
chr1    6844488 6844732 .       167     .       26.1507228133338        -1      1.3422008190137 122
chr1    7669033 7669277 .       132     .       20.7725171056509        -1      1.24639610451395        122
chr1    7764587 7764831 .       125     .       19.6798804866856        -1      1.24285094687436        122
chr1    8933714 8933958 .       209     .       32.8047056184361        -1      1.45819850201426        122
chr1    9488912 9489156 .       300     .       47.0204903548675        -1      1.50753840495668        122

-k オプションを複数指定すると、左から優先でファイルがソートされます。1列目は文字列で2列目が数値なので、-k2,2の方だけに-n オプションがついています。
これが冒頭で示した、BEDtoolsで要求されるソートになります。

ピーク強度でソート

ピークが強い順にtop 1000ピークを抽出したい、というような場合は、以下のようになります。
MACS2の場合はq値がピーク強度を表しているので、9列目でソートします。

$ sort -k9,9nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak |head
chr1    156643538       156643782       .       584     .       91.4369629380883        -1      1.78532983501005        122
chr1    205263766       205264010       .       746     .       116.906226943541        -1      1.78532983501005        122
chr1    214489937       214490181       .       556     .       87.1055294122947        -1      1.78532983501005        122
chr1    30447209        30447453        .       623     .       97.5499934901461        -1      1.78532983501005        122
chr10   65389313        65389557        .       649     .       101.626596581112        -1      1.78532983501005        122
chr10   80917404        80917527        .       957     .       149.808548515806        -1      1.78532983501005        57
chr10   92079336        92079580        .       569     .       89.1329212789991        -1      1.78532983501005        122
chr11   65264905        65265149        .       547     .       85.6244766911001        -1      1.78532983501005        122
chr11   86451724        86451968        .       565     .       88.4367426823928        -1      1.78532983501005        122
chr11   9634736 9634980 .       586     .       91.8571732244745        -1      1.78532983501005        122

同じq値をもつピークが(何故か)多くあり、同値でうまくソートできませんでした。

では、signalValueを表す7列目を追加で指定してみましょう。

$ sort -k9,9nr -k7,7nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head
chr10   80917404        80917527        .       957     .       149.808548515806        -1      1.78532983501005        57
chr22   39570815        39571059        .       910     .       142.528642559654        -1      1.78532983501005        122
chr18   9643685 9643754 .       849     .       133.01812851997 -1      1.78532983501005        55
chr6    20811494        20811738        .       847     .       132.721884288658        -1      1.78532983501005        122
chr8    103918617       103918861       .       784     .       122.843985314485        -1      1.78532983501005        122
chr3    13692098        13692342        .       768     .       120.292626372352        -1      1.78532983501005        122
chr1    205263766       205264010       .       746     .       116.906226943541        -1      1.78532983501005        122
chr3    58035007        58035251        .       701     .       109.733323372861        -1      1.78532983501005        122
chr17   73781085        73781329        .       694     .       108.721733687985        -1      1.78532983501005        122
chr19   3968693 3968937 .       680     .       106.565969826384        -1      1.78532983501005        122

無事、思い通りのソートになりました。 sort -k9,9nr -k7,7nr wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak | head -n1000 とすれば、top1000ピークが手に入ります(ヘッダ行が含まれない場合。)

注意点

ソートする文字列の中に記号や日本語など複雑な文字列が含まれている場合のソート順は環境に依存します。複雑なファイルのソートの場合、場合によっては使ったOSなどでsortの結果が違うかもしれないので、ファイル比較などをする際は統一環境でソートしておくほうが無難です。

あけましておめでとうございます

本年も私とこのブログをよろしくお願いします。

と言いつつ、これからこのブログをどうしようかなあ、などと考えています。。 このブログを始めた当初はNGS解析のノウハウを記した日本語のWebサイトはRNA-seq関連くらいしかなく、ChIP-seqはMACSの使い方くらいしかなかったので、自分のツールの紹介も兼ねてこのページを始めたのですが、最近はQiitaにアドベントカレンダーができるくらいたくさんの方が解説ブログなど出されているので、相対的に需要は減った気がします。

書きたいトピックはまだ色々とあるものの、何しろ自ラボの運営に一生懸命になっており、また(ありがたいことに)総説や講演依頼も年に複数いただくので、ネタはそちらで使ってしまう。というのもあり、このブログでしかできないことは減っている感じです。こんな感じの雑文ならいくらでも書けるのですが、需要はある…?ない…?ないよな。

今後も覚え書きや新しいツールにトライしてみた。系の記事は書きたいと思っているので、もしご興味があれば今後ともお付き合いください。

DROMPAplusを公開しました

少し前になりますが、DROMPA3のアップデートであるDROMPAplusを正式にリリースしました。

DROMPAplus: a pipeline tool for ChIP-seq analysis — DROMPAplus 1.4.0 documentation

今回はその記事です。

Reference

4月に以下の論文を公開しました。ChIP-seq解析の総説なのですが、中でDROMPAplusについて紹介しています。 この論文を「DROMPAplusのpaper」ということにしようと思っているので、今後もしDROMPAplusを論文で使用してくださる場合はこちらの論文をcitationしていただけると嬉しいです。

www.sciencedirect.com

DROMPAplusの特長

DROMPA3はCで書かれていましたが、DROMPAplusはC++でスクラッチから書き直しました。 そのためDROMPA3の機能のうちまだ未実装のものもありますが、新しく追加されている機能もあります。

ここでは現時点での最新版(version 1.8.1)について記載しています。

Fragment lengthを自動で推定するようになった

DROMPA3はfragment lengthをpredefinedの固定長にしていましたが、DROMPAplusは内部でSSPを動かすことでsingle-end readでもfragment lengthを自動で推定します。 サンプルのマップリード数が極めて少ない場合などにはうまく推定できなくなりますので、得られた推定長が妥当かスタッツファイルで確認するようにしてください。

外部のwigファイルを入力にできるようになった

DROMPA3は入力のwig/Bigwig/BedGraphを読み込む時に一緒に総リード数などをまとめたスタッツファイルを読み込むため、parse2wigで生成されたwigファイルしか実質的に入力にできませんでした。 DROMPAplusではスタッツファイルのない(外部ツールで生成された)wigファイルを入力にした場合、初回時に自動でスタッツファイルを新たに生成し、2回目以降はそれを読み込むようになります。

2サンプルをoverlayできるようになった

DROMPA3でも隠れオプションで存在していたのですが、DROMPAplusで正式な機能になりました。

f:id:rnakato:20200622120724j:plain

こんな感じ。昨年出した血管内皮論文では、H3K4me3とH3K27acをoverlayすることでエンハンサー・プロモーター領域を視認できるようにしました。 (透明度もオプションで変更できます)

Hi-Cループを可視化できるようになった

f:id:rnakato:20200622121046j:plain

こんな感じ。 DROMPA3ではChIA-PETループ(Mango) しか対応していませんでしたが、DROMPAplusではHICCUPS(Juicertools)のループフォーマットにも対応しています。なお、色はループの重要度を表しています。

負の値を可視化できるようになった

図は省略。read distribuionは正の値しか持たないため負の値は対象にしていませんでしたが、Compartment scoreなどのBedGraphファイルを可視化するために負の値も可視化できるようにしました。

染色体構造の可視化

f:id:rnakato:20200622121704j:plain

この図上部の染色体の図のことです。GVでこれを表示したいとずっと思っていたのですが、Ideogramという名前で保存されていることがわかったので、この度オプションで可視化できるようにしました。

正規化後のDistributionをWigで出力できるようになった

DROMPAplusで可視化されている ChIP/Input enrichment, p-valueなどは全てwig/Bigwig/BedGraph形式で出力可能です。 これらはスムージングやspike-in normalizationなどユーザが指定したオプションを反映したかたちで出力されるため、他のツールを利用した下流解析が楽になりました。

PROFILEコマンドで .tsv ファイルを出力するようになった

PROFILEでaveraged profileを可視化する時は今までpdfとRスクリプトしか出力していませんでしたが、一緒に元データの.tsvファイルも出力するようにしました。そのため出力ファイルのサイズが大きくなっていますが、このtsvファイルを使って色々下流解析をすることが可能です。

これに伴いDROMPA3のHEATMAPコマンドは廃止し、drompa.heatmap.py という名前のPythonスクリプトを新たに追加しました。このスクリプトに上のtsvファイルを読み込ませることで今まで通りHeatmapを生成できます。

Install

詳細はこちら:

https://drompaplus.readthedocs.io/en/latest/content/Install.html#

ソースファイルからインストールする時は以下のようにします。git cloneする時に --recursiveオプションをつけることに注意してください。

DROMPAplusはCMakeを使っているのですが、内部で利用しているhtslibライブラリがMakefileで動いているためMakefileでラップしています。 以下のようにするとインストール可能です。

(version1.5.0ではMacコンパイルエラーが出ましたが、最新バージョンでは解消されているはずです。もしエラーが出ましたらご一報を。)

git clone --recursive https://github.com/rnakato/DROMPAplus
cd DROMPAplus
make

Docker

Dockerイメージも提供しています。以下のようにすれば実行可能です。

docker pull rnakato/ssp_drompa
docker run -it --rm rnakato/ssp_drompa drompa+

終わりに

言語をC++に変更したので拡張の自由度が広がりました。今後もどんどん機能を追加していく予定ですので、リクエストあればお知らせください。

旧版のDROMPA3も引き続き利用可能ですが、アップデートはバグ修正など必要最小限に留める予定ですので、現在DROMPA3を利用されている方は是非DROMPAplusへの移行をお願いします。

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、とするのがよいかもしれませんね。