Palmsonntagmorgen

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

SSHその3　ssh の設定を保存 (.ssh/config)

Linux 解析環境構築

SSHの記事その３です。ここまでの流れは以下の記事を参照してください。

rnakato.hatenablog.jp

rnakato.hatenablog.jp

「その２」の記事で、 .ssh/configに以下を記入することでサーバ間のssh, scpをする際に秘密鍵、パスワードの入力が不要になると述べました。

ForwardAgent yes

.ssh/は秘密鍵など、sshのためのファイルが格納されている場所です。ここにconfigという名前のファイルを作成すると、その中に記載されている設定がデフォルトで適用されるようになり、sshコマンドを打つ時にオプションを指定する必要がなくなります。

以下では .ssh/configファイルの利用例を紹介します。

接続切れを防ぐ

ServerAliveInterval 100

sshはサーバから応答がないと自動で接続が切れてしまいますが、このように記載しておくと接続が切れるまでの時間を延ばすことができます。

参考：ssh接続が自動切断される場合の回避方法 - セキュリティ

接続先でXウィンドウを使う

接続先のサーバでGUI（ブラウザなど）を起動したい場合は通常sshに -Xオプションを追加しますが、configに以下のように記載することでオプションの指定が不要になります。

X11Forwarding yes

参考：X11 forwarding — 京大マイコンクラブ (KMC)

秘密鍵の場所を指定

秘密鍵がデフォルトのパス（ ~/.ssh/id_rsaなど）でない場所にある場合、configに以下のように記載することで秘密鍵の場所を指定することができます。

IdentityFile <PATH_to_id_rsa>

サーバに接続するポートを指定

サーバに接続する際の設定をサーバ別に指定することもできます。以下はserver_aという名前のサーバに接続する時にデフォルトのポートを12345にするという設定です。

Host server_a
  User <ユーザ名>
  Hostname server_a
  Port 12345
  IdentityFile <PATH_to_id_rsa>

IdentityFile を指定することで、「このサーバに接続する時のみこの秘密鍵を使う」などのような指定方法も可能です。

多段sshログイン

サーバが複数台あり、特定のサーバ（ここではserver_humidaiとします）を踏み台としてその他のサーバ（ここではserver_aとします）に接続する場合、ローカルPC→server_humidai→server_aとsshしなければなりません。scpなどを使う時は面倒ですね。

ここでconfigに以下のように記載すると、直接 ssh server_aとすることでログイン可能になります。（ForwardAgent yesによって公開鍵パスワードの入力を省略しています）

ForwardAgent yes

Host server_humidai
  User <ユーザ名>
  Hostname server_humidai
Host server_a
  User <ユーザ名>
  Hostname server_a
  ProxyCommand ssh -W %h:%p server_humidai

まとめ

sshは毎日のように使うコマンドですので、configにできる限り設定を書いておいてコマンド実行を楽にするのがよいと思います。自分なりのconfigファイルをぜひ作成してみてください。

SSHその2　ssh のパスワード入力を省略 (ssh-agent)

Linux 解析環境構築

諸事情によりまたしばらく基礎的な内容になります。

クラウドサービス含め、サーバに接続する時にはセキュリティ面を考えて公開鍵方式のSSHを利用するのが一般的です。 公開鍵については過去に記事にしています。

rnakato.hatenablog.jp

しかしサーバに接続するたびにパスワードを入力するのは面倒ですね。だからと言って、短く簡単なパスワードにしてしまう（あるいはパスワードなしにしてしまう）のはセキュリティ面で大きな問題になります。

そんな時にはssh-agentを使うと、パスワードの入力を省略することができるようになります。

参考記事

Linuxの場合

Linux (Windows WSL含む) では、以下のコマンドで起動します。

$ eval `ssh-agent`

次に、以下のコマンドで公開鍵のパスワードをssh-agentに登録します。

$ ssh-add  <秘密鍵のファイル名>
(秘密鍵のパスワード入力)

こうすると、以後サーバログイン時に公開鍵のパスワード入力が不要になります。

ただこの方法では、ターミナルを立ち上げるたびに自分でssh-agentを起動しなければなりません。 ターミナル起動時（またはログイン時）にssh-agentが自動起動されるようにするには、ホームディレクトリにある .bashrc (または.bash_profile)に以下の内容を追記します。

if [ -f ~/.ssh-agent ]; then
    . ~/.ssh-agent
fi
if [ -z "$SSH_AGENT_PID" ] || ! kill -0 $SSH_AGENT_PID; then
    ssh-agent > ~/.ssh-agent
    . ~/.ssh-agent
fi
ssh-add -l >& /dev/null || ssh-add

ここでの記述方法には色々あるのですが、ひとまずこのままコピー&ペーストしてしまって大丈夫です。

Macの場合

上の参考記事によれば、Mac OS X Leopard 〜 El Capitan まではssh-agentのかわりにKeychainというプログラムが自動的に実行されており、ssh-agentのマニュアル起動は不要でしたが、macOS Sierra では自動では実行されない設定になったそうです。 自動起動のため、Sierra では以下の設定を行います。

（実機で確認していないため、訂正あればコメントをお願いします。）

  • ~/.ssh/config に以下を追記（ファイルがなければ新規作成）
Host *
  AddKeysToAgent yes
  UseKeychain yes
  IdentityFile <秘密鍵のパス>

なお最後の行 IdentityFile <秘密鍵のパス> は秘密鍵の場所を毎回 -iオプションで指定しなくてもよい設定にするためのもので、ssh-agentとは直接関係はありません。

サーバ上にも設定を反映する：ForwardAgent

これはMac/Linux共通です。

使っているサーバが複数ある場合、例えばローカルPCからサーバAにssh接続し、サーバAからサーバBにsshでログインする、あるいはscpでデータを送りたいというケースが発生します。 そうなるとサーバ上にも秘密鍵を置く必要がありますが、それはセキュリティ的に問題です。 ここでForwardAgentという仕組みを使うと、ローカルPCのssh設定をそのままログインしたサーバに反映させることができます。

具体的には ~/.ssh/config に以下を追記するだけでOKです。（ファイルがなければ新規作成）

ForwardAgent yes

こうするとローカルのssh-agent設定がサーバに引き継がれますので、サーバ間のssh, scpをする際に秘密鍵、パスワードの入力が不要になります。

注意点

ssh-agentを起動している間はパスなしでサーバにログインできてしまうため、共用PCなどで設定する際は注意が必要です。 他の人が使う可能性がある環境、職場で席を長く離れる場合などは ssh-add -D で登録した秘密鍵を削除するなどするとより安心です。

論文・グラント申請書作成のポイントリンク集

その他 Tips

随時追加します。Tweetはクリックするとスレッドが読めます。

Academic writing

アカデミックマナーの心得 - 東京大学大学院 情報学環・学際情報学府

Writing in the Sciences | Coursera

Grant writing

note.com

https://www.pnas.org/doi/10.1073/pnas.2315735121

jabberwocky.weecology.org

英語論文作成ツール集

ジャーナルに論文を通すまで

  • Natureに論文出版した人の記録。

toshikisato.wixsite.com

エピゲノム論文にありがちな間違いについて

Conference abstract

Recommendation letter

iwasakiichiro.info

fetchChromSizes を使ってgenome tableファイルを作成

NGS Genome

以前書いた以下のエントリでは、genome配列ファイルからgenome tableファイルを作成する方法を紹介しました。

rnakato.hatenablog.jp

UCSCから提供されているfetchChromSizes を使うと、この作業をより簡便に行うことが可能です。

fetchChromSizesのインストール

UCSCサイトから fetchChromSizes を直接ダウンロード可能です。

wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/fetchChromSizes
chmod +x fetchChromSizes

またはbiocondaでもインストール可能です。

conda install -c bioconda ucsc-fetchchromsizes

fetchChromSizesの利用

例えばhg38のgenometableを作成するには以下のようにします。

$ fetchChromSizes hg38 > genometable.hg38.txt
INFO: trying CURL  for database hg38
url: http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes
$ head genometable.hg38.txt
chr1    248956422
chr2    242193529
chr3    198295559
chr4    190214555
chr5    181538259
chr6    170805979
chr7    159345973
chrX    156040895
chr8    145138636
chr9    138394717

ログに表示されている通り、このコマンドはUCSC database で公開されている chrom.sizesファイルhttp://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes　をダウンロードしているだけです。なので、以下のコマンドと本質的に同じことをしています。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes -O genometable.hg38.txt
--2021-05-24 20:15:40--  http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes
...
2021-05-24 20:15:40 (163 MB/s) - genometable.hg38.txt へ保存完了

$ head genometable.hg38.txt
chr1    248956422
chr2    242193529
chr3    198295559
chr4    190214555
chr5    181538259
chr6    170805979
chr7    159345973
chrX    156040895
chr8    145138636
chr9    138394717

URLを直接指定しなくてよいぶん、fetchChromSizesの方が簡便ですね。

fetchChromSizesの応用

fetchChromSizesは結果を標準出力に表示しますので、出力をパイプすることができます。

上記のコマンドはコンティグ配列などを全て含んだgenometableを生成しますが、以下のコマンドでコンティグ配列分を除去したgenometableを作成できます。

# grep -vを使う場合
$ fetchChromSizes hg38 | grep -v _ > genometable.hg38.txt
# awkを使う場合（より厳密）
$ fetchChromSizes hg38 | awk -vOFS="\t" '($1!~/_/)' > genometable.hg38.txt
$ cat genometable.hg38.txt
chr1    248956422
chr2    242193529
chr3    198295559
chr4    190214555
chr5    181538259
chr6    170805979
chr7    159345973
chrX    156040895
chr8    145138636
chr9    138394717
chr11   135086622
chr10   133797422
chr12   133275309
chr13   114364328
chr14   107043718
chr15   101991189
chr16   90338345
chr17   83257441
chr18   80373285
chr20   64444167
chr19   58617616
chrY    57227415
chr22   50818468
chr21   46709983
chrM    16569

以下のようにするとbed形式にできます。

$ fetchChromSizes hg38 \
  | awk -vOFS="\t" '($1!~/_/){ print $1, "0", $2 }' \
  > hg38.bed
$ cat hg38.bed
chr1    0       248956422
chr2    0       242193529
chr3    0       198295559
chr4    0       190214555
chr5    0       181538259
chr6    0       170805979
chr7    0       159345973
chrX    0       156040895
chr8    0       145138636
chr9    0       138394717
chr11   0       135086622
chr10   0       133797422
chr12   0       133275309
chr13   0       114364328
chr14   0       107043718
chr15   0       101991189
chr16   0       90338345
chr17   0       83257441
chr18   0       80373285
chr20   0       64444167
chr19   0       58617616
chrY    0       57227415
chr22   0       50818468
chr21   0       46709983
chrM    0       16569

結果をソートしたい場合は以下のようにしましょう。

$ fetchChromSizes hg38 \
   | awk -vOFS="\t" '($1!~/_/){ print $1, "0", $2 }' \
   | sort -k1,1 -k2,2n \
   > hg38.bed
$ cat hg38.bed
chr1    0       248956422
chr10   0       133797422
chr11   0       135086622
chr12   0       133275309
chr13   0       114364328
chr14   0       107043718
chr15   0       101991189
chr16   0       90338345
chr17   0       83257441
chr18   0       80373285
chr19   0       58617616
chr2    0       242193529
chr20   0       64444167
chr21   0       46709983
chr22   0       50818468
chr3    0       198295559
chr4    0       190214555
chr5    0       181538259
chr6    0       170805979
chr7    0       159345973
chr8    0       145138636
chr9    0       138394717
chrM    0       16569
chrX    0       156040895
chrY    0       57227415

【DROMPAplus】Dockerを使ったparse2wig+

Docker ChIP-seq DROMPA NGS

そういえばDROMPAplusの紹介記事を全然書いていないことに気づきました。 最近はDockerイメージからのDROMPAplusの使い方をたびたび質問いただくので、今日はDockerを使ったDROMPAplusの使い方を紹介します。

なお、旧バージョンのDROMPA3とDROMPAplusの違いについては以下を参照してください。

rnakato.hatenablog.jp

本記事でやること

以下の記事ではDROMPA3に含まれるparse2wigの使い方を紹介しました。

rnakato.hatenablog.jp

この記事を docker parse2wig+ に切り替えてみましょう。 なお、Dockerそのもののインストールについては↓の記事を参照してください。

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

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

rnakato/ssp_drompa という名前でDockerイメージを提供しています。 以下のコマンドでイメージをローカルにダウンロードしてください。

 $ docker pull rnakato/ssp_drompa

ちなみに本イメージにはDROMPAplusのほか、SSPとDROMPA3も含まれています。

 $ docker run --rm -it rnakato/ssp_drompa ls
DROMPA3  DROMPAplus  SSP

BAMファイルダウンロード

BAMファイルのダウンロードは特に変更はありません。

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwTfbs/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam

Genome table作成

自分でGenome tableを作成する方法は以下の記事を参照してください。

genome tableを作成する - Palmsonntagmorgen

本Dockerイメージでは、/home/DROMPAplus/data/genometable/ ディレクトリに代表的なGenome tableファイルが含まれています。

 $ docker run --rm -it rnakato/ssp_drompa ls /home/DROMPAplus/data/genometable/
genometable.dm6.txt   genometable.hg38.txt  genometable.sacCer3.txt
genometable.hg19.txt  genometable.mm10.txt

実行

ローカルの場合

ダウンロードしたBAMファイルをparse2wig+でパースする場合、ローカルにインストールしたDROMPAplusなら以下のコマンドで実行可能です。

 $ parse2wig+ -i wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf --gt genometable.txt 
======================================
parse2wig+ version 1.8.8

Input file: wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
Output prefix: parse2wigdir+/Ctcf
Output format: BIGWIG
Binsize: 100 bp
Genome-table file: /work/Database/UCSC/hg19/genome_table
Single-end mode: fragment length will be estimated by strand-shift profile
PCR bias filtering: ON
        10000000 reads used for library complexity

Total read normalization: NONE
Number of threads: 1
======================================
Parsing wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam...

DROMPA3の時と比較して、--gt がハイフン２つになっていることに注意してください。

Dockerの場合

一方Dockerでは同じコマンドを打ち込むとエラーになります。

 $ docker run --rm -it rnakato/ssp_drompa parse2wig+ -i wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf --gt genometable.txt 
Error: Could not open genometable.txt.

genometable.txtが見つかりませんと出ます。dockerコマンドはカレントディレクトリを自動でマウントしないため、オプションで明示的にマウントしなければコマンドからgenometable.txtと入力BAMファイルが見えません。

マウントするには -vオプションを追加します。ここでは例として、カレントディレクトリをコンテナ内の /mnt ディレクトリにマウントします。（長くなるので２行に分割します）

 $ docker run --rm -it -v $(pwd):/mnt rnakato/ssp_drompa parse2wig+ \
   -i /mnt/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf --gt /mnt/genometable.txt 

======================================
parse2wig+ version 1.8.8

Input file: /mnt/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
Output prefix: parse2wigdir+/Ctcf
Output format: BIGWIG
Binsize: 100 bp
Genome-table file: /mnt/genometable.txt
Single-end mode: fragment length will be estimated by strand-shift profile
PCR bias filtering: ON
        10000000 reads used for library complexity

Total read normalization: NONE
Number of threads: 1
======================================
Parsing /mnt/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam...

$(pwd) がカレントディレクトリを示しています。他にも、/Users をマウントしたければ -v /Users:/mnt でOKです。

コンテナ内の/mntディレクトリにマウントしているので、ホスト側のカレントディレクトリにあるファイルは/mnt以下に存在するように見えます。よって、コマンドの入力BAMファイルとgenometable.txtの前に/mnt/をつける必要があります。

出力ディレクトリオプションをつける

実は上記のコマンドにはまだ問題があります。parse2wig+はデフォルトでカレントディレクトリに parse2wigdir+を作成し、その中に結果を出力しますが、コンテナ内のカレントディレクトリはホスト側から見えません。このままでは出力結果を手に入れることができません。

そこでここでは --odirオプションを追加して /mnt/parse2wigdir+ に結果が出力されるようにします。

 $ docker run --rm -it -v $(pwd):/mnt rnakato/ssp_drompa parse2wig+ \
   -i /mnt/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf --gt /mnt/genometable.txt \
   --odir /mnt/parse2wigdir+
... (コマンド完了)
 $ ls parse2wigdir+
Ctcf.100.bw  Ctcf.100.tsv

無事ファイルが生成されました。

おまけ：Docker内genometableを利用する

上で述べたように、本Dockerイメージにはgenometable.txtが含まれています。これを使えば別途genometableを用意することなくparse2wig+を実行できます。

 $ docker run --rm -it -v $(pwd):/mnt rnakato/ssp_drompa parse2wig+ \
   -i /mnt/wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam -o Ctcf --odir /mnt/parse2wigdir+ \
   --gt /home/DROMPAplus/data/genometable/genometable.hg19.txt
... (コマンド完了)
 $ ls parse2wigdir+
Ctcf.100.bw  Ctcf.100.tsv

まとめ

Dockerからparse2wig+を起動する方法を紹介しました。 Dockerからコマンドを実行する場合、「カレントディレクトリのファイルが見えない」という点で困ってしまう方が多いようです。-vコマンドを使ったマウントは最初はホスト側とコンテナ側のディレクトリの関係を考えるのがややこしいですが、慣れてしまえば平気になります。

なお、dockerコマンドを使って生成されたファイルは所有者権限がrootになってしまうという問題があります。これを回避する方法を次回紹介します。

２つのファイルの共通行を抽出する

NGS Linux

前回に続き、今回も小ネタです。

２つの遺伝子リストファイルがあって、２ファイル間で共通する遺伝子名がいくつあるか調べたいとします。

$ ls
genelist1.txt  genelist2.txt
$ cat genelist1.txt
RNA5SP516
AC092299.1
EML3
ZNF492
AC093063.1
DDX17
SRIP1
MIR3167
IGHD3-16
NRDC
ASPHD2
BNIP3L
KRT18P60
$ cat genelist2.txt
MED21
ETS2
MTCO3P39
AC008739.5
RNU1-101P
AC073548.1
MGST2
ASPHD2
BNIP3L
KRT18P60
RNA5SP516
AC092299.1
EML3

grepコマンドを使う方法

こういうファイルの場合はgrepを使うと一発でカウントできます。 ここでは２ファイルで完全一致する行を抽出しています。 （-iオプションを追加すると大文字・小文字を区別しなくなります。）

$ grep -x -f genelist1.txt genelist2.txt
ASPHD2
BNIP3L
KRT18P60
RNA5SP516
AC092299.1
EML3

遺伝子数をカウントしたい時はwcをつけましょう。

$ grep -x -f genelist1.txt genelist2.txt | wc -l
6

参考：

２つのファイルの共通行を抽出する方法 - Qiita

特定の列だけを対象にして抽出したい時

遺伝子名以外の情報も同じ行に含まれていて、遺伝子名だけで完全一致をカウントしたい時、上記のgrepは使えません。

$ cat genelist1.txt
RNA5SP516       ENSG00000201000 chrX    +       142411052       142411171       142411052       142411171       1       142411052,      142411171,      rRNA_pseudogene rRNA_p
                RNA5SP516-201   ENST00000364130
AC092299.1      ENSG00000282416 chr19   -       197309  198066  197309  198066  1       197309, 198066, processed_pseudogene    processed_pseudogene    AC092299.1-201  ENST00000632679
EML3    ENSG00000149499 chr11   -       62602217        62612775        62602474        62612457        22      62602217,62602758,62603148,62603728,62603943,62604112,62605112,62605641,62605854,62606062,62606957,62607665,62608200,62608541,62608735,62608961,62609353,62609628,62610878,62611086,62611424,62612435,        62602678,62602889,62603247,62603816,62604041,62604201,62605180,62605773,62605980,62606214,62607099,62607821,62608296,62608652,62608805,62609132,62609477,62609696,62610992,62611344,62611596,62612775,     p
                protein_coding  EML3-202        ENST00000394773
ZNF492  ENSG00000229676 chr19   +       22634323        22667671        22653399        22665265        4       22634323,22653306,22653919,22663799,    22634474,22653433,2265
                protein_coding  protein_coding  ZNF492-201      ENST00000456783
AC093063.1      ENSG00000268088 chr19   +       39693924        39696258        39693924        39696258        3       39693924,39694463,39696142,     39693960,39694674,3969
        transcribed_unprocessed_pseudogene      transcribed_unprocessed_pseudogene      AC093063.1-201  ENST00000599317
DDX17   ENSG00000100201 chr22   -       38483439        38506337        38485934        38506237        13      38483439,38487872,38492055,38493709,38494020,38494629,38494885
                                                        38486440,38488115,38492115,38493771,38494131,38494802,38495046,38495937,38498150,38498573,38499499,38501280,38506337, p               protein_coding  DDX17-202       ENST00000396821
SRIP1   ENSG00000248660 chr4    +       58103146        58103434        58103146        58103434        1       58103146,       58103434,       processed_pseudogene    proces
                SRIP1-201       ENST00000489235
MIR3167 ENSG00000266215 chr11   -       126988457       126988542       126988457       126988542       1       126988457,      126988542,      miRNA   miRNA   MIR3167-201  ENST00000579623
IGHD3-16        ENSG00000211917 chr14   -       105895633       105895670       105895633       105895670       1       105895633,      105895670,      IG_D_gene       IG_D_g
        IGHD3-16-201    ENST00000390577
NRDC    ENSG00000078618 chr1    -       51789190        51878805        51789235        51878615        33      51789190,51789567,51790532,51790899,51791577,51792045,51792376,51794471,51794822,51798248,51800555,51803813,51805509,51806793,51809314,51810280,51811993,51814034,51814550,51814692,51816311,51818065,51819799,51821497,51823663,51825286,51
                                                                51789433,51789657,51790649,51790990,51791661,51792098,51792424,51794610,51794854,51798411,51800683,51803964,51805561,51806913,51809401,51810404,51812098,51814089,51814609,51814813,51816389,51818135,51819873,51821555,51823786,51825382,51827869,51834170,51836212,51836438,51837620,51840
                protein_coding  protein_coding  NRDC-202        ENST00000354831
ASPHD2  ENSG00000128203 chr22   +       26429259        26445015        26433615        26443206        4       26429259,26433391,26442458,26443096,    26429486,26434501,2644
                protein_coding  protein_coding  ASPHD2-201      ENST00000215906
BNIP3L  ENSG00000104765 chr8    +       26383196        26505636        26383196        26505319        6       26383196,26391242,26395229,26407999,26408226,26505306,  263832
                                                        protein_coding  protein_coding  BNIP3L-207      ENST00000523949
KRT18P60        ENSG00000215208 chr12   -       65418730        65420027        65418730        65420027        1       65418730,       65420027,       processed_pseudogene p
                        KRT18P60-201    ENST00000312876
$ cat genelist2.txt
MED21   ENSG00000152944 chr12   +       27022557        27030673        27022579        27028461        4       27022557,27026419,27027346,27028284,    27022621,27026534,2702
                protein_coding  protein_coding  MED21-201       ENST00000282892
ETS2    ENSG00000157557 chr21   +       38805182        38824955        38810034        38822889        11      38805182,38805970,38810034,38813002,38814272,38814780,38817007
                                        38805616,38806120,38810106,38813114,38814392,38814981,38817091,38818646,38819766,38821704,38824955,     protein_coding  protein_coding
        ETS2-201        ENST00000360214
MTCO3P39        ENSG00000268736 chr4    +       49246020        49246538        49246020        49246538        1       49246020,       49246538,       processed_pseudogene p
                        MTCO3P39-201    ENST00000514812
AC008739.5      ENSG00000279404 chr19   -       20926776        20928714        20926776        20928714        1       20926776,       20928714,       TEC     TEC     AC0087
                ENST00000623035
RNU1-101P       ENSG00000212473 chr8    -       70077905        70078032        70077905        70078032        1       70077905,       70078032,       snRNA   snRNA   RNU1-1
        ENST00000391171
AC073548.1      ENSG00000279861 chr19   -       47487637        47489519        47487637        47489519        1       47487637,       47489519,       TEC     TEC     AC0735
                ENST00000624088
MGST2   ENSG00000085871 chr4    +       139665767       139740745       139666019       139704148       6       139665767,139678542,139695196,139703454,139704015,139740211, 1
                                                                protein_coding  protein_coding  MGST2-207       ENST00000616265
ASPHD2  ENSG00000128203 chr22   +       26429259        26445015        26433615        26443206        4       26429259,26433391,26442458,26443096,    26429486,26434501,2644
                protein_coding  protein_coding  ASPHD2-201      ENST00000215906
BNIP3L  ENSG00000104765 chr8    +       26383196        26505636        26383196        26505319        6       26383196,26391242,26395229,26407999,26408226,26505306,  263832
                                                        protein_coding  protein_coding  BNIP3L-207      ENST00000523949
KRT18P60        ENSG00000215208 chr12   -       65418730        65420027        65418730        65420027        1       65418730,       65420027,       processed_pseudogene p
                        KRT18P60-201    ENST00000312876
RNA5SP516       ENSG00000201000 chrX    +       142411052       142411171       142411052       142411171       1       142411052,      142411171,      rRNA_pseudogene rRNA_p
                RNA5SP516-201   ENST00000364130
AC092299.1      ENSG00000282416 chr19   -       197309  198066  197309  198066  1       197309, 198066, processed_pseudogene    processed_pseudogene    AC092299.1-201  ENST00000632679
EML3    ENSG00000149499 chr11   -       62602217        62612775        62602474        62612457        22      62602217,62602758,62603148,62603728,62603943,62604112,62605112,62605641,62605854,62606062,62606957,62607665,62608200,62608541,62608735,62608961,62609353,62609628,62610878,62611086,62611424,62612435,        62602678,62602889,62603247,62603816,62604041,62604201,62605180,62605773,62605980,62606214,62607099,62607821,62608296,62608652,62608805,62609132,62609477,62609696,62610992,62611344,62611596,62612775,

cutコマンドで特定列だけ抽出してもよいのですが、中間ファイルができてしまいますし、他の列の情報を以後の解析で用いたい場合は使えません。

ここでは私のscriptに含まれている combine_lines_from2files.pl を使ってみましょう。

スクリプトをダウンロードする時は以下のコマンドを実行してください（要gitコマンド）。

$ git clone https://github.com/rnakato/script_rnakato.git

スクリプトを実行します。-1, -2でファイル名を指定し、-a, -bで対象にする列数（ここでは0列目）を指定します。遺伝子名が冒頭にないファイルでも大丈夫です。

$ combine_lines_from2files.pl -1 genelist1.txt -2 genelist2.txt -a 0 -b 0
ASPHD2  ENSG00000128203 chr22   +       26429259        26445015        26433615        26443206        4       26429259,26433391,26442458,26443096,    26429486,26434501,2644
                protein_coding  protein_coding  ASPHD2-201      ENST00000215906 ASPHD2  ENSG00000128203 chr22   +       26429259        26445015        26433615        264432
        4       26429259,26433391,26442458,26443096,    26429486,26434501,26442572,26445015,    protein_coding  protein_coding  ASPHD2-201      ENST00000215906
BNIP3L  ENSG00000104765 chr8    +       26383196        26505636        26383196        26505319        6       26383196,26391242,26395229,26407999,26408226,26505306,  263832
                                                        protein_coding  protein_coding  BNIP3L-207      ENST00000523949 BNIP3L  ENSG00000104765 chr8    +       26383196     2
        26383196        26505319        6       26383196,26391242,26395229,26407999,26408226,26505306,  26383230,26391426,26395302,26408103,26408376,26505636,  protein_coding
        protein_coding  BNIP3L-207      ENST00000523949
KRT18P60        ENSG00000215208 chr12   -       65418730        65420027        65418730        65420027        1       65418730,       65420027,       processed_pseudogene p
                        KRT18P60-201    ENST00000312876 KRT18P60        ENSG00000215208 chr12   -       65418730        65420027        65418730        65420027        1    6
                65420027,       processed_pseudogene    processed_pseudogene    KRT18P60-201    ENST00000312876
RNA5SP516       ENSG00000201000 chrX    +       142411052       142411171       142411052       142411171       1       142411052,      142411171,      rRNA_pseudogene rRNA_p
                RNA5SP516-201   ENST00000364130 RNA5SP516       ENSG00000201000 chrX    +       142411052       142411171       142411052       142411171       1       142411
        142411171,      rRNA_pseudogene rRNA_pseudogene RNA5SP516-201   ENST00000364130
AC092299.1      ENSG00000282416 chr19   -       197309  198066  197309  198066  1       197309, 198066, processed_pseudogene    processed_pseudogene    AC092299.1-201  ENST00
                AC092299.1      ENSG00000282416 chr19   -       197309  198066  197309  198066  1       197309, 198066, processed_pseudogene    processed_pseudogene    AC0922
                ENST00000632679
EML3    ENSG00000149499 chr11   -       62602217        62612775        62602474        62612457        22      62602217,62602758,62603148,62603728,62603943,62604112,62605112,62605641,62605854,62606062,62606957,62607665,62608200,62608541,62608735,62608961,62609353,62609628,62610878,62611086,62611424,62612435,        62602678,62602889,62603247,62603816,62604041,62604201,62605180,62605773,62605980,62606214,62607099,62607821,62608296,62608652,62608805,62609132,62609477,62609696,62610992,62611344,62611596,62612775,     p
                protein_coding  EML3-202        ENST00000394773 EML3    ENSG00000149499 chr11   -       62602217        62612775        62602474        62612457        22   62602217,62602758,62603148,62603728,62603943,62604112,62605112,62605641,62605854,62606062,62606957,62607665,62608200,62608541,62608735,62608961,62609353,62609628,62610878,6261
                        62602678,62602889,62603247,62603816,62604041,62604201,62605180,62605773,62605980,62606214,62607099,62607821,62608296,62608652,62608805,62609132,62609477,62609696,62610992,62611344,62611596,62612775,        protein_coding  protein_coding  EML3-202        ENST00000394773

出力されるファイルは指定された２つのファイルをマージしたものになります。 中で何をやっているか知りたい人はスクリプトを開いてみてください（適当）。

例によって、行数だけカウントしたい場合はwcをつけます。

$ combine_lines_from2files.pl -1 genelist1.txt -2 genelist2.txt -a 0 -b 0 | wc -l
6

余談

このスクリプトは.plとついている通りPerlスクリプトです。昔に書いたスクリプトで、特に修正する必要もなかったのでそのまま使っていますが、Pythonのpandasなど使えばもっと簡単にできます。

コマンドライン上でgrepのように一発でできる簡単な方法がもしありましたら教えてください。

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

Linux NGS

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を使っています。