Palmsonntagmorgen

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

記事一覧

解析環境構築

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

環境変数PATHの通し方

pyenvでPython環境を構築する

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

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

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

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

データ生成

常染色体と性染色体のみのゲノム配列ファイル 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 ヒートマップ

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

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内からでもホストマシンを壊せますが、それについてはこのエントリのスコープ外ということで。

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

みなさんgtfファイルからrefFlatに変換する時ってどうされてるんですかね?Rを使っている?
自分は自作のツール "gtf2refFlat" を使っているので、ここではそれを紹介します。

gtf形式については↓をご覧ください。

https://bi.biopapyrus.jp/rnaseq/mapping/gtf.html

gtfファイルは最も詳細な遺伝子アノテーションと思いますが、1行1遺伝子にはなっていないので使いづらいことがあります。
私が作ったツール群はだいたい遺伝子アノテーションファイルをrefFlat形式で受け付けているのですが、EnsemblはrefFlat形式を提供していないので、自作しました。

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

Ensemblサイトから、主要な染色体の遺伝子のみ格納した Homo_sapiens.GRCh38.97.chr.gtf.gz をダウンロードします。

 # ダウンロード
 $ wget ftp://ftp.ensembl.org/pub/release-97/gtf/homo_sapiens/Homo_sapiens.GRCh38.97.chr.gtf.gz
 # 解凍
 $ gunzip Homo_sapiens.GRCh38.97.chr.gtf.gz
 # 確認
 $ head Homo_sapiens.GRCh38.97.chr.gtf
#!genome-build GRCh38.p12
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.27
#!genebuild-last-updated 2019-03
1       havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
1       havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; tag "basic"; transcript_support_level "1";
1       havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       havana  exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";

ChIPseqToolsのインストール

gtf2refFlatはChIPseqToolsの中に含まれています。 以下はgitがインストールされている前提です。gitについては以下の記事を参照ください。

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

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

ChIPseqTools/bin の中に gtf2refFlat が生成されていればOKです。 ヘルプを表示してみましょう。

 $ ChIPseqTools/bin/gtf2refFlat

Options:
  -g [ --gtf ] arg      Gene file
  -n [ --name ]         Output name instead of id
  -u [ --unique ]       Only output one transcript per one gene (default: all
                        transcripts)
  -h [ --help ]         Print this message

transcript単位のrefFlatを出力

gtf2refFlatはデフォルトではisoformを個別に格納したrefFlatを出力します。

 $ ChIPseqTools/bin/gtf2refFlat -g Homo_sapiens.GRCh38.97.chr.gtf \
   > Homo_sapiens.GRCh38.97.transcript.refFlat
 # 確認
 $  head -n5 Homo_sapiens.GRCh38.97.transcript.refFlat
ENSG00000210196 ENST00000387461 chrMT   -       15956   16023   15956   16023   1       15956,  16023,  Mt_tRNA Mt_tRNA
ENSG00000198804 ENST00000361624 chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding protein_coding
ENSG00000209082 ENST00000386347 chrMT   +       3230    3304    3230    3304    1       3230,   3304,   Mt_tRNA Mt_tRNA
ENSG00000210077 ENST00000387342 chrMT   +       1602    1670    1602    1670    1       1602,   1670,   Mt_tRNA Mt_tRNA
ENSG00000210127 ENST00000387392 chrMT   -       5587    5655    5587    5655    1       5587,   5655,   Mt_tRNA Mt_tRNA

1列目はgene ID, 2列目がtranscript IDです。一番右の列には遺伝子のタイプが2列表示されていますが、右から2列目がgene type, 最右列がtranscript typeです。coding geneでもtranscriptがnoncodingであることはあり得ますので、このようになっています。

Ensembl IDでなく transcript name で出力したい場合は -n オプションをつけます。

$ ChIPseqTools/bin/gtf2refFlat -n -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat
$ head -n5 Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat
MT-TP   MT-TP-201       chrMT   -       15956   16023   15956   16023   1       15956,  16023,  Mt_tRNA Mt_tRNA
MT-CO1  MT-CO1-201      chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding  protein_coding
MT-TL1  MT-TL1-201      chrMT   +       3230    3304    3230    3304    1       3230,   3304,   Mt_tRNA Mt_tRNA
MT-TV   MT-TV-201       chrMT   +       1602    1670    1602    1670    1       1602,   1670,   Mt_tRNA Mt_tRNA
MT-TA   MT-TA-201       chrMT   -       5587    5655    5587    5655    1       5587,   5655,   Mt_tRNA Mt_tRNA

トップにミトコンドリアが来ているのがアレですが、とにかく出力できています。

Gene単位のrefFlatを出力

transcript単位ではなくgene単位でrefFlatを出力したい場合は -u オプションをつけます。遺伝子名にしたい場合は、transcriptの時と同様に -n オプションをつけます。

 $ ChIPseqTools/bin/gtf2refFlat -u -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.gene.chr.refFlat
 $ head -n5  Homo_sapiens.GRCh38.97.gene.chr.refFlat
ENSG00000270016 ENST00000602595 chr15   -       30607695        30608193        30607695        30608193        1      30607695,        30608193,       lncRNA  lncRNA
ENSG00000271983 ENST00000607458 chr15   -       80693216        80693707        80693216        80693707        1      80693216,        80693707,       lncRNA  lncRNA
ENSG00000222139 ENST00000410207 chr15   +       80663772        80663878        80663772        80663878        1      80663772,        80663878,       snRNA   snRNA
ENSG00000259175 ENST00000558208 chr15   -       80554609        80562944        80554609        80562944        3      80562862,80560056,80554609,      80562944,80560167,80556933,     lncRNA  lncRNA
ENSG00000184140 ENST00000328882 chr15   +       101803509       101806887       101805720       101806658       2      101803509,101805687,     101803545,101806887,    protein_coding  protein_coding
 # 遺伝子名を出力する場合
 $ ChIPseqTools/bin/gtf2refFlat -u -n -g Homo_sapiens.GRCh38.97.chr.gtf > Homo_sapiens.GRCh38.97.gene.chr.name.refFlat
 $ head -n5 Homo_sapiens.GRCh38.97.gene.chr.name.refFlat
AC026150.3      AC026150.3-201  chr15   -       30607695        30608193        30607695        30608193        1      30607695,        30608193,       lncRNA  lncRNA
AC023302.1      AC023302.1-201  chr15   -       80693216        80693707        80693216        80693707        1      80693216,        80693707,       lncRNA  lncRNA
RNU6-380P       RNU6-380P-201   chr15   +       80663772        80663878        80663772        80663878        1      80663772,        80663878,       snRNA   snRNA
AC108451.1      AC108451.1-201  chr15   -       80554609        80562944        80554609        80562944        3      80562862,80560056,80554609,      80562944,80560167,80556933,     lncRNA  lncRNA
OR4F6   OR4F6-201       chr15   +       101803509       101806887       101805720       101806658       2       101803509,101805687,    101803545,101806887,    protein_coding  protein_coding

遺伝子ごとに出力する場合、gtf2refFlatはgtfに含まれるtranscriptの中で最もcanonicalなもの (Consensus CDS (CCDS) タグを持つtranscriptの中で最も遺伝子長が長いもの)を出力します。左から2列目がそのtranscriptのID(または名前)になります。

おまけ:coding geneだけを出力したい場合

gawkを用いてcoding gene だけを出力してみます。一番右のアノテーションが"protein_coding"となっている行だけを出力するというワンライナーです。

$ cat Homo_sapiens.GRCh38.97.transcript.name.chr.refFlat \
  | gawk '{ if ($13=="protein_coding") print }' \
  > Homo_sapiens.GRCh38.97.transcript.coding.name.chr.refFlat 
$ head -n5 Homo_sapiens.GRCh38.97.transcript.coding.name.chr.refFlat
MT-CO1  MT-CO1-201      chrMT   +       5904    7445    5904    7445    1       5904,   7445,   protein_coding  protein_coding
MT-ND2  MT-ND2-201      chrMT   +       4470    5511    4470    5511    1       4470,   5511,   protein_coding  protein_coding
MT-ND6  MT-ND6-201      chrMT   -       14149   14673   14149   14673   1       14149,  14673,  protein_coding  protein_coding
MT-CO3  MT-CO3-201      chrMT   +       9207    9990    9207    9990    1       9207,   9990,   protein_coding  protein_coding
MT-ND4  MT-ND4-201      chrMT   +       10760   12137   10760   12137   1       10760,  12137,  protein_coding  protein_coding

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

WIndows上でLinuxをエミュレートするWindows Subsystem for Linux (WSL)はDockerに不完全な対応だったのですが、完全対応の「WSL 2」がいつのまにか使えるようになっていたので、試してみました。

forest.watch.impress.co.jp

Windows Insiderに登録

WSL2がリリースされている「Windows Insider Preview Build 18917」の提供はWindows Insider ProgramのFastリング向けなので、通常のWindows Updateではインストールされません。"Windows Insider Program"に登録して、開発版のUpdateをインストールできるようにしましょう。

(注:このリリースは開発者向けはテスト版であるため、Windowsの予期せぬエラーが今後起きる可能性があります。「よくわからない」という方は、正式リリース(10月予定)まで待っておいた方が無難です。)

Windows Insider Program | 最新の Windows 機能を入手する

f:id:rnakato:20190707090010j:plain:w400

登録後、Windowsの「更新とセキュリティ」の画面でWindows Insider Program を選択し、Insiderの設定にFastを選択します。

登録後、あらためてWindows Updateを実行すると"Windows 10 Insider Preview18932"がインストールされました。

f:id:rnakato:20190707091223j:plain:w500

Hyper-Vの有効化

WSL2にはHyper-Vが必要です。Hyper-VWindowsにもともとインストールされていますが、有効化されていません。以下の手順で有効化しましょう。

(注:Hyper-Vを有効化するとVMWare が使えなくなりますので、VMWare を利用している人は注意。)

  • Windows の「アプリと機能」の一番下の関連設定にある「プログラムと機能」のリンクをクリック
  • 左のメニューの「Windows の機能の有効化または無効化」を選択
  • "Hyper-V" のボックスにチェックを入れてOKをクリック

(注2:もしHyper-Vの「Hyper-V プラットフォーム」がグレーアウトして有効化できない場合は、BIOSファームウェアで仮想化サポートが無効になっています。BIOS画面を立ち上げ、仮想化サポートを有効にしましょう。)

参考:

Windows 10 にHyper-Vの機能をインストールする (Windows Tips)

Hyper-Vを有効にする時のBIOSの設定 (HP EliteBook Folio 9470m) - メンチカツには醤油でしょ!!

WSL2のインストール

WSL1が既にインストールされている環境で、Windows Powershellを管理者権限で起動し、以下のコマンドを実行します。

Enable-WindowsOptionalFeature -Online -FeatureName VirtualMachinePlatform

Windowsを再起動した後、あらためてWindows Powershellを管理者権限で起動し、以下のコマンドを実行します。

wsl --set-version Ubuntu-18.04 2

これはUbuntu-18.04の場合のコマンドです。Ubuntu-18.04でないディストリをWSLにインストールしている場合は、適宜変更してください。
(現在インストールされているディストリは wsl cat /etc/issue または wsl -l -v で確認可能です。)

なお、"Windows 10 Insider Preview18932"がインストールされていない状態で上記コマンドを実行すると、「--set-versionというオプションはありません」というエラーになります。

確認

アップデートが完了すると(手元の環境では数分で完了しました)、WSL2が使えるようになります。使用感はWSL1と特に変わりありません。
Dockerをインストールしたところ、問題なく動くようです。Dockerについてはまた日を改めて書く予定です。

(7/9追記:アップデートに必要な時間は、どうもWSL上にどのくらいの量のファイルを置いているかに依存するようです(コピーとかしている?)たくさんファイルがある状態でやると数時間くらいかかりました。)

(8/3追記:WSL2にした時にVcXsrv が使えなくなる問題の解決法: https://github.com/microsoft/WSL/issues/4106
/etc/resolv.conf で表示されるnameserverのIPアドレスをDISPLAY変数に格納すればよいです。)

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

たかが公開鍵、されど公開鍵。「公開鍵 生成」でググるとたくさん出てくるんだけど、やり方が色々ありすぎて新人に適当に検索させるとあれこれ迷ってしまうので、ここにまとめておきます。

参考記事

公開鍵とは

サーバやGitHubなどのリモートサービスにsshでログインする際、通常のパスワード形式を用いると、パスワードが何らかの理由で漏洩した場合に他人にログインされる危険があります。そしてパスワードは常に漏洩の危険があります。
これに対して公開鍵形式の場合は、パスワードに関連付けられた「鍵と錠前」を生成します。パスワードを知っており、かつ鍵を持っている人のみがログイン可能であるため、よりセキュアになります。
この「鍵と錠前」のことを秘密鍵と公開鍵と呼びます。この二つはセットになっており、同時に生成されます。
公開鍵は錠前なので、ログインしたいサーバ上に置きます。公開鍵は他人に見られてもかまいません。
秘密鍵は鍵なので、自分だけが持っておくものです。他人に見られてはいけません。  

公開鍵の作成方法

ここではLinuxでのやり方に統一します。暗号化方式はRSAを用います。Windows10の場合はWSL、Macの場合はMacターミナル上で作業してください。
Windows10でWSLを未インストールの方は以下を参考にインストールしてください。

【WSL入門】第1回 Windows 10標準Linux環境WSLを始めよう:ITの教室 - @IT

既存の公開鍵が存在するかチェック

公開鍵はデフォルトでは~/.sshに生成・保管されます(~/はホームディレクトリ)。 既に作成済みの公開鍵を上書きしないよう、まずは公開鍵を持っているか調べてみましょう。

$ ls ~/.ssh

RSA形式の場合、秘密鍵id_rsa、公開鍵はid_rsa.pubという名前で保存されています。これらのファイルが無いか、または~/.sshディレクトリ自体が存在しない場合は、公開鍵はまだ生成されていません。

公開鍵・秘密鍵の生成

以下のssh-keygenコマンドで公開鍵と秘密鍵が生成されます。 鍵長はデフォルトの2048で良いのですが、ここではよりセキュアな4096を指定することにします。

$ ssh-keygen -t rsa -b 4096
Generating public/private rsa key pair.
Enter file in which to save the key (/home/rnakato/.ssh/id_rsa):
Enter passphrase (empty for no passphrase):
Enter same passphrase again:
Your identification has been saved in /home/rnakato/.ssh/id_rsa.
Your public key has been saved in /home/rnakato/.ssh/id_rsa.pub.
The key fingerprint is:
SHA256:6KfApJsjSs/tj6FPQioHS2MWk4fYlnJXaCp9uowg8Ac rnakato@RyuHomePC
The key's randomart image is:
+---[RSA 4096]----+
|     ..          |
|..o.o.           |
|o*=+.            |
|o+E..  .         |
|oB =. . S        |
|*o*+..           |
|==o+oo. .        |
|=.*o=.oo         |
|o.o=o=o.         |
+----[SHA256]-----+

パスフレーズの入力を求められますが、ここで入力したものがサーバにログインする際のパスワードになります。 サーバ上にあるアカウントのパスワード(sudoなどで利用するもの)とは別であることに注意してください。

パスワードは5文字未満だとエラーになります。しかし空パスワードはOKのようです(謎ですが…)。 当然ですが、空パスワードや容易に推測可能な単純なパスワードは危険ですので、十分複雑なパスワードを指定してください。

終了後再び$ ls ~/.sshを実行し、id_rsaid_rsa.pubが生成されていれば成功です。

公開鍵を使ったログイン

公開鍵の設定自体はサーバサイドで行われるものなので、ユーザの設定は不要です。 サーバ管理者に公開鍵を送付すれば準備OKです。

Linuxシェル上だと、以下のような感じでログインするでしょう。

$ ssh rnakato@hostname

何も指定しなかった場合、秘密鍵~/.ssh/id_rsaが自動的に利用されます。秘密鍵の場所が違ったり、ファイル名が異なる場合は以下のように明示的に指定しましょう。

$ ssh rnakato@hostname -i <秘密鍵>

公開鍵は複数のサーバで共用可能です。その場合すべて同じ秘密鍵とパスワードを用いてログインすることが可能です。 逆に秘密鍵を無くしたりパスワードを忘れて作り直したという場合は、サーバ上の公開鍵もすべて新しいものに更新する必要があります。

秘密鍵パーミッション

秘密鍵は「誰にも見られない」ことが暗黙に想定されているので、秘密鍵パーミッションがオープン過ぎる場合、以下のようなエラーになり、ログインできません。

$ ssh rnakato@hostname
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@         WARNING: UNPROTECTED PRIVATE KEY FILE!          @
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Permissions 0777 for '/home/rnakato/.ssh/id_rsa' are too open.
It is required that your private key files are NOT accessible by others.
This private key will be ignored.
bad permissions: ignore key: /home/rnakato/.ssh/.ssh/id_rsa
Permission denied (publickey).

この場合は以下のようにパーミッションを400に変更すればログインできるようになります。

$ chmod 400 ~/.ssh/id_rsa

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

マップデータの形式にはSAM, BAMの他にCRAMという形式があります。

https://www.ga4gh.org/news/cram-compression-for-genomics/

CRAMはBAMと比べて更に高圧縮率だそうです。
今まであまり使う機会が無かったのですが、Twitterで Ewan Birneyさんが強く推していたので、今日はCRAMを試してみようと思います。

SAMtoolsのインストール

今はSAMtools でCRAMも扱えるんですね。便利ですね。
SAMtoolsのインストールは↓を参照してください。

SAMtoolsとリダイレクト - Palmsonntagmorgen

手元のSAMtoolsのバージョンは1.9です。

$ samtools --version
samtools 1.9
Using htslib 1.9
Copyright (C) 2018 Genome Research Ltd.

ゲノムデータの準備

CRAMの変換にはゲノム配列データが必要になります。ゲノム配列作成は以下を参照してください。

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

また、ゲノムに indexが付加されていることが推奨なので、faidxコマンドでindexを作成しておきます。

$ samtools faidx genome_hg19.fa

genome_hg19.fa.fai が作成されます。

bamのダウンロード

bamファイルは何でもいいのですが、今回は以下からダウンロードした GSM733643_hg19_wgEncodeBroadHistoneK562Pol2bStdAlnRep2.bam を使います。

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM733643

(リンク先下部の"Supplementary file" からBAMがダウンロードできます)

BAM -> CRAMの変換

SAMtoolsを使ってBAMからCRAMへ変換してみます。上でindexを作成しましたが、-Tで指定するのはゲノム配列であることに注意してください。 ここではtime コマンドを追加して実行時間を計測しています。
bamのファイル名が長いと見づらいので、ここでは仮にPol2.bamとします。

$ time samtools view -C Pol2.bam -T genome_hg19.fa > Pol2.cram

real    1m46.909s
user    1m37.133s
sys 0m5.343s

変換に時間がかかるという噂がありましたが、シングルコアでも1分ほどで終わりました。

$ ls -lh
合計 906M
-rwxrwxrwx 1 rnakato rnakato 618M  410 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 288M  410 17:15 Pol2.cram

確かに、ファイルサイズが半分以下になっています。情報量が同じでサイズが半分以下なのは良いですね。

CRAM -> BAMの変換

ファイルが可逆変換か確認するために、CRAMからBAMに変換してみます。 元のBAMファイルと同一のファイルが生成されることを期待しています。

$ time samtools view -b Pol2.cram > Pol2.bam2
real    1m19.342s
user    1m17.821s
sys 0m1.135s

BAMに戻す時にはゲノムを指定する必要はないようです。

$ ls -lh
合計 1.6G
-rwxrwxrwx 1 rnakato rnakato 618M  410 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 683M  410 17:17 Pol2.bam2
-rwxrwxrwx 1 rnakato rnakato 288M  410 17:15 Pol2.cram

何故だか元bamファイルよりファイルサイズが大きくなっています。

3つともsamファイルに変換して調べたところ、Pol2.bam2はPol2.bamと比べてカラムが2つ増えていました。詳細をきちんと確認していませんが、このあたりは元SAM・BAMファイルをどのツールで生成したかによって結果が変わるかもしれません。 なお、Pol2.bam2とPol2.cramは同一でした。

CRAMのソート

$ samtools sort Pol2.cram -O cram -@ 12 > Pol2.sort.cram

問題なくソートできました。 -O cramで出力をCRAMに指定するのを忘れるとbamになってしまうので注意。

bowtie2からのリダイレクト

最後に、bowtie2でマッピングした結果をリダイレクトして直接ソート済CRAMを生成することにチャレンジ。

まず元fastqのダウンロード(何でもよいです)

$ fastq-dump SRR227359

bowtie2でマッピングしてみます。bowtie2のデフォルト出力はSAM形式です。

# SAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq > SRR227359.sam
# sorted BAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort > SRR227359.sort.bam
# sorted CRAM
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq | samtools sort -O cram - -T genome_hg19.fa > SRR227359.sort.cram
[bam_sort_core] merging from 5 files and 1 in-memory blocks...
[E::cram_get_ref] Failed to populate reference for id 0
samtools sort: failed to write header to "-"

CRAMはsamtools sortに直接投げるとエラーになりました。 修正して以下のように2段階のリダイレクトにします。

# sorted CRAM(修正)
$ bowtie2 -x bowtie2-indexes/UCSC-hg19 -p 24 SRR227359.fastq \
  | samtools view -C - -T genome_hg19.fa \
  | samtools sort -O cram \
  >  SRR227359.sort.cram

今度はうまく行きました。

BAMの時と比較しても実行時間はそれほど変わりません。計算時間的な意味でのオーバーヘッドはほとんど気になりません。

考察

SAMtoolsを使うとCRAMはほとんどストレスなく取り扱えます。ゲノム配列を指定するのがちょっと面倒なくらいです。 計算時間はほとんど変わらず、ファイルサイズは小さくなるので、有用性は高いです。
あとは各種解析ツールが入力にCRAM形式をどの程度受け付けるかによりますね。今のところは対応しているツールはそれほど多くないと思うので、「しばらく使う予定は無いが消すことはできない」というような大量のBAMファイルがある場合にCRAM形式にしておくと容量が削減できてうれしい、というところでしょうか。

あとはpaired-endの時にどうなるのかなどチェックした方がいいかもしれません。