2020-10-28

シングルセル解析③：ダブレット推定手法について

Omics Bioinfo Single-cell

ダブレットとは

シングルセル解析で現在主流になりつつあるのは10x Genomics Chromiumのようなドロップレットベースの方法である。この方法はドロップレットの中で磁気ビーズと細胞を封入し、mRNAの逆転写・増幅を行う (下図b)。ビーズにはポリT領域と、ビーズごとにユニークなバーコード配列がついたオリゴDNA (図c) が生えており、この上で細胞由来のmRNAが逆転写・増幅される仕組みである。この時にバーコードがついた状態でmRNAは逆転写・増幅されるので、配列解析時に、同じバーコードを持つリードは同一細胞に由来すると判別することができる。

引用元：Zheng et al. “Massively parallel digital transcriptional profiling of single cells” Nat Commun. 2017 Jan 16;8:14049. doi: 10.1038/ncomms14049. CC-BY 4.0

一つのドロップレットの中にどのように一つの細胞を閉じ込めているのか？これは機械的に制御されているわけではない。ドロップレットの中に封入される細胞の数は確率的に決まり、ポアソン分布に従う。例えば、一つのドロップレットあたり平均0.5個の細胞が封入されるような濃度で実験を行えば、確率的に約60%, 30%, 8%のドロップレットには0個, 1個, 2個の細胞が封入される *1。すなわち細胞の濃度をある程度下げることで、ほとんどの細胞は単独でドロップレットの中に封入されることとなり、問題なくシングルセルレベルでの計測をすることができる、という算段である。

しかし実際には数%のドロップレットには複数個 (多くは2個) の細胞が入ってしまう (下図A)。この場合は、同じビーズ上でmRNAが処理され、同一のバーコードが付けられてしまう。したがって、配列解析時に異なる細胞由来のmRNAであることが区別できないため、複数個の細胞の遺伝子発現が、一つの細胞として、混ざった状態で検出されることになる。これがいわゆるダブレットと呼ばれるアーティファクトである。

ダブレットは解析結果を歪めかねない。一見するとユニークな細胞集団に見えて、実は単なるアーティフィシャルであることも多い (下図B)。

引用元：Wolock et al. “Scrublet: computational identification of cell doublets in single-cell transcriptomic data” bioRxiv. 2018 https://doi.org/10.1101/357368 CC-BY 4.0

Chromiumのマニュアル*2を見ると、１サンプルあたり10000細胞を解析する条件では、約7.6%がダブレットになるようだ。これは結構大きな影響がありそうなものだが、多くの論文ではあまりケアされていないように感じる。ダブレットを除くために、検出遺伝子数 (もしくはUMI数) が極端に多い細胞を解析から除外する前処理をよく目にするが、これは理想的とは言えない。検出遺伝子数は細胞種間で不均一であり、活性化した細胞や癌細胞は多くの遺伝子を発現しているからだ。

そんな中で2019年ごろから、よりエレガントにダブレットを解析から除外するための方法論が提案されてきた。ここではそれらの方法についてレビューする。

Scrublet, DoubletFinder

最初に報告された2つの手法 Scrublet*3とDoubletFinder*4。2019年4月の同時期にCell Systemsに公開されたこれらの手法はいずれも似たアプローチを提案している。基本的には以下の流れ。

入力: UMIカウント行列 (細胞 x 遺伝子)
1からランダムに細胞ペアをN組選出
2のペアのカウントデータをミックス仮想ダブレットデータを作成
3を1と混ぜる
Highly-variable geneの抽出 & PCAで次元圧縮 *5
主成分得点を使ってkNN graph *6 を描画
6の情報をもとにスコアを計算

Scrubletは以下の式でスコアが求まる*7。ダブレットはこのスコアが高くなる。

$score_{i}=\frac{q_{i} \hat{\rho} / r}{ 1 - \hat{\rho} - q_{i}(1 - \hat{\rho} - \hat{\rho} / r) }$
$q_{i} = \frac{k_{d}(i) + 1}{k_{adj} + 2}$
$\hat{\rho}$ : 予想されるダブレットの割合
$r = M_{d}/M_{obs}$
$k_{adj}$ : k
$k_{d}(i)$ : 細胞iのkNNに含まれる仮想ダブレットの数
$M_{obs}$ : 全細胞数
$M_{d}$ : 作成した仮想ダブレットの数

$\hat{\rho}$ はユーザーが解析の際に指定する。（このパラメータが結果に影響することがgithubに記載されていた）

シングルセルデータにはダブレットとシングレットが混ざっているのでスコアの分布は二峰性を示す。この分布をもとにユーザーがカットオフをマニュアルで設定してダブレットを特定することが推奨されている。

github.com

DoubletDecon

2019年に報告された手法。Scrublet/DoubletFinderが $\hat{\rho}$ などのパラメータをユーザーが指定しなければならないことを問題視 *8。仮想ダブレットを作成するプロセスは共通しているが、発現変動遺伝子の情報を使って、過度に細胞を除きすぎないようにするプロセス (Rescue step) を設けている点が特徴。

入力: 正規化済み発現行列 (細胞 x 遺伝子)
前処理 (低発現遺伝子除外 & 質の低い細胞を除外)
クラスタリング *9
3のクラスタごとに代表発現プロファイル (Centroid) を計算
Cluster merging: 類似したクラスタをマージする *10
仮想ダブレットを作成: 5でできたクラスタのペアごとにランダムに細胞を選んで仮想ダブレットを作成
4をリファレンスに各細胞/仮想ダブレットの遺伝子発現プロファイルをDeconvolute。下図Eのような組成を計算する
7で得られた組成が仮想ダブレットの組成と類似している細胞は除外
8で除外した細胞を対象に3と同様の方法でクラスタリング
3と9の全クラスタ間で発現変動遺伝子を抽出
Rescue step: 特異的遺伝子が検出されたクラスタは9からレスキューする

引用元：DePasquale et al. “DoubletDecon: Cell-State Aware Removal of Single-Cell RNA-Seq Doublets” bioRxiv. 2019 doi: https://doi.org/10.1101/364810 CC-BY 4.0

イントロではScrublet/DoubletFinderのパラメータ設定を問題視していたが、DoubletDeconも5, などのステップでad hocにユーザーがパラメータを決めなければならない様子。 Rescue stepは2つの細胞が混ざってできたダブレットであれば特異的に発現が動いている遺伝子はないという仮定に基づいている。この考え方は実装も楽だし、普遍的に使えそう。

github.com

SCDS

2020年に公開された手法。Co-expression based doublet scoring (cxds)とBinary classification based doublet scoring (bcds)という2つのスコアをもとにダブレットを推定する。

Co-expression based doublet scoring (cxds)

ある遺伝子 $i$ が検出される細胞の割合を $p_i$ とする。仮に2つの遺伝子 $i$ と $j$ の発現が独立している場合は、2つのうちどちらかが検出される細胞の数 $n_{ij}$ は次の二項分布に従う。これを利用すると2つの遺伝子の発現の排他性 (互いにExclusiveか) を検定することができる。

$n_{ij} \sim Bin(n, p_{i}(1 - p_{j}) + p_{j}(1 - p_{i}))$

cxdsは以下のように求められる。通常は排他的に発現しているはずの遺伝子ペアが同時に発現しているほどスコアが高くなる仕組み。

$cxds(i) = \sum_{k}\sum_{j} B_{ki}B_{ji}S_{ij}$
$B_{ki}$ : 細胞iにおいて遺伝子kが検出されていれば1, そうでなければ0のバイナリ行列
$S_{ij}$ : 上記の二項検定で得られたp-valueを-log10変換したもの (遺伝子iとjの排他性)

Binary classification based doublet scoring (bxds)

入力: UMIカウント行列 (細胞 x 遺伝子)
1からM個のHighly-variable geneを選出 (default: M=500)
2からランダムに細胞ペアをN組選出
3のペアのカウントデータをミックス仮想ダブレットデータを作成
xgboostで仮想ダブレットと実際の細胞を分類する学習機を作る
仮想ダブレットへの分類確率をスコア (bxds) とする

github.com

まとめ

今回はダブレット推定をベースにした解析のアイデアがあったので世に出てきている手法を調べてみた*11。仮想ダブレットを作ってkNNやxgboost等で判別する手法が多い印象 *12 。性能や使いやすさについては追って検証したら追記する。

*1:λ=0.5のポアソン分布

*2:https://kb.10xgenomics.com/hc/en-us/articles/360001378811-What-is-the-maximum-number-of-cells-that-can-be-profiled-

*3:https://www.sciencedirect.com/science/article/pii/S2405471218304745

*4:https://www.sciencedirect.com/science/article/pii/S2405471219300730

*5:Scrubletは1のみでPCA & その主成分へ3を後から射影

*6:k-nearest neighbor: 最近傍k個の細胞を結んだグラフ

*7:この式が導出された理論も論文に記載されていたがいまいち理解できず

*8:パブリックデータの場合はわからない場合も多い

*9:論文ではICGSを使用

*10:Batch effectなどで同一細胞が異なるクラスタとして検出される可能性を考慮しているらしい

*11:ダブレット推定問題自体にはあまり興味はない

*12:DNNで学習する手法 Soloも最近提案されてた https://www.sciencedirect.com/science/article/pii/S2405471220301952

2020-10-14

シングルセル解析②：CITEseqデータの正規化手法について

Omics Bioinfo Single-cell

CITEseqとは

シングルセルトランスクリプトームと表面タンパク質発現量*1を同一細胞から計測する技術*2。概要としては以下のような手法である。

オリゴDNAを付与した抗体を使って細胞を標識。このオリゴDNAにはポリA配列*3、標的タンパク毎にに異なるバーコード配列が含まれる
ドロップレットの中に磁気ビーズ*4と細胞を封入
磁気ビーズ上でオリゴDNAとmRNAを増幅
シーケンシング
配列がどの抗体・細胞に由来していたかはバーコード情報から判別可能。これをもとに表面タンパク発現の定量を行う

引用元：Matula et al. 2020. “ Single-Cell Analysis Using Droplet Microfluidics.” Adv Biosyst. 2020 Jan;4(1):e1900188. doi: 10.1002/adbi.201900188.. CC-BY 4.0

10x Genomicsよりキットがリリースされており*5、CITEseq用のオリゴDNAを付与した抗体もBiolegend社から販売されている*6ため、少しづつではあるが報告例が増えてきている。また今年はIdo Amitのグループが細胞内タンパク質を抗体標識後、FACSで任意の細胞を分取しscRNA-seq解析を行う手法INs-seqを提案している*7。本手法は同時標識できるタンパク質の数は限定されているが、今後CITEseqのようにDNAバーコード技術を利用することで、より多くの細胞内タンパク質をシングルセルレベルで同時計測する技術が現れるものと予想される。

CITEseqデータの正規化

scRNA-seqデータの正規化手法は複数提案されているが、多くの研究では、Smart-seq2のようなFull-lengthの手法であればTPM (Transcript per million)、ChromiumのようなUMIベースの手法であればCP10k (Count per 10 thousands) などのシンプルな正規化手法が使われていると思う。これらの正規化手法が使われる背景には「一つの細胞内におけるmRNAの総量はほぼ等しい」という仮定が存在する。

一方でCITEseqはプロテオミクスのような網羅的計測でなく、せいぜい百種類の表面タンパク質のみが対象となる。この場合は「計測した表面タンパク質の総量が細胞間でほぼ等しい」という仮定は明らかに現実と乖離していることが予想されるため、TPMやCP10kのような合計値でのシンプルな正規化は行うことができない。このような場合にはどのように正規化をするのだろうか？色々と調べているうちに疑問に思ったので以下にメモを残す。

CITEseq 原著論文*8の場合

StoeckiusらのCITEseq原著論文ではCLR (Centered log-ratio) transformationを適用している。細胞i・タンパク質jのリードカウント数をXijとすると、変換後の値Yijは以下のようになる。

$Y_{i.} = clr(X_{i.}) = (ln(\frac{X_{i1}}{g(X_{i.})}), ln(\frac{X_{i2}}{g(X_{i.})}), ... , ln(\frac{X_{in}}{g(X_{i.})}))$
$g(X_{i.})$ ： $X_{i.}$ の幾何平均

この変換は組成データ(e.g. パーセント) によく用いられる。組成データが持つ、和が一定の値になる制約 (定数和制約) による不都合*9を解消するために使われるらしい。このあたりの話は太田らの資料が大変参考になった。

マテメソを読むと"we treated this data type as compositional data"と書いてある。結局「計測した表面タンパク質の総量が細胞間でほぼ等しい」という仮定を置いているということだと思うが、そのロジックは理解できない。明らかに間違った値へと変換される場合もあると思う*10。

Seuratの場合

Seuratでは関数NormalizeDataにCITEseqデータ用の正規化手法も実装されている。Vignetteによると、ここでもCLRが推奨されているが、"This is a slightly improved procedure from the original publication"と若干の改良がなされたことも記載されている*11。どのような改良がされているのだろうか？見当たらないのでソースコードを見てみた。

github.com

2326行目〜

      'CLR' = CustomNormalize(
        data = object,
        custom_function = function(x) {
          return(log1p(x = x / (exp(x = sum(log1p(x = x[x > 0]), na.rm = TRUE) / length(x = x)))))
        },
        margin = margin, # default margin = 1
        verbose = verbose
        # across = across
      )

3173行目〜

CustomNormalize <- function(data, custom_function, margin, verbose = TRUE) {
  # margin <- switch(
  #   EXPR = across,
  #   'cells' = 2,
  #   'features' = 1,
  #   stop("'across' must be either 'cells' or 'features'")
  # )
  norm.data <- myapply(
    X = data,
    MARGIN = margin,
    FUN = custom_function)
  if (margin == 1) {
    norm.data = t(x = norm.data)
  }
}

該当する箇所を上に示しました。どうやらデフォルトではFeatureごとにlog1p(x = x / (exp(x = sum(log1p(x = x[x > 0]), na.rm = TRUE) / length(x = x))))という変換を加えているようだ。

原著論文では細胞ごとに全タンパクの幾何平均を求めてCLR変換をしていたわけだが、ここではタンパク質ごとに全細胞の幾何平均を求めてCLR変換をしていた。根拠は不明。そしてlog1p (log(x+1))を使っているがこれも謎。同じところを気になった人はいるようだがSeuratチームの回答はありませんでした*12。

DSB

新しい手法DSBが提案されたようなのでチェック。手法の論文はまだbioRxivでしか公開されていないが、すでにDSBを使ってCITEseqデータを解析した事例がNature Medicineに出ていた*13。

mattpm.github.io

本手法はバックグラウンド(空のドロップレット)の情報を使う。確かにこれであればCLRのような組成データの仮定がないので問題が少なそう。

$Y_{ij} = \frac{log(X_{ij} + P) - \mu}{\sigma}$
$P$ ：pseudocount to stabilize variance of small values (default 10)
$\mu$ ：バックグラウンドにおけるタンパクjの平均リードカウント数
$\sigma$ ：バックグラウンドにおけるタンパクjのカウント数のSD

結論

CITEseqデータの正規化手法はまだ発展途上という印象。CLRベースの手法はスタンダードとして使われているが、結局のところ問題は多そう。内部標準を設けることが難しいことを考えるとこれが限界なのだろうか？

*1:<100種類程度のタンパク質を同時計測可能

*2:https://www.nature.com/articles/nmeth.4380

*3:TotalSeqBは別のCapture sequenceが付与されている

*4:ビーズごとに異なるバーコード配列を持つ

*5:https://www.10xgenomics.com/products/single-cell-immune-profiling

*6:https://www.biolegend.com/en-us/totalseq

*7:https://www.sciencedirect.com/science/article/abs/pii/S0092867420308096

*8:https://www.nature.com/articles/nmeth.4380

*9:例えばパラメータ間の相関を見る場合組成データはパラメータ間が独立でないので(一方が増えると他方が減るので)相関が議論できない

*10:例えば、解析対象の表面タンパク質をほとんど発現していない細胞があるならば、どうなるだろうか

*11:will release more advanced version, と書いてあるのでSatijaらもCLRには問題を感じているのだろう

*12:https://github.com/satijalab/seurat/issues/1268

*13:https://www.nature.com/articles/s41591-020-0769-8

2020-10-12

Local BLASTの結果からLCA (Lowest common ancestor)で由来生物種を推定する

きっかけ

最近仕事でLocal BLASTの結果をパースして配列が由来する生物種を推定する機会があったのだが、なかなか思うような機能を実装したツールが存在していなかった*1ので作ってみた。

準備

以下のインストール・ダウンロードを行う

NCBI BLAST+
BLAST database files
- 以下のDBはBLAST用のフォーマットに変換されたものがNCBIによって提供されているのでそれをダウンロード
  - RefSeq protein
  - NCBI nr
Python 3
- sqlite

使い方

スクリプト

https://github.com/tmaruy/blast2taxonomy

1. SQLiteでDBを作成する

git clone https://github.com/tmaruy/blast2taxonomy.git
cd blast2taxonomy
chmod +x preprare_taxdb.sh
./prepare_taxdb.sh

上記のコマンドで次の2つのテーブルを含むDBファイル (taxonomy.db) が作られる。トータルで30GB強のサイズとなるので注意。

accession2taxid: NCBIのアクセッション番号とTaxonomy IDの対応表
- ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/から作成
taxonomy: Taxonomy IDと生物種名の対応表
- ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zipから作成

accession2taxid

SELECT * FROM accession2taxid LIMIT 5;

accession|taxid
A0A009IHW8|470    
A0A023GPI8|232300 
A0A023GPJ0|716541 
A0A023GS28|37565  
A0A023GS29|37565

taxonomy

SELECT * FROM taxonomy LIMIT 5;

taxid|category|name|t_domain|t_phylum|t_class|t_order|t_family|t_genus|t_species
1|no rank|root|||||||
2|superkingdom|Bacteria|Bacteria||||||
6|genus|Azorhizobium|Bacteria|Proteobacteria|Alphaproteobacteria|Rhizobiales|Xanthobacteraceae|Azorhizobium|
7|species|Azorhizobium caulinodans|Bacteria|Proteobacteria|Alphaproteobacteria|Rhizobiales|Xanthobacteraceae|Azorhizobium|Azorhizobium caulinodans
9|species|Buchnera aphidicola|Bacteria|Proteobacteria|Gammaproteobacteria|Enterobacterales|Erwiniaceae|Buchnera|Buchnera aphidicola

2. BLAST+

タブ形式で結果を出力 (-outfmt 6)

blastx -db /path/to/refseq_protein \
       -query test/input.fasta \
       -outfmt 6 \
       -out test/blast.txt

3. 由来生物種推定

-i : Input file, result of blast search (made at step 2)
-db : DB file taxonomy.db (made at step 1)
-o : Output file
-evalue : Threshold of evalue (Ignore hits if their evalues are above this threshold)
-lca : Top N hits used for calculation of lowest common ancestors (default N=10)
-lca_threshold : Returns classification if shared by X % of top hits (default X=50)

#git clone https://github.com/tmaruy/blast2taxonomy.git
python blast2taxonomy_lca.py -i test/blast.txt \
                             -db taxonomy.db \
                             -o test/tax_lca.txt \
                             -evalue 0.01 \
                             -lca 10 \
                             -lca_threshold 50

*1:すでにサポートが終了したものは見つかった https://github.com/frederikseersholm/blast_getLCA

2020-09-26

シングルセル解析①: reticulateでRとPythonを両方使う

Bioinfo Omics Single-cell

RとPythonどちらを使う？

10x GenomicsのChromiumが販売され始めてからシングルセル解析は身近なツールになり多くの研究者に使われるようになっています*1。

それに伴い解析のためのツールも続々と開発されており、誰でもデータを利用しようと思えばできる環境がすでにできつつあります*2。そんな贅沢な状況の中で少し厄介なのは、開発に使われる言語がRとPythonで二分されていることです。

基本的な解析*3は、RではSeurat、PythonではScanpyですべて行うことが出来ます。CITEseq等のマルチオミクスデータへの対応はSeuratの方が優れている一方で、Scanpyで利用できるよう実装されている解析手法もあり *4、どちらもうまく活用したいのが本音です。

また、Trajectory inference・バッチ補正・ダブレット推定・細胞種アノテーションなどの工程は独立したR/Pythonライブラリとして公開されていることが多いです。そして、これらのライブラリを試そうと思う時に実装されている言語が自分が使っているものと違うと、少し面倒くさい訳です。

そんな中で最近はreticulateパッケージを使ってRとPythonをシームレスに使う方法が楽だと感じているのでメモします。

reticulateでRとPythonを両方使う

RパッケージreticulateはRとPythonをシームレスに実行するための機能が実装されています。具体的には以下のようなことが出来ます。(※ Pythonにもrpy2という類似した機能を持つライブラリが存在していますが、後述の理由によりscRNA-seq解析においてはreticulateの方が問題がなさそうでした)

RとPythonの間で変数を共有する
R/Rstudio上でPythonコードを走らせる

library(reticulate)
use_python("/Users/maruyamatooru/opt/anaconda3/envs/single_cell/bin/python") # 使いたいpythonのパスを指定する

a = list(c(1,2), c(3,4))

np = import("numpy")
mat = np$array(a)
mat

class(mat)

     [,1] [,2]
[1,]    1    2
[2,]    3    4

[1] "matrix"

これと同様の方法でR上でScanpyなどのPythonライブラリを操ることができるようになります。Seuratで前処理・クラスタリングした結果をScanpy等のPythonライブラリ上で使ってみましょう。

library(Seurat)

# Seuratでデータを読み込む
m = Read10X("public_dataset/10x/Parent_NGSC3_DI_PBMC_filtered_feature_bc_matrix/filtered_feature_bc_matrix/")
sc_seurat = CreateSeuratObject(m, assay="RNA")

# ミトコンドリアリード率・検出遺伝子数で細胞を選出
sc_seurat[["percent.mt"]] = PercentageFeatureSet(sc_seurat, pattern = "^MT-")
sc_seurat = subset(sc_seurat, subset = nFeature_RNA > 200 & percent.mt < 10)

# 正規化
sc_seurat = NormalizeData(sc_seurat, normalization.method = "LogNormalize", scale.factor = 10000)

# PCA & UMAP & Clustering
sc_seurat = FindVariableFeatures(sc_seurat, selection.method = "vst", nfeatures = 2000)
sc_seurat = ScaleData(sc_seurat)
sc_seurat = RunPCA(sc_seurat, features = VariableFeatures(sc_seurat))
sc_seurat = FindNeighbors(sc_seurat, dims = 1:30)
sc_seurat = FindClusters(sc_seurat, resolution = 0.5)
sc_seurat = RunUMAP(sc_seurat, dims = 1:30)
sc_seurat

An object of class Seurat 
36601 features across 8295 samples within 1 assay 
Active assay: RNA (36601 features, 2000 variable features)
 2 dimensional reductions calculated: pca, umap

以上の工程で出来上がったSeuratオブジェクトをPythonのAnnDataクラスへと変換してみます。

# AnnDataへ変換
sc = import("scanpy")
anndata = import("anndata")
sc_scanpy = anndata$AnnData(t(sc_seurat[["RNA"]]@data), 
                            obs=sc_seurat@meta.data,
                            var=data.frame(Symbol=rownames(sc_seurat)))
rownames(sc_scanpy$var) = rownames(sc_seurat)
sc_scanpy$obsm = list(pca=sc_seurat@reductions$pca@cell.embeddings[,1:2],
                      umap=sc_seurat@reductions$umap@cell.embeddings[,1:2])
sc_scanpy$var["highly_variable"] = sapply(sc_scanpy$var_names$tolist(), function(x) x %in% sc_seurat[["RNA"]]@var.features)
sc_scanpy

AnnData object with n_obs × n_vars = 8295 × 36601
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.5', 'seurat_clusters'
    var: 'Symbol', 'highly_variable'
    obsm: 'pca', 'umap'

このような流れでscanpy.external API *5やその他のPythonライブラリと連携すれば一つのスクリプトで一括して解析を行うことが可能です。

Python側でRを使う場合は？

前述の通りPythonにもrpy2という類似したライブラリが存在しており、これを使うことでRをpython スクリプトに組み込むことが可能です。ただ私の環境ではSeuratがrpy2ではうまく動きませんでした。*6同様の現象は多数報告されているみたいで解決される見込みはなさそうです。*7

最近公開されたanndata2riというpythonライブラリを使うことでRのSeuratオブジェクトをAnnDataへ変換することができるようです。ただこちらも内部でrpy2が使われているようなので動くかどうかは？です。*8

*1:2018年以降シングルセルの文献が増えていることが報告されている https://www.biorxiv.org/content/10.1101/742304v1

*2:シングルセルデータ解析用ソフトウェア専門のデータベースが公開されるほど https://www.scrna-tools.org/

*3:データ前処理・クラスタリング・DEG解析・可視化など

*4:Dana Pe'er, Sarah Teichmann, Fabian Theis

*5:https://scanpy.readthedocs.io/en/stable/external/index.html

*6:rpy2でimportしようとするとクラッシュしました

*7:https://github.com/satijalab/seurat/issues/2079

*8:まだ試していないだけです

2020-09-20

Common Workflow Language (CWL)で解析の再現性を高める

CWL Bioinfo

Common Workflow Language (CWL) とは

バイオインフォマティクスでは複数のツールを組み合わせたパイプラインを構築してデータを処理することが多いと思います。例えばRNA-seq解析であれば以下のような流れが考えられます。

FASTQファイルの準備 (bcl2fastq)
シーケンスクオリティーの評価 (e.g. FASTQC)
マッピング (e.g. hisat2)
遺伝子発現量の推定 (e.g. Stringtie)

ここで問題になりがちなのが再現性です。

各ソフトウェアのバージョン、マッピングに用いるリファレンスゲノムのバージョン、GTFファイルの種類...など種々の要因を正しく管理しなければ、完全に結果を再現することは困難です。2年前にシェルスクリプトで構築したパイプラインを動かすと結果も微妙に変わった、というような経験のある人も多いのではないでしょうか。

このような問題を解決する上で役に立つのがCWLをはじめとするワークフロー言語です。

CWLは何が良いのか？

上述したRNAseqのパイプラインをCWLで構築しようとすると次のような構成になります。

f:id:auroratummy:20200920043540p:plain

用意するファイルは大きく分類すると以下の3種類です。

各ステップの要件を記述したcwlファイル (例: bcl2fastq.cwl)
パイプライン(ステップの並び)を記述したcwlファイル (例: RNAseq.cwl)
パイプラインへ与える引数 (e.g. 入力ファイル、パラメータ) を記述したyamlファイル (例: RNAseq.yml)

cwlファイルごとに異なるDocker imageを指定することが出来ます。したがって、作成したcwlファイルを別環境に持っていってもすぐに動かすことが出来ますし、依存関係の問題で動かすことが困難なソフトウェアもワークフローに楽に組み込むことが出来ます。

CWLの書き方

CWLの細かな文法に関しては公式ドキュメントに詳しく記載されています。特に

Common Workflow Language User Guideは良く出来た入門ページだと思います。（日本語訳を選択することも可能です）

しかしながらこのUser Guideには記載されていない便利なオプションが色々とありましたので備忘録も兼ねて記録しておこうと思います。以下にメモを残したワークフローはこちらに置いてあります。

github.com

CommandLineTool (各ステップのcwlファイル)

ワークフローを構成する各ステップの詳細を記載したcwlファイルを用意します。基本的には以下のような形式です。もっとオプションは色々あるようですが、よく使いそうなものを記載しました。

gistaa0012b7250665b4706e14f23f20fe02

デバッグの際には実際にcwlファイルを走らせてみると良いでしょう。cwlを走らせるには、cwlへの入力についての情報を記載したyamlファイルを用意する必要があります。

gist9813be5a7303d4c8f2fed973b9475bdb

cwlファイルとymlファイルを指定してcwltoolを起動します。この際に何らかの問題があればエラーが表示されるのでそれを元にデバッグしましょう。何も問題なければcwlファイルのoutputにて指定されているファイルが出力されるはずです。

Workflow (ワークフローのcwlファイル)

各コマンドのcwlファイルが完成した後はワークフローのcwlファイルを作成します。

gist04f8be640433260141838466ce17eb28

こちらも同様にインプットの情報を記載したyamlファイルを作成します。cwltool等で実行すればワークフローが走ります。

gistab68128bb380cf9ab115b2fe49bc4140

2020-04-05

このブログの目的

f:id:auroratummy:20200405234333j:plain

自己紹介

博士 (工学)
オミクスデータを利用した生物学研究が私の学生時代の一貫するテーマです。
大学院では発生生物学・微生物生態学・遺伝学などの研究課題に携わり、多様なデータを扱ってきました。
新しいタイプのデータを生み出すバイオテクノロジーの話が好きです。
学位取得後に国内製薬企業へ研究職として入社しました。
製薬企業でもバイオインフォマティシャンとして創薬標的の探索研究や新しい解析手法の検討・開発を行っています。楽しいです。
今年から海外赴任が始まりました。

このブログで何をしたいか？

新しい技術・研究に関する知識の整理
フォローできていない研究トピックを学ぶことが目的です。ブログを開設して、一定期間のうちに1つの記事を書く、というノルマを設定すれば自ずと目的が達成できるのではと考えました。
自分の考えを整理する
キャリア・スキルアップ・海外生活のことなど、技術・研究以外のトピックについても自分の整理する場にしようと思っています。

Aurora blog

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

シングルセル解析③：ダブレット推定手法について

ダブレットとは

Scrublet, DoubletFinder

DoubletDecon

SCDS

Co-expression based doublet scoring (cxds)

Binary classification based doublet scoring (bxds)

まとめ

シングルセル解析②：CITEseqデータの正規化手法について

CITEseqとは

CITEseqデータの正規化

CITEseq 原著論文*8の場合

Seuratの場合

DSB

結論

Local BLASTの結果からLCA (Lowest common ancestor)で由来生物種を推定する

きっかけ

準備

使い方

スクリプト

1. SQLiteでDBを作成する

accession2taxid

taxonomy

2. BLAST+

3. 由来生物種推定

シングルセル解析①: reticulateでRとPythonを両方使う

RとPythonどちらを使う？

reticulateでRとPythonを両方使う

Python側でRを使う場合は？

Common Workflow Language (CWL)で解析の再現性を高める

Common Workflow Language (CWL) とは

CWLは何が良いのか？

CWLの書き方

CommandLineTool (各ステップのcwlファイル)

Workflow (ワークフローのcwlファイル)

このブログの目的

自己紹介

このブログで何をしたいか？