Aurora blog

バイオインフォ・バイオテクノロジーにまつわる情報

R:オミクスデータのためのクラス (SummarizedExperiment, MultiAssayExperiment)

はじめに

最近一つの実験で複数種類の高次元データを取得する機会や、ATAC-seqやChIP-seqなどのパラメータがゲノム領域に対応するデータ*1を扱う機会が増えてきた。これまで単一のデータをRで解析する際にはdata.frameクラスやmatrixクラスでデータを扱ってきたが、不便さを感じる場面が出てきたので、少しずつ別のクラスへと移行していこうと考えている。R/BioconductorではSummarizedExperimentMultiAssayExperimentなどの便利なクラスが提供されている。ここではその機能についてまとめようと思う。

SummarizedExperiment

SummarizedExperimentクラスは下図のような構造を有する。遺伝子発現などの行列データに加えて、行 (遺伝子名やピーク) と列 (サンプル) の情報、実験背景等のメタデータを一つのオブジェクトでまとめて扱うことができる。

Biobaseライブラリが提供するExpressionSetクラス*2も類似した機能を有し、広く使われているが、パラメータの情報をGRangesクラスで格納できるためATAC-seq等のデータへ親和性が高いことや、異なる方法で前処理した行列データ *3 を一つのオブジェクトへまとめて格納できるなどの優位性がある。個人的にはSummarizedExperimentを使おうと思う。

f:id:auroratummy:20210222073633p:plain

引用元:http://bioconductor.org/packages/release/bioc/vignettes/SummarizedExperiment/inst/doc/SummarizedExperiment.html#anatomy-of-a-summarizedexperiment

SummarizedExperimentオブジェクトの例

library(SummarizedExperiment)
data(airway, package="airway") # test data in airway package
airway
## class: RangedSummarizedExperiment 
## dim: 64102 8 
## metadata(1): ''
## assays(1): counts
## rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
## rowData names(0):
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(9): SampleName cell ... Sample BioSample

数値データを取得

assay(airway, "counts") %>% head()
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003        679        448        873        408       1138       1047        770        572
## ENSG00000000005          0          0          0          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587        799        417        508
## ENSG00000000457        260        211        263        164        245        331        233        229
## ENSG00000000460         60         55         40         35         78         63         76         60
## ENSG00000000938          0          0          2          0          1          0          0          0

行(遺伝子)の情報を取得

rowRanges(airway)
## GRangesList object of length 6:
## $ENSG00000000003
## GRanges object with 17 ranges and 2 metadata columns:
##        seqnames            ranges strand |   exon_id       exon_name
##           <Rle>         <IRanges>  <Rle> | <integer>     <character>
##    [1]        X 99883667-99884983      - |    667145 ENSE00001459322
##    [2]        X 99885756-99885863      - |    667146 ENSE00000868868
##    [3]        X 99887482-99887565      - |    667147 ENSE00000401072
##    [4]        X 99887538-99887565      - |    667148 ENSE00001849132
##    [5]        X 99888402-99888536      - |    667149 ENSE00003554016

列(サンプル)の情報を取得

colData(airway)
## DataFrame with 6 rows and 9 columns
##            SampleName     cell      dex    albut        Run avgLength Experiment    Sample    BioSample
##              <factor> <factor> <factor> <factor>   <factor> <integer>   <factor>  <factor>     <factor>
## SRR1039508 GSM1275862   N61311    untrt    untrt SRR1039508       126  SRX384345 SRS508568 SAMN02422669
## SRR1039509 GSM1275863   N61311      trt    untrt SRR1039509       126  SRX384346 SRS508567 SAMN02422675
## SRR1039512 GSM1275866  N052611    untrt    untrt SRR1039512       126  SRX384349 SRS508571 SAMN02422678
## SRR1039513 GSM1275867  N052611      trt    untrt SRR1039513        87  SRX384350 SRS508572 SAMN02422670
## SRR1039516 GSM1275870  N080611    untrt    untrt SRR1039516       120  SRX384353 SRS508575 SAMN02422682
## SRR1039517 GSM1275871  N080611      trt    untrt SRR1039517       126  SRX384354 SRS508576 SAMN02422673

特定のサンプル・遺伝子を抽出

airway[1:3, 1:3]
## class: RangedSummarizedExperiment 
## dim: 3 3 
## metadata(1): ''
## assays(1): counts
## rownames(3): ENSG00000000003 ENSG00000000005 ENSG00000000419
## rowData names(0):
## colnames(3): SRR1039508 SRR1039509 SRR1039512
## colData names(9): SampleName cell ... Sample BioSample

オブジェクトの作成

# 遺伝子発現データ (pasillaパッケージのデータを使用)
count_file = system.file("extdata", "pasilla_gene_counts.tsv", package="pasilla", mustWork=TRUE)
count = read.table(count_file, sep="\t", header=T, row.names=1)
cpm = sweep(count, 2, colSums(count), "/") * 1000000
# サンプルデータ
annot_file = system.file("extdata", "pasilla_sample_annotation.csv", package="pasilla", mustWork=TRUE)
annot = read.csv(annot_file, header=T, row.names=1)
# SummarizedExperiment オブジェクトを作成
data = SummarizedExperiment(assays = list(count=count, cpm=cpm), colData = annot)
data

bioconductor.org

MultiAssayExperiment

SummarizedExperimentクラスは複数の行列データを格納することが可能だが、各データの行数・列数が等しくなければならないという制約がある。当然だがオミクスデータは種類によって次元数は異なり、パブリックのマルチオミクスデータは一部のデータが一部のサンプルで未計測であることも多く、多種類のデータを扱う場合はSummarizedExperimentクラスでは対応できないことが多い。MultiAssayExperimentは多種類データに向けてデザインされたデータ型である。

f:id:auroratummy:20210222092644p:plain

引用元:https://www.bioconductor.org/packages/release/bioc/vignettes/MultiAssayExperiment/inst/doc/MultiAssayExperiment.html

www.bioconductor.org

TCGAやcBioPortalのデータをMultiAssayExperimentクラスで利用可能な形でダウンロードできるRパッケージcuratedTCGAData cBioPortalDataも公開されている。

bioconductor.org

www.bioconductor.org

MultiAssayExperimentオブジェクトの例

library(MultiAssayExperiment)
library(curatedTCGAData)
brca = curatedTCGAData("BRCA", c("RNASeqGene", "miRNASeqGene", "Mutation"), FALSE)
brca

マルチオミクスデータへアクセス

各オミクスデータはSummarizedExperimentクラスやRaggedExperimentクラスで格納される

experiments(brca)
# experiments(brca)[[1]] # 1番目のデータへアクセス
# brca[ , , 1] # 1番目のデータへアクセス
## ExperimentList class object of length 3:
##  [1] BRCA_miRNASeqGene-20160128: SummarizedExperiment with 1046 rows and 849 columns
##  [2] BRCA_Mutation-20160128: RaggedExperiment with 90490 rows and 993 columns
##  [3] BRCA_RNASeqGene-20160128: SummarizedExperiment with 20502 rows and 878 columns

サンプルデータへアクセス

colData(brca)[1:5, 1:5]
## DataFrame with 5 rows and 5 columns
##                 patientID years_to_birth vital_status days_to_death days_to_last_followup
##               <character>      <integer>    <integer>     <integer>             <integer>
## TCGA-A1-A0SB TCGA-A1-A0SB             70            0            NA                   259
## TCGA-A1-A0SD TCGA-A1-A0SD             59            0            NA                   437
## TCGA-A1-A0SE TCGA-A1-A0SE             56            0            NA                  1321
## TCGA-A1-A0SF TCGA-A1-A0SF             54            0            NA                  1463
## TCGA-A1-A0SG TCGA-A1-A0SG             61            0            NA                   434

サンプルとデータの対応表

colDataに記述されているサンプル名と、各オミクスデータの列名の対応表。一つのサンプルから同じ種類のオミクスデータを複数取得している場合にも対応可能

sampleMap(brca) %>% head()
## DataFrame with 6 rows and 3 columns
##                      assay      primary                      colname
##                   <factor>  <character>                  <character>
## 1 BRCA_RNASeqGene-20160128 TCGA-A1-A0SB TCGA-A1-A0SB-01A-11R-A144-07
## 2 BRCA_RNASeqGene-20160128 TCGA-A1-A0SD TCGA-A1-A0SD-01A-11R-A115-07
## 3 BRCA_RNASeqGene-20160128 TCGA-A1-A0SE TCGA-A1-A0SE-01A-11R-A084-07
## 4 BRCA_RNASeqGene-20160128 TCGA-A1-A0SF TCGA-A1-A0SF-01A-11R-A144-07
## 5 BRCA_RNASeqGene-20160128 TCGA-A1-A0SG TCGA-A1-A0SG-01A-11R-A144-07
## 6 BRCA_RNASeqGene-20160128 TCGA-A1-A0SH TCGA-A1-A0SH-01A-11R-A084-07

特定のサンプルのデータだけを抽出

サンプルTCGA-E9-A54XにはRNA-seqデータがないのでmiRNAと変異データしか出力されていない

brca[ , "TCGA-E9-A54X" , ]
## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes.
##  Containing an ExperimentList class object of length 2:
##  [1] BRCA_miRNASeqGene-20160128: SummarizedExperiment with 1046 rows and 1 columns
##  [2] BRCA_Mutation-20160128: RaggedExperiment with 90490 rows and 1 columns
## Features:
##  experiments() - obtain the ExperimentList instance
##  colData() - the primary/phenotype DataFrame
##  sampleMap() - the sample availability DFrame
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment
##  *Format() - convert into a long or wide DataFrame
##  assays() - convert ExperimentList to a SimpleList of matrices

その他にも便利な関数が多数

# すべてのデータが揃っているサンプルを抽出
complete.cases(brca)
intersectColumns(brca)

# すべてのデータに共通する行 (変数) を抽出
intersectRows(brca)

# 同じサンプル由来のデータをまとめる (デフォルトは平均値でまとめる)
mergeReplicates(brca)

# longフォーマットへの変換
longFormat(brca)

*1:ピークコール後のデータを想定

*2:https://www.bioconductor.org/packages/release/bioc/vignettes/Biobase/inst/doc/ExpressionSetIntroduction.pdf

*3:リードカウントデータ/TPM/RPKM/CPM等