オミクス解析入門 ― トランスクリプトーム解析の実践例 ―

1. はじめに

生命現象を分子レベルで理解するためのアプローチとして、オミクス解析が広く用いられています。オミクス解析とは、生体内の分子を網羅的に計測・解析する手法の総称であり、RNA、タンパク質、代謝物などの対象分子に応じた各種解析（トランスクリプトーム解析、プロテオミクス解析、メタボロミクス解析など）に分類されます。いずれも個別の分子ではなく、生体システム全体を分析することで、疾患メカニズムの解明やバイオマーカー探索に貢献しています。

中でも、発現している RNA を網羅的に定量するトランスクリプトーム解析は、次世代シーケンサー (NGS) の発展とともに広く利用されているオミクス解析の一つです。RNA-seq は mRNA 由来の配列情報を大規模に取得することで、数万種類の遺伝子の発現量を同時に定量できる強力な手法です。

RNA-seq で得られたデータからは、目的に応じて様々な情報を引き出すことができます。代表的な解析として、以下のようなものがあります。

• 発現変動解析 (DEG 解析) — 疾患群と健常群など、条件間で発現量に統計的有意差がある遺伝子を同定する解析であり、「どの遺伝子が変動しているか」を明らかにするための基本的手法です。
• エンリッチメント解析 — 変動遺伝子群 (DEGs) がどのような生物学的プロセスやシグナル伝達経路に関与しているかを評価し、個々の遺伝子を上位の機能カテゴリに集約することで、対象となる事象の分子機構を俯瞰的に理解するための解析です。

本稿では、公開データベースである NCBI GEO (Gene Expression Omnibus) から取得した RNA-seq データを用いて、これらの解析を段階的に実施し、トランスクリプトーム解析の流れを見ていきます。

2. 実行環境とデータセット

2.1. 実行環境

本稿の解析は、以下のバージョンの Python および主要ライブラリを用いて実施しました。

2.2. RNA-seq データについて

RNA-seq で得られる発現量データには、マッピング直後の生カウントと、それを補正した正規化データがあります。これらは解析の目的や手法に応じて適切な入力形式が異なるため、それぞれの違いを理解しておくことが重要です。

生カウント (Raw Counts) は、各遺伝子にマッピングされたリード数そのものです。遺伝子長やライブラリ サイズの補正を行っていないため、そのままではサンプル間や遺伝子間の比較には適しません。一方、一部の統計的検定手法では、内部で正規化や分散推定を行う設計になっているため、生カウントを入力として使用します。

正規化データは、遺伝子長やライブラリ サイズの違いを補正し、発現量を比較しやすくした指標です。代表的なものとして以下があります。

• RPKM / FPKM — リード数を遺伝子長およびライブラリ サイズで補正した発現量指標です。遺伝子ごとに独立して補正が行われるため、同一サンプル内での遺伝子間比較には適していますが、サンプルごとにスケールが異なるため、サンプル間比較には注意が必要です。RPKM はシングルエンド、FPKM はペアエンドの RNA-seq に対応しています。
• TPM (Transcripts Per Million) — 遺伝子長で補正した後、サンプル全体で再スケーリングを行う指標です。このため各サンプルの総和が一定 (100 万) となり、サンプル間での比較に適した値が得られます。

正規化データは、機械学習や可視化など発現量の相対的な違いを比較する用途で用いられます。一方、本稿の DEG 解析では、DESeq2 [1] の入力として生カウントを使用しています。

2.3. 使用データセット

本稿では、NCBI GEO で公開されている GSE62944 データセットを使用しました。これは TCGA (The Cancer Genome Atlas) プロジェクトで収集された 24 種のがんの RNA-seq データを含む大規模データセットです [2]。

本稿では BRCA (乳がん) の腫瘍組織と正常組織を対象に、発現変動解析とエンリッチメント解析を行います。

以下のコードでは、生データから BRCA サンプルを抽出し、DEG 解析用のデータセットを準備しています。

import pandas as pd

# 解析対象とするがんの種類を指定します。今回は BRCA を指定します。
TARGET = "BRCA"

# ダウンロードした RNA-seq データのファイル名を指定します。
type_file_normal = "GSE62944_06_01_15_TCGA_24_Normal_CancerType_Samples.txt.gz"
type_file_tumor = "GSE62944_06_01_15_TCGA_24_CancerType_Samples.txt.gz"
counts_file_normal = "GSM1697009_06_01_15_TCGA_24.normal_Rsubread_FeatureCounts.txt.gz"
counts_file_tumor = "GSM1536837_06_01_15_TCGA_24.tumor_Rsubread_FeatureCounts.txt.gz"

# BRCA 患者の正常サンプルの ID を取得します。
df_type_normal = pd.read_csv(
    type_file_normal,
    sep=r"\s+",
    header=None,
    compression="gzip",
    engine="python"
)
normal_ids = df_type_normal[df_type_normal[1] == TARGET][0].tolist()

# BRCA 患者の腫瘍サンプルの ID を取得します。
df_type_tumor = pd.read_csv(
    type_file_tumor,
    sep=r"\s+",
    header=None,
    compression="gzip",
    engine="python"
)
tumor_ids = df_type_tumor[df_type_tumor[1] == TARGET][0].tolist()

print(f"BRCA サンプル数: 正常={len(normal_ids)}, 腫瘍={len(tumor_ids)}")

# 生カウント データを読み込みます。
df_counts_normal = pd.read_csv(
    counts_file_normal,
    sep=r"\s+",
    index_col=0,
    compression="gzip",
    engine="python"
)
df_counts_tumor = pd.read_csv(
    counts_file_tumor,
    sep=r"\s+",
    index_col=0,
    compression="gzip",
    engine="python"
)

# BRCA 患者の ID と一致するカウント データを抽出して、
# 正常サンプルと腫瘍サンプルを結合します。
df_brca_counts = pd.concat([
    df_counts_normal[normal_ids].T,
    df_counts_tumor[tumor_ids].T
])

print(f"遺伝子数: {df_brca_counts.shape[1]}")

# メタデータを作成します。`condition` 列に "normal" または "tumor" を格納します。
df_brca_meta = pd.concat([
    pd.DataFrame("normal", index=normal_ids, columns=["condition"]),
    pd.DataFrame("tumor", index=tumor_ids, columns=["condition"])
])

# Pickle ファイルとして保存します。
df_brca_counts.to_pickle("BRCA_counts.pkl")
df_brca_meta.to_pickle("BRCA_metadata.pkl")

3. 発現変動解析 (DEG 解析)

3.1. DEG 解析とは

RNA-seq データ解析の最初のステップとして広く行われるのが、DEG 解析 です。これは、2 つの条件間 (例: 腫瘍 vs 正常) で統計的に有意な発現変動を示す遺伝子を同定する手法であり、疾患に関連する遺伝子候補を抽出する出発点となります。

DEG 解析の結果からは、以下のような情報が得られます。

• 有意に変動した遺伝子のリスト — バイオマーカー候補や創薬ターゲットの探索に利用
• 発現変動の方向と大きさ (log2 fold change) — 発現の増減の程度を定量的に評価

本稿では PyDESeq2 ライブラリ [3] (DESeq2 の Python 実装、MIT License) を使用しました。DESeq2 は RNA-seq の生カウントを対象とした統計手法であり、各遺伝子の発現量を負の二項分布に基づいてモデル化し、ライブラリ サイズの違いを内部で補正したうえで検定を行います。

3.2. 解析結果

以下のコードでは、生カウント データと腫瘍または正常のメタデータを PyDESeq2 に入力し、遺伝子ごとの log2 fold change と補正後 p 値を算出しています。

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats

# フィルタリングに使用するリード数の下限を指定します。
MIN_READ_COUNT = 10

# 計算に使用する CPU の数です。
N_CPUS = 8

# 生カウント データとメタデータを読み込みます。
df_counts = pd.read_pickle("BRCA_counts.pkl")
df_meta = pd.read_pickle("BRCA_metadata.pkl")

# 全サンプルにおける合計リード数が `MIN_READ_COUNT` 未満の遺伝子は除外します (フィルタリング)。
genes_to_keep = df_counts.columns[df_counts.sum(axis=0) >= MIN_READ_COUNT]
df_counts = df_counts[genes_to_keep]

# DEG 解析用のオブジェクトを作成します。
inference = DefaultInference(n_cpus=N_CPUS)
dds = DeseqDataSet(
    counts=df_counts,
    metadata=df_meta,
    design="~condition",
    refit_cooks=True,
    inference=inference
)

# DEG 解析を実行し、統計量を計算します。
dds.deseq2()
stats = DeseqStats(
    dds,
    contrast=["condition", "tumor", "normal"],
    inference=inference
)
stats.summary()

# MA プロットを描画します。
stats.plot_MA(s=20, save_path="BRCA_MA_plot.png")

# 結果を CSV ファイルとして保存します。
stats.results_df.to_csv("BRCA_deg_results.csv")

DEG 解析の結果、腫瘍組織と正常組織の間で、16,690 個の遺伝子が有意 (補正後 p 値 < 0.05) な発現変動を示しました。

ここで、有意遺伝子数の解釈には注意が必要です。本データセットは症例数が非常に多く (腫瘍 1,119 例、正常 113 例) 検出力が高いため、ごく小さな発現差でも統計的に有意となりえます。すなわち、「統計的に有意」であることは「発現変動そのものが大きい」ことと同義ではありません。加えて、腫瘍組織と正常組織では細胞組成 (腫瘍細胞、脂肪細胞、免疫細胞など) そのものが大きく異なるため、観察された差のすべてが腫瘍細胞固有の発現変動を意味するわけではなく、組織全体としての細胞構成の違いも含まれている点に留意する必要があります。

まず、以下の MA プロットで発現変動の全体像を俯瞰します。横軸は各遺伝子の平均発現量 (log スケール) を表し、右にいくほど発現量の多い遺伝子です。縦軸は log2 fold change であり、正の値の遺伝子ほど腫瘍組織で高発現、負の値の遺伝子ほど正常組織で高発現であることを表します。

MA プロットから、本データセットでは有意遺伝子 (赤) の割合が非常に高いことが確認できます。ただし、補正後 p 値は「変動の確からしさ」を示すに過ぎず、個々の遺伝子がどの程度大きく変動しているかを表すものではありません。今回のデータセットは検出力が高いため、有意遺伝子の中には発現変動の小さい遺伝子も多く含まれることが想定され、有意性だけを見ていては個々の遺伝子の生物学的な重要性を判断できません。そこで、統計的有意性と発現変動の大きさの両方の分布を分析するために Volcano plot を確認します。

以下のコードにより、各遺伝子を 5 群に色分けした Volcano plot を描画できます。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# 補正後 p 値の閾値です。
PADJ_THRESHOLD = 0.05

# log2 fold change の絶対値の閾値です。
LOG2FC_THRESHOLD = 2.0

# DEG 解析の結果を読み込みます。
df = pd.read_csv("BRCA_deg_results.csv", index_col=0)

# 補正後 p 値または log2 fold change が欠損している遺伝子を除外します。
df = df.dropna(subset=["padj", "log2FoldChange"])

# -log10(補正後 p 値) を計算します。
# 補正後 p 値が 0 の行はゼロ除算を避けるため最小非ゼロ値で置き換えます。
min_nonzero_padj = df.loc[df["padj"] > 0, "padj"].min()
padj_clipped = df["padj"].clip(lower=min_nonzero_padj)
neg_log10_padj = -np.log10(padj_clipped)

# log2 fold change の値を取得します。
log2fc = df["log2FoldChange"]

# 有意性と発現変動の大きさで分類します。
# 有意かつ変動の大きい群は、さらに上昇または下降で分類します。
is_sig = df["padj"] < PADJ_THRESHOLD
is_large = log2fc.abs() > LOG2FC_THRESHOLD
mask_both_up = is_sig & is_large & (log2fc > 0)
mask_both_down = is_sig & is_large & (log2fc < 0)
mask_sig_only = is_sig & ~is_large
mask_large_only = ~is_sig & is_large
mask_none = ~is_sig & ~is_large

# 各群の遺伝子数を表示します。
print(f"総遺伝子数: {len(df):,}")
print(f"両方満たす (上昇) (padj<0.05 & log2FC>2): {mask_both_up.sum():,}")
print(f"両方満たす (下降) (padj<0.05 & log2FC<-2): {mask_both_down.sum():,}")
print(f"有意のみ (変動小): {mask_sig_only.sum():,}")
print(f"変動のみ (非有意): {mask_large_only.sum():,}")
print(f"両方満たさない: {mask_none.sum():,}")

# Volcano plot を描画します。
fig, ax = plt.subplots(figsize=(7, 6))

ax.scatter(
    log2fc[mask_none], neg_log10_padj[mask_none],
    s=6, c="#d0d0d0", alpha=0.5, linewidths=0,
    label=f"neither ({mask_none.sum():,})",
)
ax.scatter(
    log2fc[mask_sig_only], neg_log10_padj[mask_sig_only],
    s=6, c="#c0a8d8", alpha=0.6, linewidths=0,
    label=f"padj < 0.05 only ({mask_sig_only.sum():,})",
)
ax.scatter(
    log2fc[mask_large_only], neg_log10_padj[mask_large_only],
    s=6, c="#e0b070", alpha=0.85, linewidths=0,
    label=f"|log2FC| > 2 only ({mask_large_only.sum():,})",
)
ax.scatter(
    log2fc[mask_both_down], neg_log10_padj[mask_both_down],
    s=6, c="#2e86ab", alpha=0.85, linewidths=0,
    label=f"padj < 0.05 & log2FC < -2 (down, {mask_both_down.sum():,})",
)
ax.scatter(
    log2fc[mask_both_up], neg_log10_padj[mask_both_up],
    s=6, c="#d7263d", alpha=0.85, linewidths=0,
    label=f"padj < 0.05 & log2FC > 2 (up, {mask_both_up.sum():,})",
)

# 閾値線 (点線) を描画します。
ax.axvline(LOG2FC_THRESHOLD, color="black", linestyle="--", linewidth=0.8, alpha=0.6)
ax.axvline(-LOG2FC_THRESHOLD, color="black", linestyle="--", linewidth=0.8, alpha=0.6)
ax.axhline(-np.log10(PADJ_THRESHOLD), color="black", linestyle="--", linewidth=0.8, alpha=0.6)

ax.set_xlabel("log2 fold change (tumor / normal)")
ax.set_ylabel("-log10(adjusted p-value)")
ax.set_title("Volcano plot (BRCA tumor vs normal)")

# 凡例を右上に配置します。
ax.legend(loc="upper right", fontsize=9, framealpha=0.9)

# 図を保存します。
fig.tight_layout()
fig.savefig("BRCA_volcano_plot.png", dpi=150)
plt.close(fig)

print("BRCA_volcano_plot.png を保存しました。")

Volcano plot は横軸に log2 fold change (発現変動の方向と大きさ)、縦軸に -log10 (補正後 p 値) (統計的有意性) をとり、各遺伝子の発現変動を 1 つの図で俯瞰するプロットです。図の上方にある遺伝子ほど統計的に確かな変動であり、右側ほど腫瘍組織で高発現、左側ほど正常組織で高発現であることを意味します。

Volcano plot を見ると、遺伝子の大多数は |log2FC| の小さい中央の領域に位置しており、「統計的には有意だが発現変動そのものは小さい」遺伝子 (紫) が大多数 (有意遺伝子の 88%) を占めることが分かります。一方で、|log2FC| > 2 の領域にも相当数の遺伝子が分布しており、log2 fold change が ±6 を超える大きな発現変動を示す遺伝子も確認できます。閾値 |log2FC| > 2 を満たす遺伝子を上昇側 (赤) と下降側 (青) で比べると、上昇側のほうが log2 fold change の最大値が大きく、遺伝子数も上昇 1,241 個に対して下降 754 個と上昇側のほうが多い分布になっています。

では、具体的にどのような遺伝子が大きな発現変動を示しているのでしょうか。以下の表に、腫瘍組織で発現が上昇した遺伝子と減少した遺伝子の上位 10 個をそれぞれ示します。

3.3. 結果の考察

発現上昇遺伝子の上位には、CGA (糖タンパク質ホルモン α サブユニット) [4] や COL10A1 (X 型コラーゲン) [5]、CST1 (システイン プロテアーゼ阻害剤) [6] のように、乳がんでの高発現や予後不良との関連が報告され、予後バイオマーカーとしての可能性が示唆されている遺伝子が見られます。また、MAGEA3 や MAGEA6 はがん精巣抗原と呼ばれ、正常組織ではほとんど発現せず腫瘍で特異的に活性化される遺伝子群です [7]。このように、発現が上昇した遺伝子の上位には腫瘍細胞で特異的に発現が上昇するものが多く含まれていることがわかります。

一方、発現減少遺伝子の上位には、LEP (レプチン)、ADIPOQ (アディポネクチン)、PLIN1 (ペリリピン)、GPD1、CIDEC など、脂肪細胞に特徴的な遺伝子が集中しています。乳房の正常組織は脂肪細胞を豊富に含むことが知られており [8]、本結果は腫瘍組織における脂肪細胞成分の減少を反映していると考えられます。また、GLYAT (グリシン N-アシル転移酵素) は、乳がん細胞で発現が低下することが報告されています [9]。

このように、DEGs の上位遺伝子は「腫瘍組織に特徴的な発現変化」と「正常な脂肪組織成分の減少」という 2 つの生物学的変化を示唆しています。次章では、こうした個々の遺伝子レベルの変動をパスウェイに集約し、より俯瞰的な理解を目指します。

4. エンリッチメント解析

4.1. DEGs から生物学的意味を読み解く

DEG 解析によって「どの遺伝子が変動しているか」は分かりますが、多数の有意遺伝子を個別に解釈するのは現実的ではありません。そこで、DEGs が特定の生物学的パスウェイ (経路) に偏って存在するかどうかを統計的に検定するのがエンリッチメント解析です。

この解析により、個々の遺伝子レベルの変動を「代謝経路」「シグナル伝達」「細胞周期」などの上位の生物学的プロセスに集約でき、疾患がどのような分子メカニズムを介して進行しているかを俯瞰的に理解することができます。たとえば、がん関連の遺伝子が多数変動しているだけでなく、それらが特定の代謝経路に集中していれば、その経路が疾患に関連している可能性が示唆されます。

本稿では、Python の gget ライブラリ [10] (BSD 2-Clause License) を使用して WikiPathways [11] データベースに対してエンリッチメント解析を行いました。WikiPathways は、研究コミュニティが共同で編集・キュレーションを行うオープン ソースのパスウェイ データベースです。

次節では、あらかじめ選択した注目遺伝子セットが特定のパスウェイに偏って含まれるかを統計的に評価する over-representation analysis (ORA) を実施します。前章の Volcano plot で見たとおり、有意遺伝子には発現変動の小さいものが多数含まれます。パスウェイの偏りを生物学的に意味のある形で読み取るには、「統計的に確からしい」だけでなく「発現変動そのものが大きい」遺伝子に絞って入力することが望ましいと考えられます。そこで本解析では、有意性 (補正後 p 値 < 0.05) に加えて発現変動の大きさ (|log2 fold change| > 2) の条件も課し、両方を満たす遺伝子 1,995 個を ORA の入力としました。これは前章の Volcano plot で赤 (上昇) および青 (下降) に塗られた遺伝子群に相当します。

さらに、腫瘍で発現が上昇した遺伝子群と下降した遺伝子群では、背景にある生物学的意味が大きく異なることが予想されます。両者を同一セットにまとめて ORA を行うと、検出されたパスウェイが上昇・下降のどちらに由来するかが区別できなくなってしまいます。そこで本解析では、上昇遺伝子群 (1,241 個) と下降遺伝子群 (754 個) を分けて個別に ORA を実行し、腫瘍で活性化している経路と抑制されている経路をそれぞれ分析します。

なお、エンリッチメント解析ではバックグラウンド遺伝子リストの定義が結果に影響することに注意が必要です。本解析では、DEG 解析の対象となった遺伝子集合 (全サンプル合計リード数 10 以上の遺伝子) をバックグラウンドとしています。また、ORA は解析者が設定した閾値でフィルタリングした遺伝子セットを入力とする手法のため、閾値周辺の遺伝子の情報を十分に反映できません。今回は記事を簡潔にするために ORA のみを行っていますが、全遺伝子をランキングして直接解析する GSEA (Gene Set Enrichment Analysis) [12] は、閾値を設けずに経路レベルの偏りを評価する代表的な代替手法として広く用いられています。

4.2. 解析結果

以下のコードにより、有意かつ発現変動の大きい遺伝子を上昇または下降に分け、それぞれに対して WikiPathways に対するエンリッチメント解析を実行しています。

import gget
import pandas as pd

# 補正後 p 値の閾値です。
PADJ_THRESHOLD = 0.05

# log2 fold change の絶対値の閾値です。
LOG2FC_THRESHOLD = 2.0

# DEG 解析の結果を読み込みます。
df_deg = pd.read_csv("BRCA_deg_results.csv", index_col=0)

# バックグラウンド遺伝子のリストを取得します (DEG 解析の対象となった遺伝子集合です)。
background = df_deg.index.tolist()

# 有意かつ発現変動の大きい遺伝子を抽出し、上昇または下降に分けます。
mask_sig_large = (
    (df_deg["padj"] < PADJ_THRESHOLD) &
    (df_deg["log2FoldChange"].abs() > LOG2FC_THRESHOLD)
)
df_query = df_deg[mask_sig_large]
genes_up = df_query.index[df_query["log2FoldChange"] > 0].tolist()
genes_down = df_query.index[df_query["log2FoldChange"] < 0].tolist()

print(f"ORA 入力: 上昇 {len(genes_up):,} 遺伝子、下降 {len(genes_down):,} 遺伝子")

# 上昇側と下降側でそれぞれ WikiPathways に対してエンリッチメント解析を実行します。
result_up = gget.enrichr(
    genes=genes_up,
    database="WikiPathways_2024_Human",
    background_list=background
)
result_down = gget.enrichr(
    genes=genes_down,
    database="WikiPathways_2024_Human",
    background_list=background
)

# 結果をそれぞれ CSV ファイルとして保存します。
result_up.to_csv("enrichment_results_up.csv", index=False)
result_down.to_csv("enrichment_results_down.csv", index=False)

以下の表に、上昇遺伝子群および下降遺伝子群について、補正後 p 値の小さい順に上位 10 個の WikiPathways パスウェイを示します。補正後 p 値 < 0.05 を有意とすると、上昇側で 15 個、下降側で 56 個 のパスウェイが有意となりました。特に上昇側のヒット数が少ないため、より多面的に経路を読み取りたい場合は Reactome、KEGG、Gene Ontology など別のデータベースに対する ORA を併用するのも有効と考えられます。

上昇遺伝子群 (腫瘍で高発現) 上位 10 パスウェイ

下降遺伝子群 (腫瘍で低発現) 上位 10 パスウェイ

4.3. 結果の考察

上昇側と下降側で別々に ORA を行った結果、両者で検出されるパスウェイの傾向が大きく異なりました。以下ではそれぞれを分けて解釈します。

上昇遺伝子群では、Cell Cycle や G1/S 期移行、染色体分離制御といった細胞周期関連パスウェイが有意に濃縮されました。これらのパスウェイには細胞分裂の進行を直接制御する遺伝子群が含まれており、腫瘍細胞における増殖能の亢進を強く反映していると考えられます。また、炎症性サイトカインや MMP を含むパスウェイ (Overview Of Proinflammatory And Profibrotic Mediators) の濃縮は、腫瘍微小環境における炎症応答や細胞外マトリックス再構築の活性化を示唆しています。これらの特徴は、がんにおける代表的な性質である細胞増殖の活性化に加え、腫瘍微小環境の変化を介した浸潤・進展といった生物学的過程とも一致する結果といえます。

一方、下降遺伝子群では、Adipogenesis、Fatty Acid Omega Oxidation、PPAR signaling といった脂肪細胞の分化や脂質代謝に関連するパスウェイが顕著に濃縮されました。これらには脂肪細胞機能の維持に重要な遺伝子が含まれており、腫瘍組織における脂肪細胞由来シグナルの低下を示唆しています。本結果は、腫瘍組織における脂肪細胞成分の減少 (細胞組成の変化) や腫瘍細胞における代謝リプログラミングを反映していると考えられます。

このように、上昇遺伝子群と下降遺伝子群を分けて ORA を行うことで、「腫瘍組織で活性化される細胞増殖プロセス」と「腫瘍組織で失われる正常な脂肪細胞機能や代謝経路」という、2 種類の生物学的シグナルを切り分けて読み取ることができました。ただし、本解析のような bulk RNA-seq による解析では「サンプル間の細胞組成差」と「腫瘍細胞固有の代謝変化」を分離することができないため、両者を厳密に切り分けて議論するには single-cell RNA-seq による細胞種別の発現解析などが必要となります。

5. まとめ

本稿では、乳がん患者の RNA-seq データを題材に、トランスクリプトーム解析の基本的なワークフローを実演しました。

DEG 解析では、腫瘍組織でがん精巣抗原や乳がん関連マーカー (MAGEA3/6、COL10A1、CST1 など) が大きく上昇する一方、正常組織に豊富に含まれる脂肪細胞由来の遺伝子群 (LEP、ADIPOQ、PLIN1、GPD1、CIDEC など) が大きく下降することが確認できました。続くエンリッチメント解析では、上昇遺伝子群から細胞周期や染色分体分離制御などの細胞増殖関連経路の亢進が、下降遺伝子群から脂肪細胞分化や脂質代謝の減少が検出され、正常組織と腫瘍組織の対比が「細胞増殖プロセスの活性化」と「正常な脂肪細胞機能の低下」という 2 つの軸で読み解けることを示しました。

このように段階的に解析を進めることで、網羅的な RNA-seq データから生物学的情報を抽出し、疾患に特有の発現変化やその背後にある分子機構を体系的に把握することが可能となります。RNA-seq データは、本稿で取り上げた DEG 解析やエンリッチメント解析に加え、遺伝子発現パターンに基づくがん種分類や、機械学習による発現量予測など、多様な解析に展開することが可能です。

参考資料

• [1] M. I. Love, W. Huber, S. Anders, "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2", Genome Biology, 15:550, 2014. https://doi.org/10.1186/s13059-014-0550-8
• [2] M. Rahman et al., "Alternative preprocessing of RNA-Sequencing data in The Cancer Genome Atlas leads to improved analysis results", Bioinformatics, 31(22):3666-3672, 2015. https://doi.org/10.1093/bioinformatics/btv377
• [3] B. Muzellec, M. Teleńczuk, V. Cabeli, M. Andreux, "PyDESeq2: a python package for bulk RNA-seq differential expression analysis", Bioinformatics, 39(9):btad547, 2023. https://doi.org/10.1093/bioinformatics/btad547
• [4] J. Zhou et al., "Glycoprotein hormone α-subunit promotes cell proliferation and tumorigenesis in breast cancer", Oncology Letters, 23(5):142, 2022. https://doi.org/10.3892/ol.2022.13263
• [5] M. Zhang et al., "Bioinformatics analysis of prognostic significance of COL10A1 in breast cancer", Bioscience Reports, 40(2):BSR20193286, 2020. https://doi.org/10.1042/BSR20193286
• [6] Dn. Dai et al., "Elevated expression of CST1 promotes breast cancer progression and predicts a poor prognosis", Journal of Molecular Medicine, 95:873-886, 2017. https://doi.org/10.1007/s00109-017-1537-1
• [7] A. Alsalloum, J. A. Shevchenko, S. Sennikov, "The Melanoma-Associated Antigen Family A (MAGE-A): A Promising Target for Cancer Immunotherapy?", Cancers, 15(6):1779, 2023. https://doi.org/10.3390/cancers15061779
• [8] C. Wu et al., "Cancer-Associated Adipocytes and Breast Cancer: Intertwining in the Tumor Microenvironment and Challenges for Cancer Therapy", Cancers, 15(3):726, 2023. https://doi.org/10.3390/cancers15030726
• [9] X. Tian et al., "Downregulation of GLYAT Facilitates Tumor Growth and Metastasis and Poor Clinical Outcomes Through the PI3K/AKT/Snail Pathway in Human Breast Cancer", Frontiers in Oncology, 11:641399, 2021. https://doi.org/10.3389/fonc.2021.641399
• [10] L. Luebbert, L. Pachter, "Efficient querying of genomic reference databases with gget", Bioinformatics, 39(1):btac836, 2023. https://doi.org/10.1093/bioinformatics/btac836
• [11] A. Agrawal et al., "WikiPathways 2024: next generation pathway database", Nucleic Acids Research, 52(D1):D679-D689, 2024. https://doi.org/10.1093/nar/gkad960
• [12] A. Subramanian et al., "Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles", PNAS, 102(43):15545-15550, 2005. https://doi.org/10.1073/pnas.0506580102