R (tidyverse), inner_join

「計算用の数値を含むデータ」と、「各行を説明する名前 (symbol) や説明の文章などのデータ」という2つのテーブル同士を、それぞれに含まれる共通の列をキーに連結させたいことは、よくある処理だと思います。

このようなとき、R の tidyverse では、 inner_join() 関数を用いて、簡単に結合することができます。

例えば、下記のように2つのテーブルを読み込んだとします。

library(tidyverse)

input_data <- read_tsv("input_data.tsv")
meta_data  <- read_tsv("meta_data.tsv")

テーブルの中身は、下記のようなイメージです。2つのファイルも1列目に id の列を持ちます。

> # 計算用のデータ。
> input_data
# A tibble: 3 x 3
  Id    Sample1 Sample2
  <chr>   <dbl>   <dbl>
1 id1         1       4
2 id2         2       5
3 id3         3       6
>
>
>
> # 各行を説明するデータ。アノテーションの情報など。
> meta_data
# A tibble: 4 x 3
  Id    Symbol Description
  <chr> <chr>  <chr>
1 id1   AAA    lorem ipsum a
2 id2   BBB    lorem ipsum b
3 id3   CCC    lorem ipsum c
4 id4   DDD    lorem ipsum d

2つのテーブルの共通な列である Id をもとに、テーブルを結合するには、下記のように書けます。

result_data <- inner_join(input_data, meta_data)

実行結果は、下記のようになります。

> result_data <- inner_join(input_data, meta_data)
Joining, by = "Id"
> result_data
# A tibble: 3 x 5
  Id    Sample1 Sample2 Symbol Description
  <chr>   <dbl>   <dbl> <chr>  <chr>
1 id1         1       4 AAA    lorem ipsum a
2 id2         2       5 BBB    lorem ipsum b
3 id3         3       6 CCC    lorem ipsum c

Joining, by = “Id” となっていることで、 Id 列をもとに結合されたことが分かります。この場合は、2つのテーブルに同じ名前の列があるので、自動的に認識されています。明示的にキーの列を指定する場合は、inner_join(input_data, meta_data, by = "Id") と書きます。

また、結果のテーブル (result_data) では、元の input_data に含まれない id4 の行がないことが確認できます。

 

R (tidyverse), read_tsv

tidyverse には、タブ区切りテキスト (tsv) を読み込むための関数 read_tsv() が使えます。

library(tidyverse) # パッケージを読み込む。

input_data <- read_tsv("input_data.tsv")

上記の例では、input_data.tsv という名前のファイルを読み込んで、 input_data に格納しています。(csv ファイルであれば、read_csv() 関数で読み込みます。)

> input_data
# A tibble: 3 x 3
  id    sample1 sample2
  <chr>   <dbl>   <dbl>
1 AAA         1       4
2 BBB         2       5
3 CCC         3       6

読み込んだファイルは、従来のデータフレームではなく、 tibble という形式のオブジェクトになります。tibble は下記のような特徴があります。

  • 行名 (rownames) を持たない。(列名は持っている)
  • 列ごとに推論された型を持つ。
  • オブジェクトを表示させると、上から10行が表示される。

サイズの大きいデータの場合、データの先頭だけ表示されるのは便利です。読み込みも早いので、これだけでも tidyverse を使うメリットがあると思います。

よく使うオプションとしては、下記のものがあります。*その他は、help(read_tsv) で確認。

  • skip: ヘッダーが不要だったり、1行目がコメント行だったりするとき。
    • 例:skip = 1
  • col_names: 列名を指定したいとき。
    • 例:col_names = c("id", "sample1", "sample2")
  • col_types: 型推論がうまくいかなかったとき。
    • 例:col_types = cols(id = col_character())
input_data <- read_tsv("input_data.tsv", skip = 1,
                       col_names = c("id", "sample1", "sample2"),
                       col_types = cols(id = col_character(),
                                        sample1 = col_double()))

本家のチートシートが大変参考になります。

Data Import:: CHEAT SHEET

https://github.com/rstudio/cheatsheets/blob/master/data-import.pdf
 

R (tidyverse) の使い方

R のパッケージである tidyverse を使用すると、データの読み込み、変形、出力が簡単になります。特に、これから R を始められる方は、最初から tidyverse を使用した書き方で始める方が分かりやすいと思います。(細かいところは、従来のデータフレームの知識や書式が必要になるかもしれませんが。)

tidyverse を利用したコードのイメージ。

library(tidyverse)

input_data <- read_tsv("input_data.tsv")
meta_data  <- read_tsv("meta_data.tsv")

output_data <- input_data %>% inner_join(meta_data)
write_tsv(output_data, "output_data.tsv")

これだけのコードで、タブ区切りテキスト(TSV)形式のファイルを読み込んで、同じ列名のIDをキーに2つのテーブルを結合して、TSV形式で出力することができます。

短いコードで書くことができ、読みやすい表現が可能なことで、コードのメンテナンス性が良くなります。(多少、冗長な表現であったとしても、わかりやすい方が、後でコードを見返したときに困りません。)

single cell の解析に用いられているパッケージの Seurat など、最近のパッケージは、tidyverse をベースに使用しているため、R の従来の表現から、今後、tidyverse を使った新しい表現に切り替えていくことをお勧めします。

 

NGSデータ解析時のファイルのフォーマット (2)

FASTQ, BAM (SAM) に続いて、よく利用されるのが、BED (べっど)形式のファイルです。

BED

BED は、Browser Extensible Data の略です。Browser (ブラウザ)は、UCSC のゲノムブラウザのことです。染色体(ゲノム)上の位置情報に、何らかのスコアを加えて表現するためのファイルです。

これもテキスト形式のファイルで、1行ごとに、位置情報 (chrom, chromStart, chromEnd) とスコア (name, score)、その他の情報をオプションで記述します。

最低限のデータとしては、下記のようなものです。(各項目の区切りはタブ文字)

chr1 100 200 geneA 1000
chr2 1000 2000 geneB 10

chrom は染色体の番号、chromStart は、指定したい配列の開始位置、chromEnd は終了位置です。name は、ブラウザに表示するときのラベル名、score は、表示したい任意の数字 (0-1000) です。ここまでが、BEDファイルに必須の項目です。残りのオプション部分にいろいろな情報を指定できるので、同じ BED ファイルでも様々なバリエーションがあります。

https://genome.ucsc.edu/FAQ/FAQformat.html#format1

基本的には座標を指定するものなので、プロモーターやエンハンサーの位置を示したり、ChIP-seq の結果のピークの位置とスコアを表したり、様々なことに使用されています。

BED を元にさらに拡張されたフォーマットとして、bigBed, bedGraph があります。同様に、ゲノム上の位置(区間、領域)に情報を記述するものとして、WIG とそれを拡張、圧縮した bigWig などがあります。

bedtools

ゲノム上の位置情報どうしを演算するためのツールとして、 bedtools というソフトウェアがあります。bed で指定した配列の座標から、1000 base 上流の座標を求めたり、指定した座標の両端 500 から 1000 base の座標求めることができます。

さらにリファレンスゲノム(FASTAファイル)から、求めた座標に該当する領域の配列を抽出することも可能です。具体的な配列操作のイメージは、下記の解説が分かりやすいです。

getfasta

その他の操作もイメージつきで解説されています。(このようなソフトウェアの解説は、まず、作成者の準備したドキュメントに目を通すことをお勧めします。有用なことが多く書かれています。)

http://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html

Mac であれば、brew を使ってインストールできます。

brew install bedtools

* science タップを追加する必要があります。ない場合は上記のコマンドのログにその旨が表示されますので、それに従います。


 

NGSデータ解析時のファイルのフォーマット (1)

NGSのデータ解析時に用いられるファイルの形式は様々なフォーマットがあります。名前もよく似ているものが多く、初めは混乱するかもしれません。代表的なシーケンスのデータベースである UCSC に各種フォーマットの解説があります。

https://genome.ucsc.edu/FAQ/FAQformat.html

RNA-seq において、代表的なフォーマットは、FASTQ (ふぁすときゅー)とBAM(ばむ)です。

FASTQ

FASTQ は、シーケンサーで検出されたリードの配列を保存したものです。もとは、ATGCの文字を記述するだけの FASTA(ふぁすた)というフォーマットがあり、それに加え、シーケンサーで読み取った各塩基のクオリティーのスコアを格納したものです。リードを1本1本、記録しているため、数千万行から数億行の大きなサイズのファイルです。テキストファイルですので、テキストエディタで開くと、文字が読めます。(圧縮されている場合は、一度展開する必要があります。エディタによってはそのまま読めるものもあります。)

ファイルの拡張子は、.fastq, .fq などです。圧縮されていれば、.fastq.gz, .fq.gz のようになります。(FASTA の場合は、.fasta, .fa, fa.gz など。)また、_1.fq.gz, _2.fq.gz のように番号がついているものは、ペアエンドモードでシーケンスされた結果です。この場合、1サンプルあたり、2つのFASTQファイルがあります。

BAM (SAM)

シーケンサーから出力されたリード(FASTQ) を、リファレンスとなる配列にマッピングした結果は、SAM(さむ)形式のファイルになります。1本1本のリードが、ゲノム上のどこにマップされたかを示します。FASTQ の中身も含むので、ファイルのサイズは数十ギガバイト (GB) と大きくなります。ファイルの拡張子は、.sam です。SAMフォーマットの仕様は、GitHub で確認できます。

BAM は、上記の SAM を圧縮したものです。データの中身は、SAMと同一のものです。通常、マッピングした結果として提供されるのは、このBAMファイルです。1サンプルにつき1個の BAM ファイルになります(ペアエンドでもマッピング後のファイルは1個)。ファイルの拡張子は、.bam です。圧縮といっても、単純に gzip や、bzip2 で圧縮したわけではないため、専用のツール (samtools) を使って圧縮や展開を行います。(よって、sam.gz や sam.bz2 ではありません。念のため)

BAM は圧縮されているので、テキストエディタで開いても読めません。通常は、IGVなどのゲノムブラウザに読み込ませて、結果を確認することになります。(ゲノムブラウザでは、手動でSAMに変換することなく、BAMのまま読め場合がほとんどです。)

RNA-seq の場合、マッピング後の結果を比較することになりますが、単純に BAM ファイルどうしを比較するわけではありません。BAMファイルから、各遺伝子ごとに何本のリードがあるのかをカウントする作業が必要になります。