2026/05/11（月）；初稿
- 1 はじめに
- 2 注意点１．RNA-seqとタンパク質（ペプチド）質量分析のデータの違い
- 3 注意点２．解析に使用できるタンパク質の数、疾患の種類、症例数がTCGAよりも少ない
- 4 注意点３．疾患ごとのタンパク質発現量の比較と生存率の解析方法はTCGAとほぼ同じ
- 5 使用するパッケージ
- 6 ここから使用するRのダウンロード１
- 7 使用するデータ
- 8 使用するデータを読み込んで整える
- 9 CPTACのタンパク質発現量は中央値で正規化されている
- 10 missForestによる欠損値の補完（imputation）
- 11 タンパク質発現プロファイルと臨床情報を結びつける
- 12 escape、gsva、fgseaを使ってGSEA（Gene Set Enrichment Analysis）を行う
- 13 疾患名の整理
- 14 タンパク質発現量やエンリッチメントスコアのヒストグラムを表示してみる
- 15 バイオリンプロットで疾患ごとのタンパク質発現量を比較する
- 16 生存解析
- 17 タンパク質発現量と遺伝子発現量の相関を解析する
- 18 よくわからないタンパク質発現プロファイルの約９０％が韓国と中国の研究機関発だった
- 19 まとめ
はじめに
これまでTCGA（The Cancer Genome Atlas）に登録された各種がんの遺伝子発現プロファイルと、その患者の生存率（生存期間）を解析してきた。この記事ではCPTAC（Clinical Proteomic Tumor Analysis Consortium）に登録されたデータを使ってがんや腫瘍のタンパク質発現プロファイルとその疾患の臨床情報を用いて、特定のタンパク質の発現量と患者の生存期間の解析をカプランマイヤー法で解析してみようと思う。
このCPTACは、種々のがんのタンパク質発現プロファイルを質量分析（Mass Spectroscopy；MS）を使って取得し、その患者に関わる臨床情報と共に登録した公開データベースである。言ってみればTCGAのタンパク質バージョンである（CPTACでプロジェクトは別物、解析自体も別物だけど）。未だ発展途上であり、登録されている疾患の数や症例数はTCGAほど多くないし、MSのデータのクオリティーも遺伝子発現ほど良くないだろうし、CPTAC側としてもデータのインテグリティーに問題があり、特にタンパク質発現プロファイルとそのファイルとリンクする臨床情報を結びつけることができない症例が有ったりする（正直これはCPTACではなく、こういうことをやってデータを無駄にするような大学や研究機関が間抜けで無責任なのだが。これは解析のときに述べる）が、このCPTACにはプロテオームだけでなく、ホスホプロテオーム、メタボローム、アセチローム、ユビキチノーム等のデータが登録されており、そんな欠点を補って余りある有用なデータベースである。昨今、そしてこれからのトランスレーショナル研究には必要不可欠なのではないかと思っている。また、同じ症例の遺伝子発現プロファイルがTCGAに登録されていたりするので、遺伝子発現とタンパク質発現の相関なんかも解析することができる。これを直ぐに使える形で解析したデータを持っておけば、とても便利だろうと思っている。
２０２６年３月になってようやくサーバーが正常に動くようになっていたので、後日、新しい版のデータでも別に解析してみようと思う。言うて、解析手法はここやこことほとんど同じはずであり、手法としてはそこまで難しくない。それに、データが少し古いからと言っても、症例数はそこまで急に増えるようなことはないので、この解析方法や解析したデータは結構長く使えるのではないかと思っている。
余談だが、生存率解析ではなく生存期間解析や単に生存解析の方が正しい呼称だと思う。言い方としては、例えば「生存解析（survival analysis）を行った結果、overall survivalの中央値は６ヶ月です。」とかになるし、overall survalの日本語約は「全生存期間」だったりするためだ。ただし、誰にでも通じる言い方として生存率解析と言ってしまうことも多いと思う。この記事では、タイトルにこそ生存率の解析としているが、本文中では生存期間解析や生存解析と言うことにする。
注意点１．RNA-seqとタンパク質（ペプチド）質量分析のデータの違い
個人的に重要と思うことが、RNA-seqのデータとMSのデータの違いである。これまではRNA-seqの技術が圧倒的に早くに発展、普及してきたので、この手の解析は遺伝子発現プロファイルに基づくことが多かった。そして、その勢いでMSのデータを解析すると、おそらくどこかで躓くような気がする。RNA-seqの解析から入った人間は、MSのデータの性質がRNA-seqと大きく異なる事を忘れてはならない。RNA-seqの場合、昨今ではハイスループットシークエンサーが利用されている。これは読んできた（リードした）シークエンスを、ゲノム（というかエクソン上）にマップして、ゲノム上にどのくらいマップされたのかカウントしていく。データは「カウント値」として記録されるので、二項分布（負の二項分布というヤツ）に従う。それに、ゲノム上にどのくらいマップされなければならないかをカバレッジなどの指標にしてシークエンスを計画するし、ゲノム上に似たような配列があれば、そこにマップされてしまったりと、けっこうな領域をカバーすることができる。それにデータがカウント値であることから、欠損値（欠損値というのかはわからん）が出た場合は、それをカウント値がゼロとして解析を行うことができる。一方、MSのデータは、質量分析の結果でてきたフラグメントのシグナル強度を計測し、それをデータベースに照合して解析する。これはカウント値ではない（非常に申し訳ないが、自分はこのくらいしか知識がない。chatGPTの方が詳しく教えてくれる）ので欠損しているタンパク質をカウント値がゼロと見なすことが難しいし、検出されてこない原因が単に発現が低いのか、技術的に拾ってくるのが難しかったためなのか、本当に発現していないのかよくわからない。また、欠損値をカウント値がゼロとみなして解析するにしても、欠損値が多すぎるような気もする。統計分布と欠損値の取り扱いが異なるので、RNA-seqの解析でよく利用されるDESeq2やedgeRはDEG（Differential Expressed Gene）の解析には不向き（使えないこともないが。）であり、他の適切なパッケージを利用するほうが良い。また、利用するパッケージにも依るが、解析するためには欠損値は他の値（欠損値ではない値）を参考にして補完（imputation）しなければならない。それに欠損値を補完するか自体よく考えた方が良い場合もあるように思う。言うて解析する側は、その目的、例えば、何らかの群同士を比較してDEGを出すのか、この記事での解析のように単に発現で分類して他の解析との関連をみるのかによって、適切な方法で解析を進めなくてはならない。
因みに、limmaとかなら使えると思う。limmqはもともとcDNAマイクロアレイのDEG用だったはずであり、これはcDNAマイクロアレイは蛍光シグナルを検出しているためだ。
後述するが、この記事ではDEGの検定などは行っていないので、そのためのそのようなパッケージは特に使用していない。しかし、欠損値の推定はゼロを使うよりは良いということで、missForestというRのパッケージを使って行っている。ここがRNA-seqを元にした解析と異なる点である。
注意点２．解析に使用できるタンパク質の数、疾患の種類、症例数がTCGAよりも少ない
はじめにで書いたように、CPTACは未だ展途上である。バージョンだってまだ５くらいであり、疾患の数がTCGAと比較しても少ない。特に痛いと思うのがこの点で、この解析でも最終的に１１疾患くらいしかなかった。依って、解析に必要な疾患のデータが十分でなかったり、そもそも関心のある疾患が無い可能性もある。その場合は他のCPTACだけでなく他のプロジェクトも使用する必要が出てくるが、プロジェクトを横断する場合は登録するデータの違いなども問題として出てくるに違いない（知らんけど）。それに加えて、タンパク質発現プロファイルのデータと、患者の情報をうまく結合できない症例もある。このあたりを考えて、CPTACやPDC（Proteomics Data Portal）に登録のあるデータを使用する必要がある。
注意点３．疾患ごとのタンパク質発現量の比較と生存率の解析方法はTCGAとほぼ同じ
タンパク質発現量と患者データの確認や結合の方法は、詳細は違えど結局似たようなことをやるし、疾患ごとのタンパク質発現量の比較や生存率の解析方法は以前にTCGAのデータを解析したときとほぼ同じである。
あと、この記事データのクリーニングとか整形ばっかりで本当に面白くない。９割がデータのクリーニングである。もし買う人がいたら、本当に注意したほうが良い。本当に面白くない。むしろ苦痛である。
早速以下からデータの解析を始める。
使用するパッケージ
以下がこの解析で使用するパッケージである。上記で「RNA-seqとデータ取得の方法や考え方が違うから注意」とかなんとか言っておきながら、いつもRNA-seqの解析で使用している基本パッケージを使用した。以下のコードにDEGとか初っ端に書いてあるけど、ここではDEGは行わないのでこれらは使用しない。一方、エンリッチメント解析のパッケージであるfgsea、escape、GSVA、並列化に関するパッケージforeach、doParallel、欠損値imputationにmissForest、生存率解析用にsurvminer、データの読み込みにdata.table、tidyverse、gridExtra等を使用する。以下を読み込んでおけばこの記事の解析は一通りできるはずである。ここで使用するエンリッチメント解析のパッケージであるescape、gsva、fgsea、そしてここでは解析に費やす時間の関係上実際には使用しないが（動くコードは書いた）、ssGSEA2は、本来遺伝子発現解析のために使用されると思うが、エンリッチメント解析としてやっていることはタンパク質発現解析にも利用出来るはずである。実際、ssGSEA2はPTM（Post Transcriptional Modification）の解析にも使用されているようである（やったことはないが）。
# DEG
library(baySeq)
library(edgeR)
library(DESeq2)
library(limma)
# enrichment analysis
library(clusterProfiler)
library(fgsea)
library(ssGSEA2) # devtools::install_github("nicolerg/ssGSEA2") # devtools::install_github(broadinstitute/ssGSEA2.0) and devtools::install_github(broadinstitute/ssGSEA2) did not work.
library(escape)
library(GSVA)
library(pathview)
# immune cell deconvolution
library(ConsensusTME)
library(estimate)
library(immunedeconv)
library(MCPcounter)
library(msigdbr)
library(quantiseqr)
library(TIMER) #devtools::install_github('hanfeisun/TIMER'); it will not work!! use TIMER of immunodeconv if you need.
library(xCell)
# stats
library(NSM3) # sudo apt install libgmp-dev libgmp10 libgmp3-dev
library(pwr)
library(missForest)
# basic
library(beeswarm)
library(circlize)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(doParallel)
library(foreach)
library(parallel) # Need it for use of baySeq
library(pheatmap)
library(rtracklayer)
library(rvest) # this is for extraction of html table; results of leading edge analysis.
library(xml2) # this is for extraction of html table; results of leading edge analysis.
library(openxlsx) # for read xlsx file.
library(ComplexHeatmap)
library(tidyverse)
library(data.table)
library(survminer)
library(survival)
library(gridExtra)
library(MASS)
library(coin)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(RhpcBLASctl)
library(psych) # for geometric mean
# # shiny
# library(shiny)
# library(bslib)
# Check global environment
Sys.getenv()