はじめに
最近一つの実験で複数種類の高次元データを取得する機会や、ATAC-seqやChIP-seqなどのパラメータがゲノム領域に対応するデータ*1を扱う機会が増えてきた。これまで単一のデータをRで解析する際にはdata.frameクラスやmatrixクラスでデータを扱ってきたが、不便さを感じる場面が出てきたので、少しずつ別のクラスへと移行していこうと考えている。R/BioconductorではSummarizedExperimentやMultiAssayExperimentなどの便利なクラスが提供されている。ここではその機能についてまとめようと思う。
SummarizedExperiment
SummarizedExperimentクラスは下図のような構造を有する。遺伝子発現などの行列データに加えて、行 (遺伝子名やピーク) と列 (サンプル) の情報、実験背景等のメタデータを一つのオブジェクトでまとめて扱うことができる。
Biobaseライブラリが提供するExpressionSetクラス*2も類似した機能を有し、広く使われているが、パラメータの情報をGRangesクラスで格納できるためATAC-seq等のデータへ親和性が高いことや、異なる方法で前処理した行列データ *3 を一つのオブジェクトへまとめて格納できるなどの優位性がある。個人的には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
MultiAssayExperiment
SummarizedExperimentクラスは複数の行列データを格納することが可能だが、各データの行数・列数が等しくなければならないという制約がある。当然だがオミクスデータは種類によって次元数は異なり、パブリックのマルチオミクスデータは一部のデータが一部のサンプルで未計測であることも多く、多種類のデータを扱う場合はSummarizedExperimentクラスでは対応できないことが多い。MultiAssayExperimentは多種類データに向けてデザインされたデータ型である。
TCGAやcBioPortalのデータをMultiAssayExperimentクラスで利用可能な形でダウンロードできるRパッケージcuratedTCGAData
cBioPortalData
も公開されている。
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)