初めに
この記事について
この記事はRのパッケージBioconductorに収録されているTDbasedUFEを最低限利用できるために書かれた記事です。パッケージ内で使われている数理的原理の解説は殆ど行いません。TDbasedUFEの詳細はBioconductorの該当ページや田口善弘教授の論文を参照してください。
筆者
中央大学理工学部物理学科B4で、研究室は異なりますが田口善弘教授から直に研究の指導を受けています(記事作成当時：2023/9)。本所属の研究室は素粒子論です。またR歴は半年なのでプログラムが雑かと思います。
TDbasedUFEについて
端的に言えば、RNA-seqやATAC-seqなど遺伝子解析からテンソル分解という手法を用いて特徴ある遺伝子を教師なし学習で抽出できるパッケージです。中央大学の田口善弘教授によりBioconductorに誰でも使いやすい形で収録されます。更に機能を拡張したTDbasedUFEadvというパッケージもありますが、今回は無印のTDbasedUFEのみの解説です。また工夫次第でマルチオミックスの解析についても可能です。
今回はチュートリアルということでRNA-seqの解析を用いて、2つの条件下で異なる発現をする遺伝子を探します。この記事ではデータの取得から遺伝子抽出まで行いますので、バイオインフォマティクスに精通している人はデータ取得などは飛ばしても構いません。私自身、バイオインフォマティクスについては修業中の身ですので、軽い気持ちで雰囲気だけでも知っていただければ幸いです。
データの取得から前処理
データをダウンロード
National Center for Biotechnology Information (NCBI)のGene Expression Omnibus (GEO)からデータを取得します。NCBIは実験データだけでなくタンパク質やSNPなども取りまとめている、生物学に関わるデータベース運営・ソフトウェア開発を行っている機関です(Wiki)。
今回はGEO上のデータセットであるGSE176154を解析します。選んだ理由はGEO DataSetsに"hiPS RNA-seq"と適当に検索したら上位にあった、それだけです(解析するのはES細胞ですが)。
GSE176154はヒトES由来心筋細胞(hiPS-CMs)のRNA-seqデータで、37日間増殖させた細胞と35日目にもう一度カルチベートさせた細胞の対比実験です。さて、ページの一番下にある"GSE176154_Galaxy26-..."のhttpをクリックしてデータをダウンロードします(画像)。
解凍して、タブ区切りのRNA-seqの解析結果がテキストファイルに書かれていることが分かると思います。今回の解析では、37日齢の細胞と35日目に再培養した細胞で異なる発現をする遺伝子を探します。
次にこれをRで読み込もうとしたのですが、Rにデータに欠損があるとか何とか言われて中途半端にしか読み込めませんでした。対処策としてExcelを開き、[データ]>[データの取得]>[ファイルから]>[テキストまたはCSVから]で先ほどダウンロードしたテキストファイルを開いたところ全て読み込みました。Excel君は優秀ですね。ファイル名を"hES.csv"として保存しました。
Rでデータの確認と整形
Rで先ほど作ったcsvファイルを読み込みましょう。ちなみに私はR Studioを用いて、R Markdown形式でGoogle Colaboratoryのようにチャンクごとに実行してるので、ソースコードは小出しになりますのでご了承ください。
hES <- read.csv("hES.csv", header = TRUE) # csvファイルの読み込み
head(hES) # データフレームの上部を表示
names(hES) # 列名をすべて表示
【出力結果】
| sort0 | Gene.ID | Gene.Name | ... | X1..N2148.txt | ... | |
|---|---|---|---|---|---|---|
| 1 | 1 | ENSG00000175592.9 | FOSL1 | 302.4622 | ||
| 2 | 2 | ENSG00000157927.17 | RADIL | 100.8207 | ||
| 3 | 3 | ENSG00000106366.9 | SERPINE1 | 7450.6530 | ||
| 4 | 4 | ENSG00000118523.6 | CCN2 | 46779.8200 | ||
| 5 | 5 | ENSG00000119138.4 | KLF9 | 1032.4040 | ||
| 6 | 6 | ENSG00000105894.12 | PTN | 1278.4070 |
[1] "Sort0" "Gene.ID" "Gene.name" "Gene.description"
[5] "Chr" "Start" "End" "Strand"
[9] "Length" "IGV.coordinates" "X1..N2148.txt" "X2..N2150.txt"
[13] "X3..N2152.txt" "X4..N2147.txt" "X5..N2149.txt" "X6..N2151.txt"
[17] "Base.mean" "log2.FC." "StdErr" "Wald.Stats"
[21] "P.value" "P.adj"
見覚えのあるDESeq2形式ですね。正規化された遺伝子の発現量が"X1..N2148.txt"の列にあることが分かりました。またnames()でそれぞれの列名を表示させ、11列目の"X1..N2148.txt"から16列目の"X6..N2151.txt"に遺伝子発現量があることが分かりました。
TDbasedUFEではこの遺伝子発現量と遺伝子名を使います。遺伝子名は解析に影響は無いので、SymbolやEnsemble.IDでも何でも良いです。自分の分かりやすい文字列にしましょう。
さて、用があるのは遺伝子発現量と遺伝子名だけですので、そこだけ抽出しちゃいましょう。
hES <- hES[ , c(3, 11:16) ] # 1,11~16行目を抽出
names(hES) <- c("Gene", "N2148", "N2150",
"N2152", "N2147", "N2149", "N2151") # 列名を変更
head(hES) # データフレームを表示
【出力結果】
| Gene | N2148 | N2150 | N2152 | N2147 | N2149 | |
|---|---|---|---|---|---|---|
| 1 | FOSL1 | 302.4622 | 125.4727 | 270.79910 | 17.39651 | 16.60064 |
| 2 | RADIL | 100.8207 | 104.0505 | 50.91023 | 584.15630 | 430.69450 |
| 3 | SERPINE1 | 7450.6530 | 5687.0760 | 13836.75000 | 975.11990 | 1449.79000 |
| 4 | CCN2 | 46779.8200 | 32702.4700 | 44763.09000 | 11950.48000 | 14800.40000 |
| 5 | KLF9 | 1032.4040 | 989.5003 | 754.98780 | 149.24370 | 266.53260 |
| 6 | PTN | 1278.4070 | 1649.5070 | 1287.92000 | 3341.04500 | 3449.24500 |
これで使いやすいデータフレームになりました。1列目を遺伝子名に、2列目以降を遺伝子発現量としています。またGSE176154から被検体の情報を読み取ると、"N偶数"が37日齢のES細胞、"N奇数"が35日齢に再培養したES細胞とのことなので、次は検体を条件ごとに整列させていきます。
テンソルの作成
TDbasedUFEの解析に必要なのは、Feature(遺伝子名や遺伝子座位など)、Sample(サンプル名)、そしてValue(遺伝発現量やスコアなど)の3つです。今回の場合は、Featureは表の1列目、Sampleは列名、Valueは表にある遺伝子発現量となります。また遺伝子発現量は、遺伝子とサンプルと条件の3次元データであるともいえます。つまり、$k$番目の条件の下で$j$番目のサンプルの$i$番目の遺伝子の発現量が$x_{ijk}\in \mathbb{R}^{N \times 3 \times 2}$だと見なします。これはまさに3階テンソルであり、3次元の直方体の配列を作成していきます。
今回のテンソルの概形は図の通りです。
遺伝子発現量$x_{ijk}$が3階テンソル量ですので直方体に遺伝子発現量が並びます。イメージとしては、現在2次元の表になっているデータフレームをぱたんと折って、少し厚みのある直方体にする感じです。
少し話が難しくなったところですが、これはValueを作る際の話ですので、まずは簡単なFeatureとSampleを作成しましょう。Featureは単に遺伝子名のベクトルを用意すればよく、Sampleはサンプル名を$3\times 2$に整形すればよいので、
Feature <- hES[,"Gene"] # 遺伝子名をベクトルとして抽出
Sample <- names(hES)[2:length(names(hES))] # 2列目以降のサンプル名を抽出
Sample <- matrix(Sample, ncol = 2, byrow = FALSE) # 3 x 2に整形
head(Feature) # 1.遺伝子名を表示
length(Feature) # 2.遺伝子の総数[N]を表示
Sample # 3.サンプル名のテンソルを表示
【出力結果】
[1] "FOSL1" "RADIL" "SERPINE1" "CCN2" "KLF9" "PTN"
[1] 60669
[,1] [,2]
[1,] "N2148" "N2147"
[2,] "N2150" "N2149"
[3,] "N2152" "N2151"
上の表と合わせて見ると上手く行ってそうですね。遺伝子数は$N=60669$であることが分かり、Sampleは$\mathbb{R}^{3\times 2}$の2階テンソルになってます。$k=1$列目が37日齢で$k=2$列目が35日齢です。またこの2つはValueの添え字(番地)でもあります。一般化して、データのValueが$n$階テンソルなら、Sampleは$(n-1)$階テンソルになるわけですね。
ここで、Sampleの添え字の順番を覚えてください。後にテンソル分解をする際にこの順番が大切になります。今回の場合は、一番目(左)の添え字は細胞を表し、二番目(右)の添え字は条件です。添え字が細胞⇒条件の順番であることを覚えておきましょう。大きさが3と2であることからも分かると思います。
最後にSampleに合わせてValueを整形すればよいので、雑にfor文を用いて、
Value <- array(dim = c(length(Feature), dim(Sample) ) ) # 空のテンソル作成
dim(Value) # Valueテンソルの次元を表示
# Sampleの通りに遺伝子発現量をValueに格納
for (cond in 1:dim(Value)[3]) {
for (samp in 1:length(Sample[ , cond]) ) {
Value[ , samp, cond] <- hES[ , Sample[samp, cond]]
}
}
# 2番目の条件の1番目の細胞(N2150)のValueを表示
Sample[2, 1] # サンプル名
head(Value[ , 2, 1]) # Valueの始め
【出力結果】
[1] 60669 3 2
[1] "N2150"
[1] 125.4727 104.0505 5687.0760 32702.4700 989.5003 1649.5070
上手くいってますね。Rに慣れて自身がある方は確認は最低限で十分かもしれませんが、心配性なのでチャンクごとに確認してしまいます。
以上でテンソルの作成は完了です。後はTDbasedUFEの関数を用いて解析する流れ作業ですので、山場は越えました。
TDbasedUFEの実行
さて本記事のメインです。上で作成したテンソルを分解して遺伝子を抽出してみましょう！
TDbasedUFEのインストール
Bioconductor経由で簡単にインストールすることができます。以下のコードを実行してください。環境につき1度だけで十分です。
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("TDbasedUFE")
もしここで、一部のライブラリーがインストールできずにエラーが出た場合は、該当するライブラリーをGitHub経由でインストールしてください。稀にBioconductorからのクローニングが上手く行かない場合があります(私はありました)。
次にライブラリーをインポートして下さい。これはRを起動した毎に必要です。
library("TDbasedUFE")
テンソルの下準備
はじめに、Feature, Sample, Valueを1つの変数にまとめる必要があります。特に重要な計算をしている訳でもなく、テンソル分解する前の枕詞だと思って下さい。
TENSOR <- PrepareSummarizedExperimentTensor(feature = Feature,
sample = Sample,
value = Value)
テンソル分解(HOSVD)
次に、いよいよテンソル分解を行います。HOSVDというアルゴリズムを用いてテンソルを1つのコアテンソルと次元の数だけある行列に分解します。イメージとしては、高次元の対角化だという理解で遺伝子を抽出するだけなら問題ありません。
HOSVD <- computeHosvd(TENSOR)
# HOSVD <- computeHosvd(TENSOR, dims = c(6, dim(attr(TENSOR, "value"))[-1]))
computeHosvd()でテンソル分解を実行します。上のコードで動く方は気にしなくてよいですが、環境によってはメモリ不足でエラーを吐かれます。HOSVDのテンソル分解は一意ではなく、ざっくり言うとコアテンソルのサイズが計算量に関わってきます。そのコアテンソルの(Feature方向の)大きさがコメントアウトした下の文のdims=であり、デフォルトでは10となっています。10では私の環境ではメモリ不足だったので、6にしました。エラーを吐かれた際は、皆さんのPC環境に合わせてコアテンソルの大きさを調整して下さい。
特異値ベクトルの選択
今回の解析の目的をおさらいします。2つの条件下で異なる発現をする遺伝子を見つけたいので、差が出て欲しいのは発現量$x_{ijk}$の条件の添え字である$k$です。細胞の添え字である$j$では差が出ない方が望ましいです。
今から、この条件に合致する遺伝子群を選択してきます。この遺伝子群と対応する数値ベクトルを特異値ベクトルと呼びます。というのも、HOSVDではコアテンソルの要素の数だけ特異値ベクトルの候補があります。教師無しで分解したので、その中には当然全く興味のない遺伝子群もあるという訳です。
では特異値ベクトルを選んでいきましょう。TDbasedUFEの関数selectSingularValueVectorSmall()を用います。
# 興味のあるベクトルの選択
input_all <- selectSingularValueVectorSmall(HOSVD)
すると、この様なウィンドウが出てきます。
[Next][Prev]を押して興味のあるベクトルを選びます。ここで添え字の順番が細胞⇒条件であったことを思い出してください。始めに表示されているこのベクトルは細胞を表しています。つまり、差が出ない方が興味があるので、[Select]を押して、この画像のベクトルを選択します。
すると、また同じレイアウトのウィンドウが立ち上がります。続いては条件のベクトルです。こちらは差が出て欲しいので、1番目と2番目で正負が反対であるものを選びましょう。
今回、私が選んだベクトルは以下の画像です。左の1が細胞で、右の2が条件です。とても上手く分解できた特異値ベクトルが選ばれたと思います。ちなみにこのベクトルは規格化(大きさが1)されてますので、グラフ間での数値の大きさは気にしないで大丈夫です。
ベクトルを選び終わると、特異値ベクトルが自動で選ばれます。そしてその特異値ベクトルでテンソル分解された結果が自動でプロットされます。今回は以下のような2つのグラフが結果としてプロットされました。
左のグラフは最適化された標準偏差を表すのですが、今回はチュートリアルということで説明は割愛します。右側のヒストグラムは、補正$P$値$P_i$に対して$(1-P_i)$ごとに抽出された遺伝子の数を表します(両側検定)。つまり右側の方が$P$値が低く、有意に差がある遺伝子という訳です。上の画像のような二峰性のヒストグラムが出れば成功です。左側の大きな山は興味のないありきたりな遺伝子で、右側の山が目的の遺伝子になります。ちゃんと遺伝子を分けて抽出されたことが確認できました！
遺伝子の抽出
最後は流れ作業です。以下のソースコードを実行してください。
# p値を決めて遺伝子を抽出
index <- selectFeature(HOSVD, input_all = input_all, p0 = 0.01)
# 抽出した遺伝子を表にまとめる
DEG <- tableFeatures(TENSOR, index)
head(DEG) # 結果の表示
【出力結果】
Feature p value adjusted p value
3 SERPINE1 0 0
4 CCN2 0 0
6 PTN 0 0
7 THBS1 0 0
8 VSIR 0 0
9 PCDHAC2 0 0
DGEに抽出された遺伝子が$P$値とともに記録されます。selectFeature()の引数p0=で、ご自身の目的に合わせた有意水準$P_0$を設定して下さい。今回は$P_0=0.01$で、これはデフォルト値でもあります。出力結果を見ての通り、$P$値が0.01以下の遺伝子が抽出されてます！0が並んでいるのは、コンピューターが丸めてしまうほど$P$値が小さいからです。
TDbasedUFEまとめ
以上を持ちまして、TDbasedUFEでの遺伝子の抽出が終わりました！実験データのテンソルさえ作ってしまえば、後は直感的で簡単に遺伝子が抽出できたことが分かったと思います。
TDbasedUFEのメリットとデメリットは次の通りだと思っています。
-
メリット
- テンソルさえ作るこができれば簡単に遺伝子が抽出できる
- (3次元以上の)多次元の情報を処理でき、マルチオミックスも可能
- 教師無しなので、分類ができていない集団でも使える
- 計算量が比較的軽い。一昔前のノートPCでも今回の解析は十分可能
-
デメリット
- 結果の解釈が難しい。例えば今回抽出できた遺伝子は 37日齢と35日齢で特異的に発言した遺伝子ではなく、線形代数学でいうところの写された別空間で有意に分解された遺伝子である。
- ライブラリの関数は最低限しか無いので、片側検定や条件で特異値ベクトルの自動抽出などには工夫が必要
- 量的でないデータを直接処理できない
テンソル分解は画像処理や深層学習などで名を馳せていますが、バイオインフォマティクスでも威力を発揮します。残念ながら執筆当時はあまり流行っていない手法でありますが、この記事を期にTDbasedUFEに触れて頂ければ幸いです。
おまけ：Enrichment解析
遺伝子を抽出しただけではあまり意味が無いので、最後にEnrichment解析を行います。今回はRを用いているので、同じBioconductorに収録されているclusterProfilerを用います。項目が多く動作も軽いEnrichrで解析してもよかったのですが、可視化が容易でオプションも多いclusterProfilerにしました。今回はGene Ontology解析(GO解析)のBiological Processで解析を行います。
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
# Symbolをentrez.idへ変換
Gene_Table <- select(x = org.Hs.eg.db, keys = DEG$Feature,
keytype = "SYMBOL", columns = "ENTREZID")
# Enrichment解析(有意水準=0.01)
GOBP <- enrichGO(gene = Gene_Table$ENTREZID, OrgDb = org.Hs.eg.db,
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01,
qvalueCutoff = 0.01, readable = TRUE)
GOBP2 <- setReadable(GOBP, 'org.Hs.eg.db', 'ENTREZID')
GOBP2 <- pairwise_termsim(GOBP2) # 解析結果の類似性を計算
BAR <- barplot(GOBP, showCategory=20) # 棒グラフの描画
EMAP <- emapplot(GOBP2, pie="count", cex_category = 1.5) # ネットワークを描画
library(cowplot)
library(ggplot2)
plot_grid(BAR, EMAP, ncol=2, labels=LETTERS[1:2]) #並べてプロット
【出力結果】
核と筋線維(アクチン、ミオシン)についてクラスターが見つかりました。培養を再び行うとこれらが活発に影響があるのでしょうか？
参考文献
- Taguchi Y h., Turki T. Application note: TDbasedUFE and TDbasedUFEadv: bioconductor packages to perform tensor decomposition based unsupervised feature extraction. Front. Artif. Intell., 01 September 2023
Sec. Medicine and Public Health
Volume 6 - 2023 | https://doi.org/10.3389/frai.2023.1237542. - Taguchi Y h. TDbasedUFE: Tensor Decomposition Based Unsupervised Feature Extraction. Bioconductor version: Release (3.17); 2023. https://doi.org/doi:10.18129/B9.bioc.TDbasedUFE.
- Wu T, Hu E, Xu S, et al. Clusterprofiler 4.0: a universal enrichment tool for interpreting omics data. The Innovation. 2021;2(3):100141. https://doi.org/10.1016/j.xinn.2021.100141