基礎医学生物学研究者のブログ

TCGAのデータとRを使ってがんや腫瘍の遺伝子発現量と患者の生存率の関係を解析する（TCGA, v.43, 2025）

R

2026/01/03（土）；初稿
2026/01/12（月）；追記；AnndataではなくSummerizeExperimentやMultipleAssayExperimentが使えると思う。

目次

はじめに

以前、TCGAのデータを使って、がんや腫瘍における遺伝子発現量（RNA-seqのカウント値）と患者の生存率の関係をRで解析したが、そのときは対応する正常組織（おそらくがんや腫瘍周辺の正常部位）におけるRNA-seqのカウン値をダウンロードしなかったため、それらをバイオリンプロットを用いた疾患毎の遺伝子発現量もしくは遺伝子セットのエンリッチメントスコアの比較に含めることが出来なかった。それに、TCGAのデータも２０２４年のバージョン４１であり、ちょっと時間が経っていたので、改めてダウンロードし、解析することにした。

前置きはここの記事と一緒になるので、ここでは述べないことにする。興味があれば、ぜひ読んでいただきたいと思う。また、ここに以前の解析との比較表を載せた。この表を見て必要な解析が載っている方を選んでいただきたい。

早速本題に入ろうと思う。

注意点１

以前の解析との違いは、以下である。

  1. 正常組織の遺伝子発現量を解析、特にバイオリンプロットに含めた。
  2. GDC portalからTCGA（とCPTAC）で採られたRNA-seqのカウント値をダウンロードした。
  3. ssGSEA2とfgseaは計算が遅すぎたので使用しなかった。escapeとGSVAのみでGene Set Enrichment Analysis（GSEA）を行った。
  4. MSigDBに登録されている（２０２５年１０月時点）ヒトの全遺伝子セットを使ってGSEAした。

ここでやっている解析は以前の解析とほぼ同じであり、もしかしたら以前の解析の方が良いかもしれない。特に、後ほどバイオリンプロットで各疾患毎さらに対応する正常組織の遺伝子発現量もしくは遺伝子セットのエンリッチメントスコアを比較するが、正直、結論を得るためには正常組織の症例数が足りないように思うし、当然ながら正常組織を採取できない症例や疾患もあるためだ。正常組織も解析に含めてみて、これに気がついた。今なら、「正常組織も解析に含めることは出来ますが、症例数がちょっと足りないようにも思いますし、そもそも正常組織を採っていない疾患も多いです。もしかしたら、正常組織との比較は参考程度にしかならず、それで何かを決定するのは微妙かもしれません。」と自信を持って言える。

なお、CPTACについては今後、CPTACで採られたタンパク質発現量と遺伝子発現量の相関を解析するためである。ここではデータのクリーニングのときにちょっとだけ出てくるが、この記事ではそれらを解析しない。以下に比較表を載せておく。

項目この記事以前の記事
タイトルTCGAのデータとRを使ってがんや腫瘍の遺伝子発現量と患者の生存率の関係を解析する（TCGA, v.43, 2025）TCGAのデータとRを使ってがんや腫瘍の遺伝子発現量と患者の生存率の関係を解析する（TCGA, v.41, 2024）
TCGAのバージョンv.43, 2025v.41, 2024
GDC portalからのデータのダウンロード方法ありなし
GMTファイルの作り方なしあり（ただし、Geneontologyから選んできた遺伝子セット１つだけ）
使ったGMTファイルMsigDBに登録されているヒトの遺伝子セット全部（Human MSigDB v2025.1.Hs、２０２５年１０月時点）Geneontologyから選んできて自分で作成したHRR（Homologous Recombination Repair）と、MsigDBのH: hallmark gene sets（Human MSigDB v2025.1.Hs、２０２５年１０月時点）のみ。
エンリッチメント解析GSVA
escape		ssGSEA2
escape
GVSA
fgsea
バイオリンプロットによる疾患毎の遺伝子発現量もしくは遺伝子セットのエンリッチメントの比較正常組織あり
ウィルコクソンの順位和検定あり（ただし、正常組織の例数が腫瘍組織に比べて少ないのでちょっと注意が必要と思う。）		常組織なし
ウィルコクソンの順位和検定なし
生存期間解析でのログランク検定ありあり
その他比較するデータの数（リストの要素の数）が多く、そのために少しデータの準備がすこし複雑でメモリを喰いまくる。比較するデータの数（リストの要素の数）がそんなに多くないので、データの準備がそんなに大変じゃない。

注意点２

以前の解析でも書いたが、この解析で使用するコンピューターの性能はずいぶんと高いと思う。以下の画像の通りメモリを合計で1.3TB以上も使っており、たとえRandom Access Memory（RAM）（物理メモリ）を256GB積んでいたとしても全く足りていないことがわかる。ちなみに、これはデータセットの全部がすべて整って、生存率曲線の解析とそれに続くログランク検定を行っているとき（３つめのRのコード）のメモリの使用量である。１つ目から３つ目のコードを一つのRStudioに入れようとしても、作業領域が一杯になってまともに動かなくなってしまう。また、昨今のノートPCや軽量のラップトップでは計算能力が足りない可能性があるので、要注意である。もはやサーバーのようなコンピューターが必要と言える。

GDC portalからTCGAのカウント値をダウンロードする

まずは、GDC portal（https://portal.gdc.cancer.gov/）からTCGAとCPTACのRNA-seqのカウント値をダウンロードする。ダウンロード方法については以前の記事が詳しいので、そちらを参照したほうが良いかもしれない。GDC portalは以下のようなフロントエンド（２０２５年１２月時点）である。まず「Cohort Builder」を押す。

「Cohort Builder」を押すと、以下のようなページに入る。「Program」というパネルには様々なプロジェクト（の大元）があるので、TCGAとCPTACにチェックを入れて、そのコホートに名前を付けて保存する。こうすれば一々TCGAとCPTACを選び直して…ということが無くなる。ここでは2025_10_18_DLという名前で保存した。次に、「Repository」をクリックする。

「Program」のところにあるTGCAとCPTACにチェックを入れる。
Cohort Builder値
ProgramTCGA
CPTAC
Cohort Builderでチェックを入れるところ。

次に「Repository」をクリックし、Filterのパネルのところを、以下の表のようにチェックを入れていく。なお、ダウンロードしたい症例はCartに入れてからまとめてmanifest fileなどをダウンロードするが、このとき、Cartには10,000ファイルしか入らなかったと思うので、以下のように順番１から３のように分割して、Cartに入れるファイル数を10,000より少なくしてからダウンロードすることにした。

ダウンロード１回目

まずは１回目のダウンロードが以下の表である。

Filter値順番
Accessopen1
Experimental StrategyRNA-seq1
Tissue typetumor
1
Tissue Descriptormetastatic
primary
recurrence		1
Preserved Methodffpe
frozen
oct
1
１回目の値のまとめ。

各パネルの画像は以下のような感じである。

Access、Tissue type、Tissue Descriptorのパネル。
Experimental Strategyのパネル。
Preservation Methodのパネル。
１回目のダウンロードでのファイル数、症例数、総容量。

上記を選んだら、「Add All Files to Cart」をクリックし、次に右上にある「Cart」をクリックする。

選んだファイルをカートに入れる。

カートに移動したら、Manifestファイル、サンプルに関するデータ（Biospecimen）、臨床情報に関するデータ（Clinical）、Sample Sheet、Metadataをダウンロードする。Sample SheetとMetadataは表示の「Sample Sheet」、「Metadata」のところをクリックするだけである。他は以下のようにダウンロードする。BiospecimenとClinicalに関するファイルは圧縮されているので、適宜展開する。

「Download Cart」の「Manifest」をクリックすると、Manifestファイルがダウンロードされる。
「Biospeciemen」の「TSV」をクリックすると、サンプルの情報に関するデータ一式がダウンロードされる。
圧縮されて入っているサンプルに関する情報。
「Clinical」の「TSV」をクリックすると、症例毎の臨床情報に関するデータ一式がダウンロードされる。
圧縮されて入っている臨床情報に関する情報。

ダウンロード２回目

以下がダウンロード２回目の値である。対象となるファイル選んでカートにいれたら、カートに移動して手順１の通り各種ファイルをダウンロードする。

Filter値順番
Accessopen2
Experimental StrategyRNA-seq2
Tissue typetumor
2
Tissue Descriptormetastatic
primary
recurrence		2
Preserved Methodsnap frozen
unkown		2
２回目の値のまとめ。
２回目のダウンロードでのファイル数、症例数、総容量。

ダウンロード３回目

以下がダウンロード３回目の値である。３回目は正常組織のみを選んだ。これは、後ほどバイオリンプロットにより疾患間の遺伝子発現量もしくは遺伝子セットのエンリッチメントスコアを比較するときに、正常組織の発現量も加えるためである。手順１と２と同様に、対象となるファイル選んでカートにいれたら、カートに移動してダウンロード１回目の通り各種ファイルをダウンロードする。

３回目のダウンロードでのフィルターの値。
Filter値順番
Accessopen3
Experimental StrategyRNA-seq3
Tissue typenormal3
３回目（正常組織のみ）の値のまとめ。
３回目のダウンロード（正常組織のみ）でのファイル数、症例数、総容量。

gdc-clientを使ってManifest fileに記載されたファイルをダウンロードする

以下のコードがManifest fileを使ったRNA-seqのカウント値のダウンロード方法である。gdc-clientをここからダウンロードし、インストールする必要がある。これらの方法は以前の記事にある程度詳しく書いてある。以下がgdc-clientを使ったダウンロードの方法である。

Bash
# 1
gdc-client download -m "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/1_ffp_frozen_oct/gdc_manifest.2025-10-25.154055.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_tumor" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/1_ffp_frozen_oct/gdc_manifest_1_log.txt"

# 2
gdc-client download -m "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/2_snapfrozen_unkown/gdc_manifest.2025-10-25.160033.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_tumor" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/2_snapfrozen_unkown/gdc_manifest_2_log.txt"

# 3
gdc-client download -m "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/3_normal/gdc_manifest.2025-10-25.160847.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_normal" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/3_normal/gdc_manifest_3_log.txt"

# 3 retry
gdc-client download -m "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/3_normal/gdc_manifest.2025-10-26.024514.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_normal" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/3_normal/gdc_manifest_3_retry_log.txt"

# 1 retry
gdc-client download -m "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/1_ffp_frozen_oct/gdc_manifest.2025-10-26.103837.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_tumor" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/1_ffp_frozen_oct/gdc_manifest_1_retry_log.txt"

# 2 retry

/mnt/team4tb/Dropbox/Downloads Safari/gdc_manifest.2025-10-27.003414.txt
gdc-client download -m "/mnt/team4tb/Dropbox/Downloads Safari/gdc_manifest.2025-10-27.003414.txt" -d "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/count_tumor" --log-file "/mnt/seqdata/public_data/Blog/2025 10 18 TCGA/2_snapfrozen_unkown/gdc_manifest_2_retry_log.txt"
Bash

使用するパッケージ

以下のパッケージが自分がいつもRNA-seqの解析に使っているものである。TCGAのデータを使った生存率解析のためのパッケージとして必須なものは、tidyverse、data.table、suvival、survminer、 grid.Extraである。それに加えて、この解析ではescapeとgsvaを用いる。

R
# 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)

# # shiny
# library(shiny)
# library(bslib)

# Check global environment
Sys.getenv()
R