Palmsonntagmorgen

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

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

RNA-seq Genome NGS

みなさん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 解析環境構築

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

forest.watch.impress.co.jp

qiita.com

Windows Insiderに登録

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

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

The Windows Insider Program

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

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

Hyper-Vの有効化

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

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

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

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

参考:

Windows 10 にHyper-Vの機能をインストールする : Windows | iPentec

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

Windows 10 Home で WSL 2 + Docker を使う - Qiita

WSL2のインストール

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

Enable-WindowsOptionalFeature -Online -FeatureName VirtualMachinePlatform

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

(2020/11/19追記：Windows Powershellは、ISE (x86)だと "wslというコマンドはない"というエラーになりました。ISEの方では問題なく実行できました。）

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公開鍵の生成・設定の方法

Linux 解析環境構築

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

参考記事

公開鍵とは

サーバや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_rsaとid_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形式を試す

Genome Linux NGS

マップデータの形式には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

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

$ ls -lh
合計 906M
-rwxrwxrwx 1 rnakato rnakato 618M  4月 10 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 288M  4月 10 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  4月 10 17:13 Pol2.bam
-rwxrwxrwx 1 rnakato rnakato 683M  4月 10 17:17 Pol2.bam2
-rwxrwxrwx 1 rnakato rnakato 288M  4月 10 17:15 Pol2.cram

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

３つとも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に直接投げるとエラーになりました。 修正して以下のように２段階のリダイレクトにします。

# 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の時にどうなるのかなどチェックした方がいいかもしれません。

【ご報告】PIになりました

この度4月1日付で東京大学定量生命科学研究所の講師となりました。
大規模生命情報解析研究分野という名前で研究室を主宰します。
研究室HP: http://www.iam.u-tokyo.ac.jp/nakatolab/index.html

年度末営業と研究室異動のばたばたでこのブログも全く更新できていませんでしたが、これからも引き続きゆるく続けていきたいと思っています。
なお当研究室では学生を募集していますので、こういったNGS解析に興味がありそうな学生さんがいましたら是非ご連絡くださいませ。

今後とも宜しくお願い致します。

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

RNA-seq NGS Genome

前回の記事↓の続きです。

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

前回のコマンドを最後まで完了すると、starディレクトリの中にstar/Myers_HUVEC_cell_2x75_200_1.<genes|isoforms>.results のようなファイルがサンプルごとに生成されているのがわかります。 *.genes.resultsと*.isoforms.resultsはそれぞれ遺伝子単位、transcript単位の発現量（count, TPM, FPKM）が出力されます。 transcript単位の出力は基本的には精度が十分でないことが多く、ノイズがまぎれこみやすいので、特に利用目的がないのであれば遺伝子単位のものを利用しましょう。

RSEMの出力

では、ファイルの中をのぞいてみましょう。

$ head star/Myers_HUVEC_cell_2x75_200_1.genes.results
gene_id transcript_id(s)    length  effective_length    expected_count  TPM FPKM
ENSG00000000003 ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2146.69 1926.19 1086.00 40.20   34.28
ENSG00000000005 ENST00000373031,ENST00000485971 1046.60 826.17  14.00   1.21    1.03
ENSG00000000419 ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1003.23 782.75  355.00  32.34   27.57
ENSG00000000457 ENST00000367770,ENST00000367771,ENST00000367772,ENST00000423670,ENST00000470238 3672.40 3451.91 112.23  2.32    1.98
ENSG00000000460 ENST00000286031,ENST00000359326,ENST00000413811,ENST00000459772,ENST00000466580,ENST00000472795,ENST00000481744,ENST00000496973,ENST00000498289 2167.45 1946.97 501.77  18.38   15.67
ENSG00000000938 ENST00000374003,ENST00000374004,ENST00000374005,ENST00000399173,ENST00000457296,ENST00000468038,ENST00000475472 2021.00 1800.50 9.00    0.36    0.30
ENSG00000000971 ENST00000359637,ENST00000367429,ENST00000466229,ENST00000470918,ENST00000496761,ENST00000630130 2587.83 2367.38 0.00    0.00    0.00
ENSG00000001036 ENST00000002165,ENST00000367585,ENST00000451668 2285.18 2064.69 1079.67 37.29   31.79
ENSG00000001084 ENST00000229416,ENST00000504353,ENST00000504525,ENST00000505197,ENST00000505294,ENST00000509541,ENST00000510837,ENST00000513939,ENST00000514004,ENST00000514373,ENST00000514933,ENST00000515580,ENST00000616923,ENST00000643939,ENST00000650454 2327.40 2106.94 773.00  26.16   22.30

Ensemblからダウンロードしたgtfファイルを用いているので、gene_idとtranscript_idは当然Ensemblのものになっています。 effective length はpoly(A) tailを除いた転写物の長さだそうです。 expected_count は(isoformごとに振り分けられた)マップリード数、FPKMはリード数を全マップリード数と遺伝子長で正規化したスコアで、TPMはFPKMを更にFPKMの総和で正規化した値になります。

発現量の評価にはTPM、edgeRの入力にはcountを使おう

遺伝子発現量を直接見たい場合、例えば box plot や scatter plot を描く場合には TPMを使っておけば基本的に間違いありません。 一方、edgeRやDESeq2のような発現変動比較ツールにはTPMではなくcountデータを利用することが推奨されます（マニュアル参照）。何故なら、edgeRやDESeq2の中で用いている正規分布や負の二項分布はraw countに対してモデル化するものであり、遺伝子長などで正規化してしまうとこの分布を満たさなくなってしまうおそれがあるからです。全ての正規化はedgeRやDESeq2の中で行われますので、前もって正規化しておく必要はありません。更に言えば、サンプル間比較は遺伝子ごとに行われますので、基本的に遺伝子長の正規化は必要ないのです。

なお、DESeq2は整数値しか扱えませんので、DESeq2を利用する時にはcountデータの少数を切り捨てにしておきましょう。

発現量データをマージ (count)

それでは、４つのサンプルの発現量データをマージしてみましょう。RSEMの中にある rsem-generate-data-matrixを使います。 ここでは遺伝子ごとのファイルを用いています。

$ RSEM/rsem-generate-data-matrix   \
    star/Myers_H1-hESC_cell_2x75_200_1.genes.results  \
    star/Myers_H1-hESC_cell_2x75_200_2.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_1.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_2.genes.results  \
   > GeneExpressionMatrix.tsv

得られた結果を見てみましょう。

$ head GeneExpressionMatrix.tsv 
    "star/Myers_H1-hESC_cell_2x75_200_1.genes.results"   "star/Myers_H1-hESC_cell_2x75_200_2.genes.results"   "star/Myers_HUVEC_cell_2x75_200_1.genes.results" "star/Myers_HUVEC_cell_2x75_200_2.genes.results"
"ENSG00000000003"    1086.00 2327.00 1086.00 1891.00
"ENSG00000000005"    14.00   34.00   14.00   0.00
"ENSG00000000419"    355.00  755.00  355.00  875.00
"ENSG00000000457"    112.23  212.28  112.23  145.02
"ENSG00000000460"    501.77  878.72  501.77  201.98
"ENSG00000000938"    9.00    15.00   9.00    1.00
"ENSG00000000971"    0.00    0.00    0.00    7301.00
"ENSG00000001036"    1079.67 2036.77 1079.67 1889.43
"ENSG00000001084"    773.00  1263.00 773.00  608.00

行が遺伝子、列がサンプルの行列データになっていることがわかります。タブ区切りのテキストファイルなので、ExcelやR, Pythonなど好きなもので分析することが可能です。

発現量データをマージ (TPM)

上記で出力されているデータはカウントデータです。rsem-generate-data-matrixはカウントデータにしか対応していません（手抜き？）。

TPM や RPKM のマトリックスを出力したい場合は、僕が作成した rsem-generate-data-matrix-modified を使ってください。 rsem-generate-data-matrix-modifiedは、↓の記事で紹介しているscript_rnakato の中にあります。

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

ではやってみましょう。

$ script_rnakato/rsem-generate-data-matrix-modified \
    TPM \    # count または TPM または FPKM を指定
    star/Myers_H1-hESC_cell_2x75_200_1.genes.results  \
    star/Myers_H1-hESC_cell_2x75_200_2.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_1.genes.results  \
    star/Myers_HUVEC_cell_2x75_200_2.genes.results  \
    > GeneExpressionMatrix.TPM.tsv
# 結果を表示
$ head GeneExpressionMatrix.TPM.tsv 
    Myers_H1-hESC_cell_2x75_200_1.genes.results Myers_H1-hESC_cell_2x75_200_2.genes.results Myers_HUVEC_cell_2x75_200_1.genes.results   Myers_HUVEC_cell_2x75_200_2.genes.results
ENSG00000000003 40.20   41.67   40.20   39.30
ENSG00000000005 1.21    1.91    1.21    0.00
ENSG00000000419 32.34   34.32   32.34   39.18
ENSG00000000457 2.32    2.47    2.32    1.88
ENSG00000000460 18.38   15.50   18.38   4.85
ENSG00000000938 0.36    0.30    0.36    0.02
ENSG00000000971 0.00    0.00    0.00    93.63
ENSG00000001036 37.29   34.21   37.29   36.32
ENSG00000001084 26.16   25.42   26.16   11.66

無事TPMのマトリックスを得ることができました。ちなみにgene id の" も必要ないので除去しています。

おわりに

これで無事、fastqファイルから遺伝子発現量マトリックスを生成できるようになりました。 あとはこのtsvファイルをBioconductorなどで解析すればばっちりですね。 細かい点では、gene idを遺伝子シンボルに変更したいなどの要求があると思いますが、機会があれば追記します。

【Ubuntu 18.04】 Rのバージョンを3.5.2にアップグレード

その他 解析環境構築 Linux

最近知ったツールがRの3.5以上を要求しているのだが、手元のUbuntu 18.04はRが3.4.4 だったので、3.5にアップグレードしました。

以下備忘録。 やってることは、aptで参照されるレポジトリにR3.5を含むレポジトリを追加した上で改めてRをインストールするというものです。シンプル。

既存のRをアンインストール（必要ないかも？）

$ sudo apt remove --purge r-base r-base-core r-recommended
$ sudo apt autoremove

aptのレポジトリを追加

$ sudo add-apt-repository 'deb https://cloud.r-project.org/bin/linux/ubuntu bionic-cran35/'
$ sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 51716619E084DAB9
$ sudo apt update

Rをインストール

$ sudo apt install r-base r-base-core r-recommended r-base-dev

Rを起動して

update.packages()  # 既存パッケージのアップデート

余談

最初追加したレポジトリをbionic-cran35でなくxenial-cran35 にしていたせいで libpng12-0がないですとか色々依存関係のエラーになってしまった。Ubuntuのバージョンはちゃんと確認しよう。