Gene annotation データを用意する（gtf形式）

RNA-seq解析の記事を書こうとしたらgtfファイルの部分が長くなり過ぎたので単独記事にしました。

Gene annotation

既知遺伝子の情報を記載するファイルにはいくつかの形式があり、gtf形式はそのひとつです。 よく似たgff形式というものもありますが、微妙に異なります。

gtf形式の詳細については以下の記事を参考にしてください（丸投げ）。

GTFファイル | 遺伝子アノテーションファイルの処理

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

遺伝子情報はダウンロード元のデータベースによってIDや情報の詳細度がかなり違います。 ここでは私がよく使うUCSC genome browserとEnsemblからダウンロードしてみたいと思います。

UCSC genome browser

Table browserで遺伝子情報を指定し、"get output" ボタンを押します。RefseqのHuman build hg38 であれば、以下のようになります。

Table Browser

Ensembl

Ensemblは各種データを↓のFTP siteからダウンロード可能です。

http://asia.ensembl.org/info/data/ftp/index.html

Human のGRCh38 (hg38に対応)データであれば、以下のコマンドでダウンロードできます。

$ wget ftp://ftp.ensembl.org/pub/release-94/gtf/homo_sapiens/Homo_sapiens.GRCh38.94.gtf.gz

scaffoldなどの遺伝子情報を含んだgtfファイルも別にありますが、通常は使いません。

ダウンロードした二つのファイルを解凍後、サイズを比較してみましょう。

$ gunzip *.gtf.gz
$ ls -lh *.gtf
-rwxrwxrwx 1 rnakato rnakato 1.1G Dec 13 11:17 Homo_sapiens.GRCh38.94.gtf
-rwxrwxrwx 1 rnakato rnakato  50M Dec 13 11:13 UCSC.hg38.gtf

Ensemblの方がファイルサイズが圧倒的に大きいですね。解凍すると1Gbを超えています。Ensemblの方がprocessed transcriptやpseudogeneなどの情報が多く、アノテーションも詳細に記述されていますので、遺伝子情報をもとに詳細な解析をしたい場合はEnsemblを使う方が良いと思います。

染色体名の違い

UCSC と Ensembl はgene IDが違いますが、染色体名も違います。UCSCの方は １番染色体は"chr1"となっていますが、Ensemblは単に"1"となっています。UCSCのgenomeファイルを用いて解析している場合、この違いを考慮しないツールであれば「染色体名が異なる」とエラーになり、これが地味に面倒です。
ゲノムの染色体名の方をEnsemblに揃えるやり方もありますが、ここではEnsemblの遺伝子名にchrを追加する方法を試してみましょう。

修正前

$ head Homo_sapiens.GRCh38.94.gtf
#!genome-build GRCh38.p12
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.27
#!genebuild-last-updated 2018-07
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 "processed_transcript"; 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 "processed_transcript"; 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 "processed_transcript"; 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 "processed_transcript"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";

冒頭が染色体名になっているので、行頭にchrを付加するだけでよいですね。これは以下のようなワンライナーで実行可能です。

# 冒頭のコメントアウト文を除外し、各行の行頭にchrを追加して出力
$ grep -v \# Homo_sapiens.GRCh38.94.gtf | awk '{print "chr" $0 }' > Homo_sapiens.GRCh38.94.addchr.gtf
# 内容確認
$ head Homo_sapiens.GRCh38.94.addchr.gtf
 head Homo_sapiens.GRCh38.94.addchr.gtf
chr1    havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
chr1    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 "processed_transcript"; tag "basic"; transcript_support_level "1";
chr1    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 "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    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 "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    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 "processed_transcript"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";
chr1    havana  transcript      12010   13670   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12010   12057   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12179   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12613   12697   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
chr1    havana  exon    12975   13052   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";

無事追加されました。 冒頭のコメント文を除外したくない場合は後から付け足すか、最初のコマンドでgrep部分を使わずに出力してください（後で手作業で#前のchrを取り除く）。

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

前回の記事の続きです。

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

発現量データ生成

前回の記事で発現量データを生成するのを忘れていたので、ここで生成します。stringtieを使います。

for prefix in ERR188044 ERR188104 ERR188234 ERR188245 ERR188257 ERR188273 ERR188337 ERR188383 ERR188401 ERR188428 ERR188454 ERR204916
do
    stringtie -e -B -p 16 -G stringtie_merged.gtf \
    -o ballgown/${prefix}/${prefix}_chrX.gtf \
    ${prefix}_chrX.bam
done

これでballgown/${prefix} というディレクトリが自動生成され、その中に発現量データが出力されます。
ディレクトリの中身を見てみましょう。

$ ls ballgown/ERR188044/
ERR188044_chrX.gtf  e2t.ctab  e_data.ctab  i2t.ctab  i_data.ctab  t_data.ctab

gtfファイルと、いくつかのテキストファイルができています。 eはexon, iはid, tがtranscriptの情報を表しているようですね。

比較するサンプルの確認

サンプルの情報はchrX_data/geuvadis_phenodata.csvに記述されています。 どんなサンプルなのか確認しましょう。

$ cat chrX_data/geuvadis_phenodata.csv
"ids","sex","population"
"ERR188044","male","YRI"
"ERR188104","male","YRI"
"ERR188234","female","YRI"
"ERR188245","female","GBR"
"ERR188257","male","GBR"
"ERR188273","female","YRI"
"ERR188337","female","GBR"
"ERR188383","male","GBR"
"ERR188401","male","GBR"
"ERR188428","female","GBR"
"ERR188454","male","YRI"
"ERR204916","female","YRI"

サンプルIDごとに、性別とpopulation情報が載っています。 解析内容は、population間のばらつきを除去しながら性別間の遺伝子発現変動を見る、というもののようです。

Ballgownで解析

BallgownはRで動きますので、同じディレクトリでRを起動しましょう。
以降のコマンドはR上で実行するものとします。Nature Protocolのコマンドから少し修正しています。

# ライブラリ読み込み
library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
# サンプル情報の読み込み(通常は自分で用意する)
pheno_data = read.csv("chrX_data/geuvadis_phenodata.csv")
# 先程生成した発現量データの読み込み
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)

# bg_chrXに格納されたものを確認
bg_chrX
ballgown instance with 3520 transcripts and 12 samples

bg_chrX という名前のballgownインスタンスが生成されました。

# 低発現遺伝子の除去
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)

# 発現変動遺伝子の推定（transcript単位）
results_transcripts = stattest(bg_chrX_filt,
                   feature="transcript",
                   covariate="sex",
                   adjustvars = c("population"),
                   getFC=TRUE,
                   meas="FPKM")

# 得られた発現変動遺伝子の確認
head(results_transcripts)
     feature id        fc      pval      qval
1 transcript  1 0.9347423 0.7043147 0.9373013
2 transcript  2 1.2006990 0.8716451 0.9723579
3 transcript  3 1.0173174 0.9896415 0.9990549
4 transcript  4 0.3814504 0.5200179 0.9192672
5 transcript  5 0.6014996 0.3216161 0.9192672
6 transcript  6 0.6540123 0.3264784 0.9192672

同様に、遺伝子単位の発現変動情報を抽出します。

results_genes = stattest(bg_chrX_filt,
             feature="gene",
             covariate="sex",
             adjustvars = c("population"),
             getFC=TRUE,
             meas="FPKM")

head(results_genes)
  feature         id        fc      pval      qval
1    gene    MSTRG.1 1.0824443 0.6847873 0.9196130
2    gene   MSTRG.10 0.9335310 0.8132130 0.9599962
3    gene  MSTRG.100 0.6126762 0.5687947 0.8823732
4    gene MSTRG.1000 1.2044223 0.4500794 0.8360675
5    gene MSTRG.1001 1.3799856 0.6140598 0.8936018
6    gene MSTRG.1003 0.9087684 0.5909100 0.8867180

以下ではresults_transcripts に遺伝子名の情報を付加します。

results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
                 geneIDs=ballgown::geneIDs(bg_chrX_filt),
                 results_transcripts)

head(results_transcripts)
  geneNames geneIDs geneNames.1 geneIDs.1    feature   id          fc
1         . MSTRG.2           . MSTRG.528 transcript 1668 0.030602170
2    PLCXD1 MSTRG.2        XIST MSTRG.528 transcript 1667 0.002990002
3         . MSTRG.2           . MSTRG.528 transcript 1666 0.015949119
4         . MSTRG.2           . MSTRG.528 transcript 1669 0.028359333
5         . MSTRG.3        TSIX MSTRG.527 transcript 1665 0.077438632
6    PLCXD1 MSTRG.2           . MSTRG.612 transcript 1862 7.335431029
          pval         qval
1 1.474474e-10 2.824911e-07
2 2.497711e-10 2.824911e-07
3 6.564345e-10 4.949516e-07
4 6.871521e-08 3.885845e-05
5 1.911435e-06 8.647331e-04
6 1.331176e-05 5.018534e-03

遺伝子名はid順にソートされていますが、発現変動遺伝子が見たいので、p値の小さい順にソートしなおします。

results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)

head(results_transcripts)
  geneNames   geneIDs    feature   id          fc         pval         qval
1         . MSTRG.528 transcript 1668 0.030602170 1.474474e-10 2.824911e-07
2      XIST MSTRG.528 transcript 1667 0.002990002 2.497711e-10 2.824911e-07
3         . MSTRG.528 transcript 1666 0.015949119 6.564345e-10 4.949516e-07
4         . MSTRG.528 transcript 1669 0.028359333 6.871521e-08 3.885845e-05
5      TSIX MSTRG.527 transcript 1665 0.077438632 1.911435e-06 8.647331e-04
6         . MSTRG.612 transcript 1862 7.335431029 1.331176e-05 5.018534e-03

head(results_genes)
  feature        id          fc         pval         qval
1    gene MSTRG.528 0.002409209 2.041411e-11 2.086323e-08
2    gene MSTRG.527 0.082689539 6.110906e-06 2.409310e-03
3    gene MSTRG.143 3.158100078 7.719104e-06 2.409310e-03
4    gene MSTRG.612 7.327271932 1.073484e-05 2.409310e-03
5    gene MSTRG.763 0.278924648 1.178723e-05 2.409310e-03
6    gene MSTRG.618 9.109365087 4.555636e-05 7.759766e-03

以下のコマンドを実行すると、得られたresults_transcripts, results_genesをテキストファイルに出力します。（ファイル名は任意）

write.csv(results_transcripts, "chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv", row.names=FALSE)

可視化

以下は可視化のコマンドです。得られた発現変動遺伝子についてさらに調査する時に用います。

# 色の設定
tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)

# 箱ひげ図
fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

ラベルが見切れていますが、サンプルごとの遺伝子発現量の分布が性別ごとに色分けされて箱ひげ図で表示されました。 縦軸のFPKMというのは遺伝子ごとにマップされたリード数を遺伝子長と全リード数で正規化したものです。

次に、特定の遺伝子について可視化してみましょう。 ここでは例として１２番めの遺伝子について可視化します。

ballgown::transcriptNames(bg_chrX)[12]
         12 "NR_027231" 
ballgown::geneNames(bg_chrX)[12]
         12 "LINC00685" 

おや、、、何故だかNature Protocolの例と表示される遺伝子名が違います。順番が変わったのでしょうか。
ともかく表示してみましょう。

plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
     main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
        ballgown::transcriptNames(bg_chrX)[12]),
     pch=19,
     xlab="Sex",
     ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
       col=as.numeric(pheno_data$sex))

サンプルごとの発現量のばらつきを性別で分けてプロットしているのですが、全てのサンプルで発現量が0のため、差が見えません。
この遺伝子ではない遺伝子が例になっていたはずなのですが。僕のスクリプトにどこか間違っているところを見つけたら教えてください。

plotTranscripts(ballgown::geneIDs(bg_chrX)[1729],
        bg_chrX,
        main=c('Gene XIST in sample ERR188234'),
        sample=c('ERR188234'))

遺伝子の簡単なisoform情報も表示することができます。

plotMeans('MSTRG.56', bg_chrX_filt, groupvar="sex", legend=FALSE)

ここではおそらく性別間で差がある遺伝子構造を可視化しているのですが、差がない遺伝子を指定していまっているのでよくわかりません。

おわりに

最後ミスったためしまりがないことになってしまいましたが、ともかく一通りの作業ができるようになりました。
新規転写物を探索したい時には有効なワークフローなので、機会があれば利用されると良いと思います。

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

せっかくNature Protocolの論文があるので試してみよう企画。
Nature Protocolはスクリプトがそのまま載っているので、追試に最適です。 一方、古いライブラリが指定されていると手元の環境で動かなかったり、著者のレベルによってへんてこなスクリプトになっていたりしますが…

今回使うのは以下の論文です。

www.nature.com

HISAT-StringTie-Ballgown については前回の記事を参照してください。

RNA-seqによる発現量解析 - Palmsonntagmorgen

hisat2, stringtieをインストール

hisat2, stringtie およびgffcompareも必要になりますので、コンパイル済ファイルをダウンロード、解凍します。

$ wget http://ccb.jhu.edu/software/hisat2/dl/hisat2-2.1.0-Linux_x86_64.zip
$ wget http://ccb.jhu.edu/software/stringtie/dl/stringtie-1.3.5.Linux_x86_64.tar.gz
$ wget http://ccb.jhu.edu/software/stringtie/dl/gffcompare-0.10.6.Linux_x86_64.tar.gz
$ unzip hisat2-2.1.0-Linux_x86_64.zip
$ tar zxvf stringtie-1.3.5.Linux_x86_64.tar.gz
$ tar zxvf gffcompare-0.10.6.Linux_x86_64.tar.gz

解凍したらパスを通します。 ~/.barhrcに以下のような感じで記述し（場所は自分のディレクトリに合わせる）、source ~/.barhrcを実行します。

export PATH=$PATH:$HOME/hisat2-2.1.0/:$HOME/stringtie-1.3.5.Linux_x86_64:$HOME/gffcompare-0.10.6.Linux_x86_64/

パスが通ったか確認しましょう。

$ hisat2 --version
/home/rnakato/git/binaries/hisat2-2.1.0/hisat2-align-s version 2.1.0
64-bit
Built on login-node03
Wed Jun  7 15:53:42 EDT 2017
Compiler: gcc version 4.8.2 (GCC) 
Options: -O3 -m64 -msse2 -funroll-loops -g3 -DPOPCNT_CAPABILITY
Sizeof {int, long, long long, void*, size_t, off_t}: {4, 8, 8, 8, 8, 8}

$ stringtie --version
1.3.5

$ gffcompare --version
gffcompare v0.10.6

バージョンが返ってきたら成功です。Not foundであればパスが通っていませんので.bashrcを見直しましょう。

Ballgownのインストール

BallgownはRで動きますので、R上で以下のコマンドを実行してインストールします。

> install.packages("devtools",repos="http://cran.us.r-project.org")
> source("http://www.bioconductor.org/biocLite.R")
> biocLite(c("alyssafrazee/RSkittleBrewer","ballgown","genefilter","dplyr","devtools"))

上記コマンドはNature protocolそのままですが、ballgownと共に必要になる関連ツール群もインストールされているようです。 最初にdevtoolsをインストールしていますが、これはgithub上のツールをインストールするためのツールです。 biocLiteはBioconductorからツールをインストールするコマンドで、sourceコマンドで参照元URLを指定します。

データダウンロード

次に使用するデータ群をダウンロードします。これは論文で提供されているものです。

$ wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz
$ wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/hg38_ucsc.annotated.gtf
$ tar zxvf chrX_data.tar.gz 

これはけっこう時間がかかります。
解凍されたディレクトリに何が入っているか確認しましょう。

$ ls chrX_data/*
chrX_data/geuvadis_phenodata.csv  chrX_data/mergelist.txt

chrX_data/genes:
chrX.gtf

chrX_data/genome:
chrX.fa

chrX_data/indexes:
chrX_tran.1.ht2  chrX_tran.3.ht2  chrX_tran.5.ht2  chrX_tran.7.ht2
chrX_tran.2.ht2  chrX_tran.4.ht2  chrX_tran.6.ht2  chrX_tran.8.ht2

chrX_data/samples:
ERR188044_chrX_1.fastq.gz  ERR188257_chrX_1.fastq.gz  ERR188401_chrX_1.fastq.gz
ERR188044_chrX_2.fastq.gz  ERR188257_chrX_2.fastq.gz  ERR188401_chrX_2.fastq.gz
ERR188104_chrX_1.fastq.gz  ERR188273_chrX_1.fastq.gz  ERR188428_chrX_1.fastq.gz
ERR188104_chrX_2.fastq.gz  ERR188273_chrX_2.fastq.gz  ERR188428_chrX_2.fastq.gz
ERR188234_chrX_1.fastq.gz  ERR188337_chrX_1.fastq.gz  ERR188454_chrX_1.fastq.gz
ERR188234_chrX_2.fastq.gz  ERR188337_chrX_2.fastq.gz  ERR188454_chrX_2.fastq.gz
ERR188245_chrX_1.fastq.gz  ERR188383_chrX_1.fastq.gz  ERR204916_chrX_1.fastq.gz
ERR188245_chrX_2.fastq.gz  ERR188383_chrX_2.fastq.gz  ERR204916_chrX_2.fastq.gz

ゲノム配列や遺伝子ファイル、hisat2用のインデックス、fastqデータなどが入っています。 ファイルサイズと計算量を抑えるために、chrXだけを抽出したもののようですね。

マッピング

ツールとデータが揃いましたので、早速実験しましょう。まずはhisat2によるマッピングです。

for prefix in ERR188044 ERR188104 ERR188234 ERR188245 ERR188257 ERR188273 ERR188337 ERR188383 ERR188401 ERR188428 ERR188454 ERR204916
do
    hisat2 -p 16 --dta -x chrX_data/indexes/chrX_tran \
       -1 chrX_data/samples/${prefix}_chrX_1.fastq.gz \
       -2 chrX_data/samples/${prefix}_chrX_2.fastq.gz \
       -S ${prefix}_chrX.sam
    samtools sort -@ 16 -o ${prefix}_chrX.bam ${prefix}_chrX.sam
done

論文ではわかりやすさのため１サンプルずつ順番に書き下していますが、面倒なので、ここではfor文を使って全サンプルを順にマッピングしています。
${prefix} のところにERR188044, ERR188104, ... と順に入っていく感じです。
論文ではCPU数を8に指定していましたが、早く計算したいのでここでは16にしています。
ついでにsamtoolsによるsort.bamへの変換も同時に行っています。
手元の環境では１サンプルにつき数十秒で計算終了しました。

遺伝子ファイル作成、発現量計測

次にStringtieを使って発現量を計測します。こちらもfor文で回します。

for prefix in ERR188044 ERR188104 ERR188234 ERR188245 ERR188257 ERR188273 ERR188337 ERR188383 ERR188401 ERR188428 ERR188454 ERR204916
do
    stringtie -G chrX_data/genes/chrX.gtf -o ${prefix}_chrX.gtf \
              -p 16 -l ${prefix} ${prefix}_chrX.bam
done

入力がbamファイル、出力がgtfファイルです。 完了したら、各サンプルのgtfファイルをマージしましょう。

$ stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf \
 -o stringtie_merged.gtf chrX_data/mergelist.txt

入力となる各サンプルのgtfファイルは、chrX_data/mergelist.txt内に記載されています。
ここでは著者が親切にリストを既に作ってくれていますが、通常は自分で新規作成するものです。

遺伝子ファイルの比較

論文ではoptionalとなっていますが、gffcompareを使って、生成された遺伝子リストとreferenceの遺伝子リストを比較することができます。

$ gffcompare -r chrX_data/genes/chrX.gtf -G -o merged stringtie_merged.gtf

この作業により、referenceに含まれない新規転写物を探索することができます。

次回に続きます。

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

RNA-seqによる発現量解析

本業の方で色々忙しくなっておりまして、更新の間が開いてしまいました。

今回はRNA-seqについて語りたいと思います。 RNA-seqはChIP-seqよりもメジャーなので、日本語での解説ブログも充実していますが、情報が古いものだと今だにtophat-cufflinksを使っていたりします。それじゃだめだよ！というお話。

ここではツールのレビューがメインです。RNA-seq解析のHowtoそのものについては下記サイトが大変わかりやすいですので、参照されてください。RベースでのRNA-seq解析法です。

RNA-Seq | 遺伝子発現量解析

RNA-seq

RNA-seqの目的は大きく分けて、

• 既知遺伝子の発現量を網羅的に計測し、複数サンプルで発現変動比較する
• 未知の転写物（isoform含む）を同定する

の２つがあり、前者のケースが大半かと思います。
前者の場合、既知RNAのデータベースに対してリードをマッピングします。後者の場合はゲノムに直接マッピングした後アセンブリングするか、まずリードをアセンブリングしてからゲノムにマッピングすることになります。

RNA-seqでの発現量解析

RNA-seqの一般的なプロトコルは、

• リードを遺伝子・ゲノムにマッピング（マッピング）
• 遺伝子ごとにマップリード数をカウント、正規化（発現量計測）
• サンプルごとに発現量を比較（発現変動解析）

の３ステップに分けられ、それぞれ異なるプログラムが必要です。 ステップの詳細が知りたい方はConesaらの論文*1を参照されると良いでしょう。

Tophat-cufflinks

Tophatはリードをマッピングするツール、cufflinksは発現量計測、発現変動解析、新規転写物の同定などのコマンドがセットになったプログラム集です。
tophat-cufflinksは様々な目的に利用可能であること、解説ページが豊富に公開されていたこと、1 replicateでも発現変動比較ができるといったメリットがあり、一時期はかなり普及していました。

その後tophat-cufflinksは開発・サポートが終了し、後継のhisat2またはkallistoに切り替えることが開発者より推奨されています。精度面の問題も明らかになってきており、今ではほとんど使われることもなくなったのですが、論文の査読をしているとtophat-cufflinksを使っているものに未だに出くわします。初心者ならともかく、ある程度サーベイして知識を持った人が未だに使っているのは正直勉強不足だろうと思いますが。。。

最近のツール事情

RNA-seqのツール開発は今でも群雄割拠の様相を呈しており、ツールごとの精度比較を行ったサーベイ論文がいくつか発表されています*2*3。この論文によれば、cufflinksはsingle exonの遺伝子に対して非常に精度が低くなることが報告されています。また、single-end RNA-seqでisoformレベルの発現量を推定した場合にも精度が悪いようです。また、１サンプルからの発現変動比較が可能ですが、１サンプルだと当然精度が悪いです。 従って、paired-endを使って遺伝子レベルの発現量比較をしている場合のみ、tophat-cufflinksはいちおう今でも利用可能ですが、特に使い続けるメリットはありません…査読者にも必ず突っ込まれるでしょう。

ツールの推奨組み合わせ

基本的にはマッピングー発現量計測ー発現変動解析の流れで用いるツール群はセットになっています。 検証されていない組み合わせを用いると予期せぬエラーが起きる可能性がありますので、素直に推奨組み合わせを使いましょう。 具体的には以下のような組み合わせがあります（/ で挟んでいるものはどちらも利用可という意味です）。

• (bowtie2/STAR)ーRSEMー(edgeR/DESeq2)
• Hisat2ーStringTieーBallgown
• kallistoーsleuth
• Salmonー(edgeR/DESeq2)

いわゆる正統派の「マッピングー発現量計測ー発現変動解析」という意味では、 (bowtie2/STAR)-RSEM-(edgeR/DESeq2)を使うのが良いと個人的には思います。精度は最高レベルで、single-endでも利用可、遺伝子単位にもtranscript単位にも対応しており、利用実績も多いです。bowtie2とSTARは精度はそれほど変わりませんが、STARの方が高速です（そのかわりメモリ消費がとても大きい）。edgeRとDESeq2はどちらも優れたツールなので、お好みでよいと思います。他にもいくつかツールはありますが、精度面でこれらを大きく超えることはありません。TIGAR2という国産の発現量計測ツールもありますが、ものすごく遅い割に精度はそんなに変わりません。。

Hisat2-StringTie-Ballgownは新規転写物の同定が可能で、プロトロルがNature protocolsに公開されています*4。一方、既知の遺伝子発現計測では精度はそれほどでもないです。 グラフゲノムアラインメントというコンセプトが取られており、最近はメタゲノム（細菌叢）の解析への応用が検討されています。 なお、eXpressという類似のツールがありますが、こちらはサポート終了していますので使わないようにしましょう。

kallistoとSalmonは"pseudo-alignment (偽アラインメント)"と言って、リードを遺伝子にひとつひとつマップすることなく遺伝子発現を計測可能となったツールです。そのため大変高速です（試してみるとびっくりします）。これはひとつのエポックメイキングであり、精度面でも最高レベルに近いため、「これからはpseudo-alignment だろう！」という機運はあるのですが、デフォルトではBAMファイルを生成してくれない（マッピングしないので）ことや、利用実績が少なくちょっと怖いこともあり、今のところ様子見している人が多いように感じます。なおSailfishというツールもありますが、これはSalmonの先行ツールです。
kallistoは少し癖があります。基本的にはpaired-endのみに対応で、sleuthという独自の発現変動解析ツールがペアになっており、これ以外は使わないことが推奨されています（組み合わせた場合にどうなるのかは未確認です）。そのためedgeRなどを使いたい人は使えないということになってしまいますね。
また、kallisto・Salmonともtranscript単位でしか発現量を出力してくれません。遺伝子単位での発現量マトリクスを生成したい場合はtximport*5というツールを併せて使いましょう。

まとめ

ここではマッピングー発現量解析ー発現変動解析のためのツールを簡単に紹介しました。 他にも、正規化手法やread assembling, alternative splicingなどさまざまな目的の手法が本当にたくさんあり、全てサーベイするのは難しいほどです。

ちなみに私はSTAR-RSEM-edgeRを使っています。pseudo-alignmentの高速性は大変魅力的なのですが、万一後から大問題が発覚すると怖いので、今のところ保守的に行っとこうという感じですね。。。お試しレベルではkallistoもSalmonも使ってはいます。おそらく二年後くらいにはまた様子がすっかり変わっているのでしょうね。

*1:A. Conesa et al., A survey of best practices for RNA-seq data analysis, Genome biology 2016, DOI 10.1186/s13059-016-0881-8

*2:A. Kanitz et al., Comparative assessment of methods for the computational inference of transcript isoform abundance from RNA-seq data, Genome Biology, 2015, DOI 10.1186/s13059-015-0702-5

*3:M. Teng et al., A benchmark for RNA-seq quantification pipelines, Genome Biology, 2016, DOI 10.1186/s13059-016-0940-1

*4:M. Pertea et al., Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, Nature Protocols, 2016, doi:10.1038/nprot.2016.095

*5:http://bioconductor.org/packages/release/bioc/html/tximport.html

2サンプル間ピーク比較

２つのサンプルから得られたピークセットがどのくらい重なるのか調べたい！という時の方法です。

今回はBEDtoolsを使うやり方と、拙作のcompare_bsを使うやり方を紹介します。

ピークデータのダウンロード

ピークはBED形式であれば何でもよいのですが、ここでは例としてUCSCからH1細胞のCjunとMaxのピークファイルをダウンロードします。
(.narrowPeak という拡張子のファイルですが、これらはMACS2で得られたピークファイルで、BED形式として扱えます。)

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak.gz

解凍後、ピーク数を確認してみましょう。

 $  gunzip *.gz
 $  wc -l wgEncode*
   2148 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
  11129 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
  13277 total

Cjunのピークはやや少ないようです。

bedtools

bedtoolsのintersectを使うと、共通するピーク領域を出力してくれます。その後数をカウントすればピーク数が判明します。

（Bedtoolsの詳細はこちら: BEDtoolsワンライナー覚書 - Palmsonntagmorgen）

 $ bedtools intersect \ 
    -a wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \
    -b wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
     > and.bed
 $  wc -l and.bed
377 and.bed

Cjunの2148ピークのうち、Maxと重なるものは377ピークになりました。

bedtools intersectは２つのピークが1塩基でも重なっている場合にoverlapと判定します。 また、サンプルAの１つのピークに大してサンプルBの３つのピークが重なる場合、出力されるピーク数は３つになります。コマンドに -wa オプションをつけると、1つのピークとして出力されます。
intersectで出力される領域については、以下に図解があります。
intersect — bedtools 2.29.2 documentation

compare_bs

compare_bsは以下の記事で説明したChIPseqTools の中にあります。

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

compare_bsを使うと以下のようになります。

 $ compare_bs -and \
    -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \
    -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
    > and.bed
 $ head and.bed
#file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
#file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
#num1: 2148     num2: 11129     num1_overlap: 373 (17.4%)       num1_notoverlap: 1775 (82.6%)   num2_overlap: 378 (3.4%)num2_notoverlap: 10751 (96.6%)
#peakwidth total1: 523816 bp    peakwidth total2: 3534874 bp    overlappeaks total: 70413 bp (13.44% / 1.99%)
chromosome      start   end     summit  length  enrich  intensity
chr10   80917404        80917527        80917465        123     0.00    149.81
chr22   39570815        39571059        39570937        244     0.00    142.53
chr18   9643685 9643754 9643719 69      0.00    133.02
chr6    20811494        20811738        20811616        244     0.00    132.72
chr1    205263766       205264010       205263888       244     0.00    116.91

compare_bsの出力では、冒頭にピーク数のスタッツが表示されます。 この場合、Cjunのピーク数2148、Maxのピーク数11129で、CjunのうちMaxと重なるものが373ピーク、MaxのうちCjunと重なるものが378ピークとわかります。
(CjunとMaxでそれぞれ重なったピーク数がずれているのは、ピークの重なりが1対1対応でないことに起因します。つまり１つのピークに対して複数の相手方ピークが重なる可能性があるということです。)

その下の行の #peakwidth total は、ピークの数でなく延べピーク幅で重なりをカウントしたものです。

compare_bs はデフォルトでは「サンプル１のピークのうち、サンプル２と重ならないもの」を出力し、-and オプションをつけると重なるものを出力します。 ですので、出力されるピークは bedtools intersect の -waオプションと同じく、サンプル1のピーク幅になります。

重なりの定義もbedtools と同じく、1塩基でも重なればoverlapと判定します。 -l 10000 とオプションをつけると、10kbpまで離れたピークをoverlapとして判定します。

 $  ~/git/ChIPseqTools/bin/compare_bs -and -l 10000  \
     -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak  \
     -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
     -nobs
#file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
#file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
#num1: 2148     num2: 11129     num1_overlap: 689 (32.1%)       num1_notoverlap: 1459 (67.9%)   num2_overlap: 787 (7.1%)num2_notoverlap: 10342 (92.9%)
#peakwidth total1: 523816 bp    peakwidth total2: 3534874 bp    overlappeaks total: 6727767 bp (1284.38% / 190.33%)

ここでは-nobs オプションをつけて、ピークリストを出力せずに冒頭のスタッツのみを表示しています。 10kbpまで重なりの閾値をゆるめたことで、少し重なりが増えました。
（注：overlapped peakwidthは10kbp extensionも考慮してoverlapを出すので、パーセンテージは少しおかしくなってしまいます）

まとめ

bedtools, compare_bsで得られる結果を図にまとめると以下のようになります。

得られるピークサイトはほぼ同じですが、compare_bsの優位点としては、statsを出してくれることとサンプル１と２のoverlap peak数をそれぞれ出してくれること、peak widthでもstatsを出してくれることがあります。
ピーク比較の時にはピークの幅は一定であると暗黙に仮定してしまいがちですが、幅の分布はかなり広がっているのが普通で、 １つの大きなピークが細切れになってカウントされていたりすることもあるので、ピーク数をカウントしてものを言いたい時には実際に得られたピークを確認する作業を忘れないようにしましょう。

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

NGS解析のための新規ツールは日々論文で発表されており、それらのほとんどは世界中の人が無償で利用可能なライセンス形態になっています。
今日はその中でも多くの人に利用されている「GitHub」に公開されたツールのインストール方法を紹介します。

オープンソースについて

今のBioinformatics界は極めてオープンソース志向が強く、「良いものや知見は共有されることでより良いものができる」という性善説的な発想に基づいています。bioRxivのようなopen preprintにも積極的です。
ので、新規のプログラムで論文を書く場合はそのプログラムがオープンソースで公開されていることが必須条件となっている論文誌がほとんどです。オープンにする際も、自分のラボページにzipファイルをアップするだけ、のようなあやしげな状態ではなく、GitHubやBioconductorのような広く使われバージョン管理がしっかりしているコミュニティスペースに公開されることが推奨されます。Bioinformatics誌のレビュアーの一人は、「GitHubに公開されていないプログラムはその時点でrejectにする」とツイートしていました。

一方、たとえばNatureやCellのような大きな論文内で使わている個々の解析技術はin-houseで作られたスクリプトであることが多く、そういうのはあまり公開はされません。Wetのラボはソースの公開に消極的なところが多い印象です。

GitHub

このブログでは僕のツールを使うことが多いので、僕のアカウント(rnakato)を例に以下説明します。 僕のGitHubサイトはこちらです。DROMPAやSSPなど、今まで開発したツールやスクリプトを置いています。

github.com

GitHubについて

Git・GitHubとは非常に平たく言うと、作成したソースコードをネット上に公開して誰でもダウンロードできるようにするためのしくみです。ファイルのバージョン管理機能に優れており、複数の人が同時にソースを修正するような開発に特に向いています。
（どこかの新聞記事では「設計図共有サイト」とか言っててコーヒー吹いた）
より詳しく知りたい場合は「GitHub 入門」でググってください。

GitHub上のツールをダウンロードするだけならばアカウントは作成しなくとも大丈夫です。

GitとGitHubのちがい

Gitはファイルのバージョン管理のしくみそのもので、GitHubはGitをWeb上で利用するためのサービス（Gitホスティングサービス）です。 GitHub以外にもBitbucketやGitLabなどいくつかのGitホスティングサービスがあります。

Gitのインストール

gitとコマンドを打ってみて、「そんなツールないよ」とエラーが出たら、以下のコマンドでインストールしましょう。

 $ sudo apt install git

GitHubからインストール（スクリプトの場合）

スクリプトとは、そのファイルをダウンロードすればコンパイル無しでそのまま使えるソースコードです。 私が作っているスクリプトを以下で公開しています。（自分用なので、大半は使い物にはならないと思いますが、一部は有用と思います。）

GitHub - rnakato/script_rnakato

上のリンクを開くと↓のような画面になります。画面右にある"Clone or download"をクリックするとクローン用URLが表示されますので（図の赤枠部分）、コピーしてください（URL右のアイコンをクリックすると自動でコピーされます）。

またはURL下の"Download zip"をクリックすると、圧縮ファイルとしてダウンロードすることもできます。

git clone <コピーしたURL> とコマンド実行すると、内容がクローン（ダウンロード）されます。

 $ git clone https://github.com/rnakato/script_rnakato.git
Cloning into 'script_rnakato'...
remote: Counting objects: 1674, done.
remote: Compressing objects: 100% (25/25), done.
remote: Total 1674 (delta 15), reused 17 (delta 7), pack-reused 1642
Receiving objects: 100% (1674/1674), 6.72 MiB | 3.53 MiB/s, done.
Resolving deltas: 100% (1002/1002), done.
Checking connectivity... done.
 $ ls  # ダウンロードされたフォルダを確認
script_rnakato

これで、script_rnakatoフォルダの中にあるスクリプトが全てダウンロードされました。

内容の更新

git cloneでクローンしたフォルダの内部でgit pull origin masterとコマンドを打つと、クローン後に加えられたGitHub上の更新内容を取り込むことができます。

 $ cd script_rnakato
 $ git pull origin master
From https://github.com/rnakato/script_rnakato
 * branch            master     -> FETCH_HEAD
Already up-to-date.

今はクローンしたばかりですので、当然修正はありません。

GitHubからインストール（プログラムの場合）

DROMPA3をダウンロード・インストールしてみます。 スクリプトと同じ要領ですが、C言語ソースのため、ダウンロード後にコンパイルが必要です。

GitHub - rnakato/DROMPA3: DROMPA version 3 (Peak-calling, Visualization, Normalization and QC for ChIP-seq analysis)

 $ git clone https://github.com/rnakato/DROMPA3
 $ cd DROMPA3
 $ make

GitHubからインストール（submoduleを含むプログラムの場合）

以下のChIPseqToolsはC++プログラムですが、レポジトリ内部にSSPのレポジトリを含んでいます。 これはsubmoduleという機能で入れ子構造になったGitHubレポジトリです。

GitHub - rnakato/ChIPseqTools

submoduleを含むレポジトリをクローンする場合は、--recursiveオプションをつけます。

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

また、内容を更新したい場合は以下のようにします。

 $ git pull origin master    # ChIPseqToolsレポジトリの更新
 $ git submodule foreach git pull origin master # ChIPseqTools内のSSPレポジトリの更新

まとめ

３種類のパターンに分けてGitHubからのツールのインストール方法を紹介しました。インストール方法は通常はページ下部のREADMEに書いてあるので、それを見れば大丈夫です。 NGS解析のツールの大半はGitHub上で公開されているので、GitHubを使いこなせるようになると自由度が各段にアップします。

またGitHub上に自分のアカウントを作成すれば、自分のスクリプトをアップすることもできるようになります。 GitHubは常に最新バージョンが保てる点がメリットなので、例えば職場と自宅でスクリプトを共有したいというような目的にも使えますし、 外部の人にスクリプトを送るのにも便利です。
(Dropboxでも共有はできますが、ファイル喪失の危険はありますし)

なおGitHubは無料プランだとファイル公開しか選べません。もしWeb上にアップしたいけれども公開にはしたくないという場合にはBitbucketかGitLabを利用しましょう。 （2019/1/16追記　GitHubは無料アカウントでもプライベートリポジトリを作成できるようになりました。）

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

公開されているゲノム配列は現在も更新中であるため、いくつかのバージョン (build) があります。 humanだとhg18, hg19, hg38などがあり、hg38が現時点で最新です。
NGS解析をするうえでは全ての解析データのbuildを統一する必要がありますが、「既存論文のデータの一部を用いたい」というような場合、buildが合わずに困ることがあります。 「自分はhg38で統一しているのに、既存論文のBEDファイルはhg19しか提供されてない」というような場合ですね。

そういった場合には LiftOver というツールを使うと、バージョンに合わせてBEDファイルを再生成してくれます。 詳細は以下のリンクを参照してください。

LiftOver - Genome Analysis Wiki

（話がそれますが、hg38がリリースされて既に久しいのですが、2018年の論文でも普通にhg19が使われています。 おそらくhg19で生成されたデータベース、アノテーションの類があまりに多すぎるため、 hg38にアップデートできない状態なのではないかと思います。ROADMAPもJuicerToolsもhg19ですしね。 既存データベースのデータを多く使う実験を組む際にはbuildをhg19にするのも手かもしれません。）

LiftOver (web)

hg19のBEDファイル(sample.bed)をhg38に変換する例です。
最も簡単なのは、UCSC web browser上で変換する方法です。

Lift Genome Annotations

デフォルトではhg19 -> hg38への変換になっていますので、そのまま使いましょう。 sample.bedをアップロードするか、bedの内容を直接テキストボックスに貼りつけてsubmitします。

submit すると、下にResultsができ、生成されたBEDファイルをダウンロードすることができます。

LiftOver (コマンドライン)

LiftOverはコマンドラインツールでも提供されています。 変換したいファイルが大量にあるような場合はこちらのほうが便利でしょう。

以下のようにしてバイナリファイルをダウンロードしましょう。

wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
chmod +x liftOver  # 実行権限の付与

変換にはLiftOverツールそのものの他に、変換用のchainfileが必要になります。

以下のディレクトリから、hg19 -> hg38 用のchainfile hg19ToHg38.over.chain.gz をダウンロードしましょう。

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/

ダウンロード後、解凍しておきます。

gunzip hg19ToHg38.over.chain.gz 

liftOver コマンドは、 liftOver <入力ファイル> <chainfile> <出力ファイル> <変換失敗したサイトを格納するファイル>というかたちで実行します。 今回の場合は以下のコマンドで変換完了です。

liftOver sample.bed hg19ToHg38.over.chain output.bed unlifted.bed

unlifted.bedには、何らかの理由で変換できなかったbedサイトが格納されます。変換できない理由は以下を参照してください。

https://genome.sph.umich.edu/wiki/LiftOver#Various_reasons_that_lift_over_could_fail

注意点

LiftOverを使っての変換は同一生物種の異なるbuildが想定されており、異なる生物種への変換は可能ではあるものの推奨はされていません。

また、変換失敗するサイトが発生することからもわかるように、変換は厳密なものではありません。 例えば、hg19でピークコールしたものをhg38に変換したピークリストと、hg38でピークコールしたピークリストはけっこう違います。遺伝子のアノテーションなどもかなり変わっています。
ので、大規模に統合したい場合には面倒でもfastqファイルから（条件を統一して）解析しなおす方がおすすめです。open chromatin regionのリストが欲しい、くらいであれば問題ないとは思います。