Palmsonntagmorgen

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

記事一覧

Linux一般

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

環境変数PATHの通し方

文字列の検索・置換

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

解析環境構築

pyenvでPython環境を構築する

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

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

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

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

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

データ生成

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

2bit genome を作成する

genome tableを作成する

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 複製解析(出芽酵母)

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

【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】パッケージ管理コマンドあれこれ

今日はUbuntuのパッケージ管理についてです。

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コマンドでインストールすることができます。

# 例: 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

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環境の利用が楽ちんという話

今回も解析環境構築にまつわるお話です。
結論を先に書くと、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のインストール

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

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

Mac, Windowsのインストール方法は試していないので、ここでは私がLinuxへインストールした時の方法を載せておきます。
Singularityのversion 2であればaptで入りますが、最新のversion 3は自分でコンパイルする必要があります。

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

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

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のインストール

$ VERSION=1.12.7.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のインストール

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

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

感想

DockerよりSingularityの方が楽。

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

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

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

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

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

Dockerベースで提供されるパッケージの割合が目に見えて増えてきたように感じるので、簡単なまとめ。

Dockerとは

Dockerとはマシン上に違うマシンを立ち上げるための仮想化技術です。 WindowsMac上でLinux環境を立ち上げたり、Linux上に異なる複数のLinux環境を構築したりできます。
OSの仮想化技術としては VMwareVirtualBox が有名ですが、これらのホスト型仮想化は仮想OSの立ち上げ自体にPCリソースを多く取られるため、仮想OS内でのCPUやメモリ消費に制約がかかってしまうのが難点でした。一方Docker はコンテナ型仮想化という仕組みにより仮想化時のオーバーヘッドを抑えることができます。

更なるDockerの技術詳細については以下をご覧ください。

【図解】Dockerの全体像を理解する -前編- - Qiita

何がうれしいのか?

ひとつには、複数の物理的マシン上で完全に同一の環境を構築できることがあります。
マシンごとにOSのバージョンや構成、インストールされているツールは異なるので、あるツールを利用した時にエラーが起きたとしても、それがツールのバグなのか、環境依存の問題なのか、切り分けが難しいという問題がありました。 また、開発者が新しいパッケージをリリースした際に、開発者が想定する環境とユーザの実環境に食い違いがあり、開発者がユーザ環境ごとのエラー処理に追われる(質問受けに忙殺される)ということもしばしば起こります。しかし、ユーザに完全に同一の環境を作らせることは大変難しいことでした。 Dockerは環境そのものをイメージとして保存しますので、同じイメージを利用する限り、完全に同一の環境であることが保証されます。これはユーザにとってはインストールの苦労から解放されますし、開発者もエラートラブルを回避することができ、お互いにストレスは大幅に軽減されます。 また、仮にイメージ内で設定をいじってしまい再起不能になったとしても、イメージを再インストールすることでリセットでき、マシンを壊す心配はなくなります*1

もうひとつの利点は、ローカルの環境を汚さない、という点です。
Linux上に複数のパッケージをインストールする際、それぞれのパッケージが要求するライブラリのバージョンが異なることが原因で、同時にインストールすることができないということはよくあります (過去ブログ記事参照)。また、後からインストールしたパッケージが既存のライブラリを上書き更新(アップグレード、ダウングレード)してしまったために既存のツールがエラーで動かなくなるというような経験をされた方も多いのではないでしょうか。Dockerを使うと、それぞれのツールを個別にイメージ化することでこのコンフリクトはなくなります。

3つめは、新しくマシンを購入した際、あるいはOSを再インストールした際に、毎回環境構築をやり直さなくて済むという点です。 必要な時にイメージファイルをダウンロードすればよく、イメージファイルの更新は一元管理可能なので、管理者にとって利点が大きいです。

他にもいくつかありますが、大まかにはこんなところです。

Dockerのインストール

これは山ほど記事がありますので、ここではリンクを貼るのにとどめます。

LinuxにDockerをインストールする - Qiita

Ubuntuにdockerをインストールする - Qiita

Docker for MacMacターミナル上で homebrewを用いる方法があります。

Docker for Mac: DockerをMacにインストールする(更新: 2019/7/13) - Qiita

homebrew: 【初心者向け】MacにDockerを手軽にインストールする方法! | Aprico

Docker for Windowsを使うか、WSL上でインストールします。WSLはWSL2がおすすめです。

Docker for Windowsをインストール

【Windows10】Windows PreviewにリリースされたWSL2をインストールしてみた (7/9追記あり) - Palmsonntagmorgen

Windows HomeではHyper-Vが使えませんので、Docker Toolboxをインストールするか、WSL上でインストールします。
なお、WSL2自体がHyper-Vを使っているので、Windows Homeでは今のところWSL1しか使えません。

Docerk Toolbox: windows 10 home で docker を導入するメモ - Qiita

WSL1: WSL上でDocker Engineが動くようになっていたっぽいという話 - Qiita

僕はWSL1を使っていましたが、管理者権限でWSLを起動すればけっこういけました。 Windowsに関してはWSL2の本リリースがまだであることもあり、状況は流動的です。来年くらいにはスタンダードが固まるのではと予想しています。

動作確認

まずはDockerを起動します。起動方法は上のインストールのページ参照。
Dockerを起動する際にIDとパスワードが必要になるかもしれませんが、これはDockerHubsign upして作成したIDを入力してください。DockerHubは生成されたDockerイメージを共有するための仕組みです。

ターミナル上で以下のコマンドが動けば成功です。

 sudo docker run hello-world

*1:本当はDocker内からでもホストマシンを壊せますが、それについてはこのエントリのスコープ外ということで。