Aurora blog

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

データベース調査①: cBioPortal

cBioPortalとは

がんに由来する大規模・多種類オミクスデータを整理し、誰もが探索的なデータ利用をできるようにしたプラットフォーム。解析に明るく無くてもGUIで簡単に可視化ができる環境が用意されている。

上記の文献が発表された2013年当時は"TCGA + 10 studies"のデータが登録されていたようであるが、2021年現在は300近いデータが公開されている*1。登録されているデータの種類も豊富であり遺伝子変異・遺伝子発現・miRNA発現・Methylomeなどが含まれる。

GUIで出来ることについては(少し古いが)この文献が参考になる。ある遺伝子の特性を色々な癌種で調べたい場合は大体のニーズはGUIで満たせると思う。ここではcBioPortalのデータへプログラムからアクセスする方法についてまとめる。

www.cbioportal.org

R package cBioPortalData

RパッケージcBioPortalDataを使うとcBioPortalに登録されているデータをMultiAssayExperimentクラスのオブジェクトとしてRに直接ロードすることができる。MultiAssayExperimentクラスの使い方については以前のエントリで触れているので良ければご参照ください。

auroratummy.hatenablog.com

cBioPortalDataの使い方はマニュアルに詳細に記載せれている。大まかには以下のように使う。

library(cBioPortalData)

# 利用可能なデータの確認
data("studiesTable", package = "cBioPortalData")
head(studiesTable) 

studiesTableの中身をみるとcBioPortalで扱えるデータの充実度がわかります。どの研究がどの種類のデータを含んでいるかがわからないのは少し不便。ここを参照する必要があるらしい。

DataFrame with 6 rows and 6 columns
      cancer_study_id             study_name            description                    URL
          <character>            <character>            <character>            <character>
1       paac_jhu_2014 Acinar Cell Carcinom.. Whole exome sequenci..                       
2 mel_tsam_liang_2017 Acral Melanoma (TGEN.. Whole exome sequenci..                       
3     all_stjude_2015 Acute Lymphoblastic .. Comprehensive profil..                       
4     all_stjude_2016 Acute Lymphoblastic .. Whole-genome and/or ..                       
5       aml_ohsu_2018 Acute Myeloid Leukem.. Whole-exome sequenci..                       
6           laml_tcga Acute Myeloid Leukem.. TCGA Acute Myeloid L.. http://gdac.broadins..
  pack_build api_build
   <logical> <logical>
1       TRUE      TRUE
2      FALSE      TRUE
3       TRUE      TRUE
4       TRUE      TRUE
5       TRUE      TRUE
6       TRUE      TRUE
# MultiAssayExperimentの読み込み
mae = cBioDataPack("paac_jhu_2014") # studiesTableの`cancer_study_id`列を参照
mae
A MultiAssayExperiment object of 2 listed
 experiments with user-defined names and respective classes.
 Containing an ExperimentList class object of length 2:
 [1] mutations_extended: RaggedExperiment with 2745 rows and 23 columns
 [2] mutations_mskcc: RaggedExperiment with 2745 rows and 23 columns
Functionality:
 experiments() - obtain the ExperimentList instance
 colData() - the primary/phenotype DataFrame
 sampleMap() - the sample coordination DataFrame
 `$`, `[`, `[[` - extract colData columns, subset, or experiment
 *Format() - convert into a long or wide DataFrame
 assays() - convert ExperimentList to a SimpleList of matrices
 exportClass() - save all data to files

マイクロバイオーム解析①:QIIME2を使ったメタ16S解析

QIIME2とは

マイクロバイオーム界隈のバイオインフォマティシャンが共同で開発するメタ16S解析用のツールキット((ちなみにチャイムと読む))。自分が学生の頃は旧バージョンのQIIMEしか存在していなかったように思うが、2019年頃に公開されたらしい。

基本的には、旧バージョンと同じように、コミュニティにおける解析の再現性を高めることや、新規解析手法をQIIMEのプラグインとして公開することでユーザビリティを向上させることが目的である。

QIIME1は得られる恩恵の割に、仕様を覚えるのが面倒で結局使わなかった記憶があるが、下記Nature Biotechnologyの記事によると、旧バージョンには無かった以下のような機能が追加されているらしい。メタ16S解析の最近のデータ解析手法を押さえる意味も込めて、使い方をまとめてみる。

www.nature.com

qiime2.org

解析の流れ

以下のページを参考にQIIME2で紹介されているデータ解析の流れを追う。

docs.qiime2.org

配列データの入力

FastqファイルをQIIME用のフォーマット (.qza) *2 として読み込む。配列データの形式によって若干オプションが変わるが、基本的にはManifestファイルを使った処理をすれば、どのようなファイルに対しても対応可能と思われる。

--input-formatオプションは、Paired-endの場合はPairedEndFastqManifestPhred33V2、Single-endの場合はSingleEndFastqManifestPhred33V2とする。最近はほぼ無いと思うがPhred64の場合は*Phred64V2とする。

qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path ${manifest} \
  --output-path ${prefix}.qza \
  --input-format PairedEndFastqManifestPhred33V2

Manifestファイルの形式については下記リンクを参照のこと。

docs.qiime2.org

Denoising / OTU clustering

OTUクラスタリングとは、類似した配列をOTU (Operational taxonomic unit) と呼ばれるクラスタに分類する操作のことである。慣習として類似性97%が閾値として使われることが多い。この場合、97%以上の互いに一致する配列は一つのOTUへとまとめられる。

# Quality filtering
qiime quality-filter q-score \
  --i-demux ${prefix}.qza \
  --o-filtered-sequences ${prefix}.qc.qza \
  --o-filter-stats ${prefix}.qc_stats.qza

# Dereprecation
qiime vsearch dereplicate-sequences \
  --i-sequences ${prefix}.qc.qza \
  --o-dereplicated-table ${prefix}.qc.derep_table.qza \
  --o-dereplicated-sequences ${prefix}.qc.derep.qza
 
# De novo OTU clustering
qiime vsearch cluster-features-de-novo \
  --i-table ${prefix}.qc.derep_table.qza \
  --i-sequences ${prefix}.qc.derep.qza \
  --p-perc-identity 0.97 \
  --o-clustered-table ${prefix}.otu_97_table.qza \
  --o-clustered-sequences ${prefix}.otu_97.qza

OTUクラスタリングでは、複数の種に由来する配列が一つのOTUとして誤って分類されることがある。そこで最近では、ユニークな配列 (ASV: Amplicon sequence variant) 単位で解析する、DADA2やDeblurなどの手法が提案されている (これらの手法はシーケンスエラーを考慮してユニークな配列を同定する)。 これらの手法にはメリットしかないような気がするが、いまだにOTUベースで解析している論文も多く見る状況。DADA2は低クオリティ配列の前処理を必要としない。以下のコマンドで実行可能。

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs ${prefix}.qza \
  --o-table ${prefix}.asv_table.qza \
  --o-representative-sequences ${prefix}.asv.qza \
  --o-denoising-stats ${prefix}.asv_stats.qza

生物種推定

上記の手法で同定されたASV/OTUの生物種を推定する。QIIMEではBLAST(Local alignment)*3、VSEARCH(Global alignment)、scikit-learn(Naive Bayes classifier)での生物種推定手法を提供している。

学習済みNaive Bayesモデルや、BLAST/VSEARCHに利用可能な参照配列データは以下のリンクでダウンロード可能。SILVAとGreengenes*4のそれぞれで参照データが用意されている。

# BLAST
qiime feature-classifier classify-consensus-blast \
  --i-query ${prefix}.asv.qza \
  --i-reference-reads SILVA.qza \
  --i-reference-taxonomy SILVA.taxonomy.qza \
  --o-classification ${prefix}.asv.blast_SILVA.qza

# VSEARCH
qiime feature-classifier classify-consensus-vsearch \
  --i-query ${prefix}.asv.qza \
  --i-reference-reads SILVA.qza \
  --i-reference-taxonomy SILVA.taxonomy.qza \
  --o-classification ${prefix}.asv.vsearch_SILVA.qza

# scikit-learn
qiime feature-classifier classify-sklearn \
  --i-reads ${prefix}.asv.qza \
  --i-classifier SILVA_sklearn.qza \
  --o-classification ${prefix}.asv.sklearn_SILVA.qza

docs.qiime2.org

系統樹作成

ASV/OTUの配列を使って系統樹を作成する。系統樹はUniFrac距離などの細菌叢のβ多様性を計算する際に使う。系統樹作成には、(i)MAFFT/Clustal等のツールによるマルチプルアラインメントと(ii)Fasttree等のツールによる系統樹推定の2ステップがあるが、QIIME2では一つのコマンドで実行することが可能。

qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences ${prefix}.asv.qza
  --o-alignment ${prefix}.asv.mafft.qza \
  --o-masked-alignment ${prefix}.asv.mafft_masked.qza \
  --o-tree ${prefix}.asv.fasttree.qza \
  --o-rooted-tree ${prefix}.asv.fasttree_rooted.qza

フォーマット変換

QIIME特有のqzaフォーマットを他のツールで利用可能なフォーマットへ変換する。細菌種組成などのテーブルデータはBiom、系統樹データはNewick、配列ファイルはFastaへ変換される。

qiime tools export \
  --input-path feature-table.qza \
  --output-path exported-feature-table

*1:GUIツールは完成版ではないらしい。

*2:QIIME artifactと呼ばれている

*3:あえてLocal alignmentを使う必要性はあるのだろうか

*4:Greengenesは数年前に管理が終わっている

CWL②:BioContainer / scatter / when

CWLを学ぶ part 2

以前のエントリでCWL (Common Workflow Language) について簡単にまとめたが、マニュアルやgithubを見ていると、いくつか知らない機能があることがわかったので、改めてメモを残しておこうと思う。

auroratummy.hatenablog.com

github.com

BioContainer

CWLとは直接は関係がないが、BioContainerには、多種多様なバイオインフォ系のツールをすぐに実行可能なDocker/Singularityイメージが登録されている。マイナーなツールもかなり見つかるので、まずはここを検索してみることを推奨する。

biocontainers.pro

Scatter

複数のサンプルのFASTQファイルを同じリファレンスゲノムにマッピングする場合など、同一作業を繰り返す場合に使われる。*1

ここでは複数のgzipファイルを一つずつ解凍するワークフローを考える。以下のyamlを書くと変数gzに複数ファイルを配列として格納することができる。

# scatter_example.yaml
gz:
  - class: File
    location: test_data_16s_light/ERR3668533_15_L001_R1_001.fastq.gz
  - class: File
    location: test_data_16s_light/ERR3668533_15_L001_R2_001.fastq.gz
  - class: File
    location: test_data_16s_light/ERR3668534_15_L001_R1_001.fastq.gz
  - class: File
    location: test_data_16s_light/ERR3668534_15_L001_R2_001.fastq.gz

stepsの中でscatterを使うことで、複数の入力に対して同一のステップを繰り返し実行することができる。以下の例では、配列gzに含まれる各ファイルに対して、ステップgunzipを繰り返している。この場合、ステップgunzipの出力も配列となるため、このワークフローの出力gunzippedの型にはファイルの配列 (File[]) を指定している。requirementsにScatterFeatureRequirementをつけないとエラーとなるので注意。

# scatter_example.cwl
cwlVersion: v1.2
class: Workflow
doc: scatter test

requirements:
    ScatterFeatureRequirement: {} # これが必要

inputs:
    gz:
        label: gzipped files
        type: File[]

steps:
    gunzip:
        run: components/gunzip.cwl
        scatter: gz # ここ
        in:
            gz: gz
        out: [gunzipped]

outputs:
    gunzipped:
        type: File[]
        outputSource: gunzip/gunzipped

when

CWL version 1.2から追加された機能。これを使うことでif/else的な条件分岐が可能となる。

入力seq_typese(single-end)とpe(paired-end)の場合でデータ処理のパイプラインを変える場合を想定する。この場合はstepにwhenを以下のように挟むことで条件分岐を実現することができる。この場合注意すべきは、ステップpipeline_sepipeline_peseq_typeを入力として受け取る必要がない場合でも、whenを使う場合にはinputにseq_typeを入れなければエラーが出る仕様になっているので注意。*2

# when_example.cwl
cwlVersion: v1.2
class: Workflow
doc: scatter test

requirements:
    InlineJavascriptRequirement: {}
    MultipleInputFeatureRequirement: {}

inputs:
    input_file_dir:
        label: input directory
        type: Directory
    seq_type:
        label: sequencing type ("pe" or "se")
        default: "se"
        type: string

steps:
    pipeline_se:
        in:
            seq_type: seq_type
            input_file_dir: input_file_dir
        when: $(inputs.seq_type == "se")
        run: components/pipeline_se.cwl
        out: [out]

    pipeline_pe:
        in:
            seq_type: seq_type
            input_file_dir: input_file_dir
        when: $(inputs.seq_type == "pe")
        run: components/pipeline_pe.cwl
        out: [out]

outputs:
  result:
    type: File
    outputSource:
      - pipeline_se/out
      - pipeline_pe/out
    pickValue: first_non_null 

ファイルを一箇所にまとめる

ワークフローでたくさんのファイルが出力される場合は、それを特定の名前のディレクトリに集約できるとありがたい。その方法が紹介されているページ*3があったのでメモ。以下のようにExpressionToolで実現できる。

cwlVersion: v1.0
class: ExpressionTool
requirements:
  InlineJavascriptRequirement: {} 

inputs:
    file1:
        type: File
    file2:
        type: File
    file3:
        type: File
    file4:
        type: File

outputs:
    dirs:
        type: Directory[]

expression: |-
    ${
    var dirs = [];
    dirs.push({
        "class": "Directory",
        "basename": inputs.file1.nameroot + "_1",
        "listing": [inputs.file1, inputs.file2]
    });
    dirs.push({
        "class": "Directory",
        "basename": inputs.file1.nameroot + "_2",
        "listing": [inputs.file3, inputs.file4]
    });
    return {"dirs": dirs};
    }

R Markdownで美しいレポートを作成するためのTips

R Markdown

RではR Markdownでコードを記述することで、データ解析の結果、データ解析に使われたプログラム、解析プロセスの説明等をHTMLやPDFなどのフォーマットでレポート化することできる。他人に解析結果と解析方法を共有する場合に大変便利な機能だが、ここではR Markdownで美しいレポートを作るためのTipsをまとめようと思う。

R Markdownの基本的な記法や使い方は以下のページがわかりやすいです。

kazutan.github.io

R Markdownに関する大体の情報は基本的には以下のページに網羅されているので、こちらもご参照ください。

bookdown.org

1. テンプレートを使う

rmdformatsprettydocといったパッケージが美しいデザインのテンプレートを用意しておりオススメ。rmdformatsだけでも>5種類のレイアウトが用意されている。個人的にはrobobookが気に入っている。

github.com

基本的にはR Markdownのヘッダーのoutputに以下のような内容を加えるだけでOK。個人的によく使うオプションは以下。

  • code_folding: デフォルトでスクリプト部分を隠したい (ボタンを押すと表示される) 場合は"hide"とする
  • number_sections: セクション番号 (e.g. 1. 〇〇〇〇〇) をつける場合はtrue
  • thumbnails: 画像をサムネイルとして表示する場合はtrue (大きい画像が存在する場合におすすめ)
---
title: "DOCUMENT TITLE"
date: "`r Sys.Date()`"
author: auroratummy
output:
  rmdformats::robobook:
    code_folding: hide
    number_sections: false
    thumbnails: true
---

2. テーブル

テーブルはkableExtraDT、formattableなどのパッケージを使うことで通常よりもスマートに表示することが可能です。formattableを使うとより高度な表の可視化が簡単にできるので、こちらも使ってみても良いかもしれません。*1

data(iris)
library(kableExtra)
kbl(iris) %>% 
  kable_styling(bootstrap_options = "stripe") %>%
  scroll_box(width = "500px", height = "200px")

f:id:auroratummy:20210314215650p:plain

3. タブを使う

複数の条件で同様の解析を繰り返す場合などは、何も工夫せずにR Markdownでまとめると、たくさんの似た解析結果・コードが繰り返された冗長なレポートが作成されてしまう。このような場合は{.tabset}を利用すると以下のようにタブで複数の解析結果・コードを一箇所にまとめることができる。

{.tabset}の使い方は以下の通り。レベル1 ("#") のヘッダに{.tabset}を指定すると、レベル2 ("##") のヘッダ下に記載される内容が一つのタブとしてまとめられる。レベル2 ("##") のヘッダに{.tabset}を使うと、レベル3 ("###") の内容がまとめられる。

f:id:auroratummy:20210314220450p:plain

f:id:auroratummy:20210314220207p:plain

4. カラーパレット

以下のライブラリが便利。個人的には、Nature系雑誌で使われるカラーパレットなど、デザイン性の高いパレットを提供してくれるggsciが好み。ggsciが提供する関数はggplot2にそのまま利用できるので、この点も使いやすい。

fig = ggplot(iris) +
  geom_point(aes(x=Sepal.Length, y=Sepal.Width, color=Species)) +
  theme_bw() + 
  ggsci::scale_color_nejm()
fig

f:id:auroratummy:20210314222532p:plain

5. plotlyを使う

可視化の際には基本的にplotlyを使うようにしている。特に一つの図の中に色の数がたくさんある場合は便利。ggplot2で作成した図のオブジェクトをggplotly関数に与えるだけで簡単にplotly的なインタラクティブな図にできる。

fig = ggplot(iris) +
  geom_point(aes(x=Sepal.Length, y=Sepal.Width, color=Species)) +
  theme_bw() + 
  ggsci::scale_color_nejm()
library(plotly)
ggplotly(fig)

f:id:auroratummy:20210314222337p:plain

空間オミクス解析③:Scanpy/Squidpyを使ってPythonで空間オミクスデータを扱う

Squidpyとは

Squidpyは、シングルセルオミクスデータの探索的データ解析(EDA)に使われるScanpyを開発したFabian Theisのグループがつい最近公開した空間オミクス解析のためのPythonモジュール。

squidpy.readthedocs.io

www.biorxiv.org

以前のエントリで空間オミクスデータのEDAに使えるツールを紹介したが、そのほとんどがRで実装されたものであり、Pythonで実装された目ぼしいツールは見つからなかった。空間オミクスデータを使って病理画像から遺伝子発現を予測する手法*1や、Visiumの空間オミクスデータの解像度を病理画像を使って向上させる手法*2など、機械学習をうまく利用した研究が現れていることからも、Pythonでの開発の基盤が欲しいと感じていたが、今回それがようやく現れたと言える。

auroratummy.hatenablog.com

データの取り扱い

SquidpyはScanpyと同様にAnnDataクラスでデータを扱う。anndataはRのSeuratクラスやSingleCellDataクラスと同じようにシングルセル データのために設計されたデータ型である。遺伝子発現行列の他に、サンプル (細胞/スポット) のメタデータ、計測した遺伝子のメタデータ、次元圧縮後の座標情報やkNNグラフなどの解析結果を一つのオブジェクトで取り扱うことができる。

空間オミクスデータのシングルセルデータとの大きな違いはHE画像などの病理画像情報が付随している場合があることだが、これはAnnDataのadata.uns['spatial']['{library_id}']adata.obsm['spatial']に格納される。VisiumデータをSpace Rangerで解析した場合は、遺伝子発現情報が格納されたhdf5ファイル *3 と画像データや各スポットの座標情報が格納されたspatialフォルダがあれば、Scanpyのread_visiumで画像もまとめて読み込むことができる。

!tree 
├── data
│   ├── filtered_feature_bc_matrix.h5
│   └── spatial
│       ├── aligned_fiducials.jpg
│       ├── detected_tissue_image.jpg
│       ├── scalefactors_json.json
│       ├── tissue_hires_image.png
│       ├── tissue_lowres_image.png
│       └── tissue_positions_list.csv
└── squidpy.ipynb
import scanpy as sc
adata = sc.read_visium("./data/", count_file="filtered_feature_bc_matrix.h5")
adata
AnnData object with n_obs × n_vars = 4992 × 36601
    obs: 'in_tissue', 'array_row', 'array_col'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatial'
    obsm: 'spatial'

adata.uns['spatial']['{library_id}']spatialフォルダに含まれる画像や各スポットの座標情報が辞書型で格納される。

adata.uns["spatial"]["Parent_Visium_Human_Cerebellum"].keys()
dict_keys(['images', 'scalefactors', 'metadata'])

`adata.uns['spatial']['{library_id}']["images"]には2つの解像度 (hi & low resolution) の画像データがarray型で格納される。

# adata.uns["spatial"]["Parent_Visium_Human_Cerebellum"]["images"]["hires"]
adata.uns["spatial"]["Parent_Visium_Human_Cerebellum"]["images"]["lowres"][:2, :2, :] # (x, y, RGB)
array([[[0.7019608, 0.7176471, 0.7176471],
        [0.7019608, 0.7137255, 0.7176471]],

       [[0.7019608, 0.7137255, 0.7176471],
        [0.7019608, 0.7137255, 0.7176471]]], dtype=float32)

adata.obsm['spatial']には各スポットの座標情報が格納されている。

adata.obsm['spatial'][:5, :] # (x, y)
array([[3366, 1800],
       [9312, 7772],
       [5227, 2153],
       [3598, 8871],
       [8745, 3459]])

adata.obsm['spatial']は圧縮前の画像データ ((フォルダspatialの画像データは高解像度(hires)のものも元の画像からは圧縮されている)) 上の座標 (ピクセル) データであるため、hires & lowresいずれの画像上に直接投影することができない。adata.uns['spatial']['{library_id}']["scalefactors"]には画像上にスポットを投影するための情報が格納されている。

tissue_hires_scaleftissue_lowres_scalefにはhires/lowres両画像が元の画像からどの程度縮小されたものであるか(=縮尺)が格納される。spot_diameter_fullresは元画像におけるスポットのサイズ (ピクセル数) を表す。これらの情報をもとにスポットの場所をhires/lowres画像上へ投影することができる。

adata.uns["spatial"]["Parent_Visium_Human_Cerebellum"]["scalefactors"]
{'spot_diameter_fullres': 89.55377,
 'tissue_hires_scalef': 0.150015,
 'fiducial_diameter_fullres': 144.66379000000003,
 'tissue_lowres_scalef': 0.045004502}

Scanpy/SquidpyはImaging mass cytometryやFISHベースの手法で得られた空間発現情報も扱うことができる、という触れ込みであるが、sc.read_visium的な関数はVisium以外のデータにはまだ用意されておらず、現状は手作業でadata.uns["spatial"]adata.obsm["spatial"]を用意する必要がありそうだ。*4

出来ること

ImageContainer

SquidpyではImageContainerクラスという画像情報を扱うクラスを新たに提供している。ImageContainerは画像をxarray *5 やDaskで扱いやすい&メモリ効率的な形で保存する。SquidpyはImageContainerを入力とした様々な画像解析の機能を提供している。*6

ImageContainerの作り方はとても簡単。

image = adata.uns["spatial"]["Parent_Visium_Human_Cerebellum"]["images"]["hires"]
img = sq.im.ImageContainer(image)
img
ImageContainer object with 1 layer:
    image: y (2000), x (2000), channels (3)

空間情報を利用した解析

anndataクラスでデータを扱うため、Scanpyが提供するシングルセル解析用の機能 (e.g. Embedding, Graph-based clustering, DEG analysis) はすべてそのまま適用することができる。Squidpyではこれに加えて空間オミクスデータならではの以下の解析手法を提供している。

  • 近い場所に位置している細胞種の抽出 Link
  • 特定の領域で発現している遺伝子の抽出:自己相関の指標 Moran's Iを利用 Link
  • リガンドーレセプター解析 (CellphoneDBの情報を利用): 特定のリガンド・レセプターのペアを高発現する空間を探索

外部ツールとの連携

Napariという画像用のビューアーと連携している。AnnDataとImageContentを指定して以下のコマンドを実施するとGUIが立ち上がりインタラクティブに遺伝子発現の可視化を行うことができる。Napariは可視化だけでなく、病理学で行われるようなマニュアルで画像にアノテーションをつける事もできる。マニュアルアノテーションの内容はadata.obsへ保存されるため、プログラムでの解析へシームレスにつなげることができる。Napariで作業すると自分のPC ((メモリ16GB, 2.3GHz Quad core Intel Core i7Mac Book Pro)) の動作がだいぶ遅くなったので、ある程度のメモリは必要と思われる。

%gui qt
viewer = img.interactive(adata)

f:id:auroratummy:20210227013112p:plain

そのほかにもTensorflow/Pytorchといった深層学習用のライブラリとの連携や、Tanglamやstlearnといった空間オミクスとシングルセルデータの統合解析手法もScanpy/Squidpyからそのまま使えるようにしていたり、Rで実装されているEDA用のツールと比較してもかなりよく出来ている印象。*7

*1:https://www.nature.com/articles/s41551-020-0578-x

*2:https://www.biorxiv.org/content/10.1101/2020.02.28.963413v2.full

*3:filtered_feature_bc_matrix.h5 or raw_feature_bc_matrix.h5

*4:一般的なフォーマットが存在しない現状は関数も開発できないのであろう

*5:Pandasよりも多次元配列を扱いやすいらしい

*6:最初からadata.unsの画像データをImageContainer型にしないのは何故だろうか

*7:機械学習を使った研究への展開が容易になりそうで、とてもありがたい

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等

空間オミクス解析②:解析用ライブラリについて

空間オミクス解析のライブラリ

シングルセル解析が普及した背景はChromiumの様な商用の解析装置がリリースされたことに加えて、SeuratScanpyなどの解析用ライブラリが整備されたことが大きい。これらのライブラリに付属するマニュアルを追いかけさえすれば、複雑なプログラムを実装せずに、主要な解析はほとんど実行することができる。

シングルセル解析が一般的なものになっていく中で、いま成熟し始めているのは前回のエントリーで触れた空間オミクス解析だと思う。Human Cell AtlasHuman Tumor Atlas Networkで位置情報付きのオミクスデータが取得されていることを鑑みると、今後身近なデータになっていくことは間違いない。

auroratummy.hatenablog.com

そこで空間オミクスの解析用ライブラリについて調べてみるといくつか見つかった。*1

Seurat/Scanpy

シングルセル解析でお馴染みのSeuratとScanpyであるが空間オミクス解析のための機能もいくつか実装されている。両ツールでできることはほぼ同じ。シングルセル解析の延長上で空間オミクスデータを扱うことができ、データ前処理・次元圧縮・クラスタリングにはscRNA-seqで使われている手法・関数を使用する。発現解析にはSpatialDEやTrendseekの様な特定の空間で発現上昇・低下している遺伝子を抽出する手法が実装されている。基本的にはVisiumのデータに向けて開発が進められているらしく、その他の手法で得られた空間オミクスデータに対してはまだ機能が不十分な印象を受ける。

  • データ前処理(QC・正規化)
  • 次元圧縮・クラスタリング
  • 空間的発現変動遺伝子の抽出
  • 可視化(病理画像上への発現情報の投影)

satijalab.org

scanpy-tutorials.readthedocs.io

STUtility

Spatial transcriptome (ST) の開発者であるJoakim Lundebergらのグループが開発した空間オミクス解析用のライブラリ。VisiumやSTのデータに向けて開発されている。複数の切片を組み合わせて3Dデータを作る際に画像間のアラインメントをするための機能や、画像のアノテーションインタラクティブに行うための機能が実装されている。

f:id:auroratummy:20201214224850p:plain

引用元:Bergenstrahle et al. “Seamless integration of image and molecular analysis for spatial transcriptomics workflows” BMC Genomics. 2020. doi:10.1186/s12864-020-06832-3 CC-BY 4.0

bmcgenomics.biomedcentral.com

Giotto

MERFISH等の開発に関わるGuo-Chen Yuanらのグループが開発した空間オミクス解析用のRライブラリ。Visium以外の様々な種類の空間オミクスデータや、2Dだけでなく3Dの空間データ*2にも対応していることが特徴。

  • データ前処理(QC・正規化)

  • 次元圧縮・クラスタリング

  • 空間的発現変動遺伝子の抽出
  • 可視化(病理画像上への発現情報の投影)
  • 共発現解析(同じ空間で発現する傾向にある遺伝子モジュールの同定)
  • 細胞間相互作用解析(e.g. 近傍に存在する細胞種のペアを探索、周辺にいる細胞が発現に与える影響を解析)
f:id:auroratummy:20201214203541p:plain

引用元:Dries et al. “Giotto, a toolbox for integrative analysis and visualization of spatial expression data” bioRxiv. 2020. doi:10.1101/701680 CC-BY 4.0

www.biorxiv.org

Giotto viewerというインタラクティブな可視化ツールも合わせて公開している。特に3Dデータの可視化は難しいので、このようなツールが用意されていることもユーザーとしてはありがたい。

www.youtube.com

Giottoの機能を著者らが説明したウェビナーの動画も公開されている。

www.youtube.com

まとめ

シングルセル解析におけるSeurat/Scanpyのようなツールが空間オミクス解析の領域でも現れつつある印象を受けた。空間オミクス解析は機械学習と相性が良いと思うのでPythonで扱えるライブラリも欲しい。

*1:他にいいものがあったら誰か教えてください

*2:STARmapのような3Dで得られるデータや、複数の切片のデータを統合した場合