動機
自分で開発しているから記事を書こうと思う。
自分で使いたいと思えないツールなんて、そもそもツールではない。
今後も、自分で作ったツールだからこそ、自分が使いたいと思えるような人間でありたい。
Javaをベースにして作ったツールなのだが、Pythonで使えないと、機械学習でも、画像処理でも、自分でも、使うタイミングに苦慮する（言ってしまえば、Colabで使えないと使う機会が激減する）＋、他の誰にも使ってもらえないことに気がついたので、Pythonラッパーを作った。その検証ログ。
PyPIからのインストール
!pip install radiomicsj
※バージョンは2.1.18以上にする。
※（重要）Colabの場合、インストール後、セッションを再起動する。
HowTo
何か、テーマがあったほうが面白いので、ストーリーを設定しようと思う。
以降、仮説になるが、概念的には、概ね、筋は通っているだろう。
（仮題）「メチル化特有のRadiomicsシグナルを抽出・証明する」
本記事のテーマは、思い切って「がん」を扱うことにした。
がんは、私が生きているうちに克服されるのではないか、と漠然とではあるが期待している。そう考える理由は、AIの発展により、医師や医科学研究者、現場のメディカルスタッフの思考力が底上げされ、研究や臨床のアプローチが飛躍的に進化する時代が到来すると考えているからだ（メディカルスタッフのやる気ではなく、ルーティーンの質が上がる）。
しかし、がんという疾患は依然として非常に難解であることに変わりはない。細胞の暴走は極めてダイナミックであり、時々刻々と浸潤や転移、微小環境の改変（血管新生など）を引き起こしながら進化していく。その実態は、理論以外には掴みどころがない（今どうなっているなど誰にもわからない）。
そのような中、近年のゲノミクス研究の進展により、がんの分子生物学的なメカニズムが次々と解明されてきている。中でも、がん細胞におけるDNAのメチル化（DNA methylation）は、がんの発生や進行を理解する上で非常に重要なテーマとなっている。DNAメチル化は、DNAの塩基配列を研究する分野「エピジェネティクス」の代表的なメカニズムである。正常細胞とがん細胞とでは、このメチル化のパターンが大きく異なり、がん細胞においては主に以下の「2つの異常」が引き起こされるとされている。
-
がん抑制遺伝子プロモーター領域の「高メチル化」（ブレーキの故障）
細胞には本来、異常な増殖を抑え、がん化を防ぐための「がん抑制遺伝子」が備わっている。正常な細胞では、この遺伝子の発現を制御するスイッチ部分（プロモーター領域）はメチル化されておらず、遺伝子が常に「ON」の状態で細胞を守っていると考えられている。しかし、一部のがん細胞では、この領域に異常なメチル化（高メチル化）が生じる。これによりがん抑制遺伝子のスイッチが「OFF」にされ、細胞増殖のブレーキが効かなくなり、がん化が促進されてしまうことがある。 -
ゲノム全体の「低メチル化」（ゲノムの不安定化とアクセルの暴走）
特定の遺伝子領域が高メチル化を示す一方で、ゲノム全体（大部分を占める反復配列など）のレベルで見ると、がん細胞では逆にメチル化が顕著に減少する現象が起きることが知られている。正常な細胞では、これらの領域がしっかりとメチル化されることで染色体の構造を安定させ、不要な遺伝子の発現を抑制している。しかし、特異ながん細胞では、低メチル化状態となることがある。これにより染色体が不安定化して遺伝子変異が誘発されやすくなるほか、本来は抑制されているべき「がん遺伝子」が誤って活性化し、細胞増殖のアクセルが暴走する要因となることがある。
こうしたがん細胞特有のメチル化パターンの異常は、現在の臨床現場においても重要なターゲットとなっている。
例えば、がん細胞から血中に遊離した異常メチル化DNAを検出することで、低侵襲にがんの早期発見や性質診断を行う技術（リキッドバイオプシーなど）の開発が進められている。また、がん細胞に付加された異常なメチル基を取り除く「DNAメチル化阻害薬」を用いたエピジェネティック治療も開発されており、骨髄異形成症候群や急性骨髄性白血病などの血液がんにおいてすでに実用化されている。このように、DNAメチル化の異常は、薬理学的な介入によって、元の正常な状態へ戻し得る「可逆的である可能性が高い」変化である。
そこで本記事では、こうしたミクロな分子生物学的変化（メチル化パターンの表現型）が、マクロな医用画像のテクスチャ情報にも反映されるのではないか、という仮説（Radiogenomicsの概念）に基づいて、擬似的なデータを用いてシナリオを展開する。具体的には、画像から得られる定量的な特徴情報（ヒストグラムやテクスチャの空間的配置など）を機械学習モデルに入力し、それが特定のメチル化状態の予測・判断に有用であるかを検証するプロセスを扱う。
Radiomics計算ツール「RadiomicsJ」
今回の例では、RadiomicsJを使う。
!pip install radiomicsj
※念の為、Colab環境では、ランタイムを再起動する。
疑似デジタルファントム
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
def generate_digital_phantom(size=256, radius=80, is_2hg=False, seed=None):
"""
MRIの擬似デジタルファントム画像を作成する関数
"""
# 【追加】シード値が指定された場合、乱数生成器を固定する
if seed is not None:
np.random.seed(seed)
# 背景（空気）を0とするベースキャンバスを作成
img = np.zeros((size, size))
# 紙コップ（関心領域：ROI）の円形マスクを作成
y, x = np.ogrid[-size/2:size/2, -size/2:size/2]
mask = x**2 + y**2 <= radius**2
# 【基本設定】アガロースゲル＋マイクロビーズのベース信号
base_signal = 100.0 # T2/FLAIRなどのベース輝度
noise_sigma = 5.0 # ビーズによる不均一性（テクスチャのベース）
if is_2hg:
# 【2-HG群の設定】
# 2-HGによるpH変化や水分子の相互作用が、ピクセル値の平均や
# テクスチャ（ばらつき）に微小な変化を与えるという仮説のシミュレーション
base_signal = 105.0 # 平均輝度がわずかに上昇すると仮定（+5%）
noise_sigma = 7.5 # テクスチャの複雑さ（ばらつき）が増すと仮定（+50%）
# マスク内にベース信号を配置
img[mask] = base_signal
# 正規分布ノイズを追加して「組織のザラつき（テクスチャ）」を表現
noise = np.random.normal(0, noise_sigma, (size, size))
img = img + (noise * mask)
# 1.5T MRIの空間分解能（ぼやけ）を擬似的に再現するためガウシアンフィルターを適用
img = gaussian_filter(img, sigma=1.0)
# マスク外（空気）の信号をゼロに戻し、背景ノイズを少し足す
background_noise = np.random.normal(0, 1.0, (size, size))
img[~mask] = background_noise[~mask]
return img, mask
# --- メイン処理 ---
# 1. コントロール群と実験群の画像を生成
# 【変更】それぞれ異なるシード値を設定し、独立したスキャンを再現しつつ毎回同じ結果になるようにする
phantom_control, roi_mask = generate_digital_phantom(is_2hg=False, seed=76)
phantom_2hg, _ = generate_digital_phantom(is_2hg=True, seed=7)
# 2. 画像の描画・比較
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# コントロール群の表示(メチル化なし)
im0 = axes[0].imshow(phantom_control, cmap='gray', vmin=70, vmax=130)
axes[0].set_title('Control Phantom\n(Gel + Beads)')
axes[0].axis('off')
fig.colorbar(im0, ax=axes[0], fraction=0.046, pad=0.04)
# 2-HG群の表示(高メチル化)
im1 = axes[1].imshow(phantom_2hg, cmap='gray', vmin=70, vmax=130)
axes[1].set_title('2-HG Phantom\n(Gel + Beads + 2-HG)')
axes[1].axis('off')
fig.colorbar(im1, ax=axes[1], fraction=0.046, pad=0.04)
# 差分画像の表示（肉眼では見えにくい差を強調）
difference = phantom_2hg - phantom_control
im2 = axes[2].imshow(difference, cmap='bwr', vmin=-10, vmax=10)
axes[2].set_title('Difference Map\n(2-HG - Control)')
axes[2].axis('off')
fig.colorbar(im2, ax=axes[2], fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()
# 3. Radiomics解析用に出力（オプション：必要に応じてTIFF等で保存）
# import cv2
# cv2.imwrite('phantom_control.tif', phantom_control.astype(np.float32))
# cv2.imwrite('phantom_2hg.tif', phantom_2hg.astype(np.float32))
データセットの作成（各10枚）
import os
import cv2
import pandas as pd
# ---------------------------------------------------------
# 2. データセットの作成（各クラス10枚）
# ---------------------------------------------------------
os.makedirs('dataset', exist_ok=True)
n_samples = 10
images_3d = []
masks_3d = []
labels = []
print("画像データを生成しています...")
# マスクラベルは1とすること
# Control群 (ラベル 0)
for i in range(n_samples):
img, mask = generate_digital_phantom(is_2hg=False, seed=i)
# radiomicsj 用に (z, y, x) の3次元配列へ変換
# mask は boolean から数値(1/0)の uint8 に変換
img_3d = np.expand_dims(img, axis=0)
mask_3d = np.expand_dims(mask.astype(np.uint8), axis=0)
images_3d.append(img_3d)
masks_3d.append(mask_3d)
labels.append(0)
# ローカル保存用（オプション）
cv2.imwrite(f'dataset/control_{i}.tif', img.astype(np.float32))
# 2-HG群 (ラベル 1)
for i in range(n_samples):
img, mask = generate_digital_phantom(is_2hg=True, seed=i)
img_3d = np.expand_dims(img, axis=0)
mask_3d = np.expand_dims(mask.astype(np.uint8), axis=0)
images_3d.append(img_3d)
masks_3d.append(mask_3d)
labels.append(1)
# ローカル保存用（オプション）
cv2.imwrite(f'dataset/2hg_{i}.tif', img.astype(np.float32))
RadiomicsJによる特徴量抽出
from radiomicsj import calculate_features
print("Radiomics特徴量の抽出を開始します。少し時間がかかります...")
df_list = []
for i, (img_3d, mask_3d) in enumerate(zip(images_3d, masks_3d)):
# radiomicsj を用いた計算
df = calculate_features(img_3d, mask_3d, spacing=(1.0, 1.0, 1.0))
df_list.append(df)
print("Now...", i)
# 各画像の行データ（DataFrame）を縦に結合
df_features = pd.concat(df_list, ignore_index=True)
df_features['Label'] = labels
# 不要な列をDropします
# 今回、マスクは不変のため、形態的な特徴量は除外します。
prefixes_to_drop = ('OperationalInfo_', 'Diagnostics_', 'Morphology_', 'Fractal_')
cols_to_drop = [col for col in df_features.columns if col.startswith(prefixes_to_drop)]
df_features.drop(columns=cols_to_drop, inplace=True)
print(f"メタデータ・形状特徴量など {len(cols_to_drop)} 列を削除しました。")
# NaN（欠損値）のチェックと処理
nan_count = df_features.isnull().sum().sum()
if nan_count > 0:
print(f"警告: 特徴量データに {nan_count} 個の NaN が含まれています。0で補完します。")
df_features.fillna(0, inplace=True)
else:
print("NaNチェック: 欠損値は検出されませんでした。")
ロジスティック回帰での学習と重要変数の取得
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
numeric_features = df_features.select_dtypes(include=[np.number])
X = numeric_features.drop(columns=['Label'], errors='ignore')
y = df_features['Label']
# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# ロジスティック回帰モデルの構築と学習
model = LogisticRegression(random_state=42, max_iter=1000)
model.fit(X_scaled, y)
coefficients = model.coef_[0]
# 係数の絶対値が大きい順にインデックスを取得
sorted_indices = np.argsort(np.abs(coefficients))[::-1]
print("--- 最も影響力が大きい特徴量 Top 5 ---")
for i in range(5):
idx = sorted_indices[i]
f_name = X.columns[idx]
f_coef = coefficients[idx]
print(f"{i+1}位: {f_name} (係数: {f_coef:.4f})")
print("--------------------------------------")
結果
--- 最も影響力が大きい特徴量 Top 5 ---
1位: IntensityVolumeHistogram_IntensityAtVolumeFraction90 (係数: 0.0588)
2位: IntensityBasedStatistical_Percentile10 (係数: 0.0588)
3位: GLDZM_GrayLevelVariance (係数: -0.0586)
4位: GLSZM_GrayLevelVariance (係数: -0.0586)
5位: IntensityBasedStatistical_Median (係数: 0.0583)
--------------------------------------
信号系は、デジタルファントムを作っている時点で、恣意的なので、ここでは除外します。
可視化マップの作成と表示（グローバルコントラスト適用・自動ラベル）
※計算結果にコントラストがつかず、変化が見られないときは、一旦、ランタイムを再起動して、再計算してください。
from radiomicsj.features import generate_feature_map
from radiomicsj import GLSZM, GLCM, GLDZM
print(f"\nFeature Mapを生成しています（全{n_samples * 2}枚）。少し時間がかかります...")
settings = {
"label": 1,
"useBinCount": True,
"nBins": 16,
# "useBinWidth": True,
# "binWidth": 2.0,
}
idx_control = 0
idx_2hg = n_samples
target_class = GLSZM # モデルの結果に合わせて変更してください（例: GLSZM）
target_feature = GLSZM.GrayLevelVariance # モデルの結果に合わせて変更してください（例: GLSZM.GrayLevelVariance）
# target_class = GLCM # モデルの結果に合わせて変更してください（例: GLSZM）
# target_feature = GLCM.JointEntropy # モデルの結果に合わせて変更してください（例: GLSZM.GrayLevelVariance）
top_feature_name = "GLSZM_GrayLevelVariance"
target_z = 0
fmap_slices_control = []
fmap_slices_2hg = []
# 10ペア分のマップを計算
for i in range(n_samples):
print(f" 計算中... ペア {i+1}/{n_samples}")
# Control群の計算
fmap_c = generate_feature_map(
image_np=images_3d[i],
mask_np=masks_3d[i],
mask_label=1,
spacing=(1.0, 1.0, 1.0),
feature_class=target_class,
feature_id=target_feature,
settings=settings,
filter_size=9,
d2_mode=True,
stride=2,
slice_idx=-1
)
slice_c = fmap_c[target_z, :, :].copy()
slice_c[masks_3d[i][target_z, :, :] == 0] = np.nan
fmap_slices_control.append(slice_c)
# 2-HG群の計算
idx_2hg = i + n_samples
fmap_2hg = generate_feature_map(
image_np=images_3d[idx_2hg],
mask_np=masks_3d[idx_2hg],
mask_label=1,
spacing=(1.0, 1.0, 1.0),
feature_class=target_class,
feature_id=target_feature,
settings=settings,
filter_size=9,
d2_mode=True,
stride=2,
slice_idx=-1
)
slice_2hg = fmap_2hg[target_z, :, :].copy()
slice_2hg[masks_3d[idx_2hg][target_z, :, :] == 0] = np.nan
fmap_slices_2hg.append(slice_2hg)
print("描画を準備しています...")
# 全20枚のマップデータから、外れ値を除外した統一スケール(1%〜99%)を計算
all_fmaps = np.stack(fmap_slices_control + fmap_slices_2hg)
valid_min = np.nanpercentile(all_fmaps, 1)
valid_max = np.nanpercentile(all_fmaps, 99)
feature_display_name = top_feature_name.replace("original_", "")
# 10行 x 2列の巨大な図面を作成
fig, axes = plt.subplots(n_samples, 2, figsize=(12, 4 * n_samples))
for i in range(n_samples):
# Control群の表示 (左列)
axes[i, 0].imshow(images_3d[i][target_z, :, :], cmap='gray', vmin=70, vmax=130)
# cmap='jet' を cmap='jet_r' に変更して色を反転
im_map = axes[i, 0].imshow(fmap_slices_control[i], cmap='jet_r', alpha=0.6, vmin=valid_min, vmax=valid_max)
axes[i, 0].set_title(f"Control Phantom (Pair {i+1})\n{feature_display_name} Map")
axes[i, 0].axis('off')
# 2-HG群の表示 (右列)
idx_2hg = i + n_samples
axes[i, 1].imshow(images_3d[idx_2hg][target_z, :, :], cmap='gray', vmin=70, vmax=130)
# こちらも cmap='jet_r' に変更
axes[i, 1].imshow(fmap_slices_2hg[i], cmap='jet_r', alpha=0.6, vmin=valid_min, vmax=valid_max)
axes[i, 1].set_title(f"2-HG+ Phantom (Pair {i+1})\n{feature_display_name} Map")
axes[i, 1].axis('off')
# 図面全体の右側に1つのカラーバーを配置
plt.tight_layout()
fig.subplots_adjust(right=0.88)
cbar_ax = fig.add_axes([0.90, 0.15, 0.03, 0.7])
cbar = fig.colorbar(im_map, cax=cbar_ax)
cbar.set_label(f'{feature_display_name} Value', fontsize=14)
ticks = np.linspace(valid_min, valid_max, num=5)
cbar.set_ticks(ticks)
cbar.set_ticklabels([f"{val:.4g}" for val in ticks])
plt.show()
結果
左はコントロール、右は2-HG(メチル化)の疑似データによる結果となっている。
意図して作られているデジタルファントムなので当然なのだが、コントロールのほうが疑似2-HGに比べて若干、緑のつぶが小さくなっている傾向が伺え、逆に、2-HGでは大きくなっている。
微細な違いだが、このような違いを機械的・ミクロ-マクロの境界を超えた傾向として捉えて定量化することが、Radiomics分野の得意技なのである。
おわり
このように、肉眼では捉えることが難しい信号パターンを、様々な特徴を用いて定量化し、差異を見出すことには、意味があるのだろうと思う。
Stay visionary