肺癌分子標的薬治療で代表的な遺伝子異常である融合遺伝子EML4-ALKの転写予測を行った。正常肺ではALKの発現はないが、EML4-ALK陽性肺癌では融合遺伝子の mRNAの発現がある。ヒトゲノム領域と組織を指定したAlphaGenomeの「予測」は、データが学習データに含まれているので、「予測」ではなく「検索」に近い。染色体転座のデータは、一部の培養細胞株を除いて学習データには含まれていないので、AlphaGenomeの予測能評価のテストに適している。
融合遺伝子周辺のゲノム配列作成
染色体切断点の検索
融合遺伝子の染色体切断点（５’及び３’側）はChimerDBのChimerSeqで検索する。ChimerSeqはヒトゲノム参照配列GRCh37/hg19に準拠しているため、AlphaGenomeが採用しているGRCh38/hg38に座標変換する必要がある。
ChimerSeqの検索結果が、５’切断点がchr2 : 42472827(+) 、３’切断点が chr2 : 29446394(-)の場合、次の形に整形してLiftOverでGRCh37/hg19からGRCh38/hg38に変換する。ここで(+) はmRNA鎖がforward方向、(-)はreverse方向を指している。
入力：
chr2 42472827 42472828 +
chr2 29446394 29446395 -
出力：
chr2 42245687 42245688 +
chr2 29223528 29223529 -
注：第４カラムは何でもよい。
５’切断点がchr2 : 42245687(+) 、３’切断点が chr2 : 29223528(-)に変換された。
配列取得と配列連結
mRNAの方向に注意して取得するゲノム配列領域を決定する。症例によって例外があるが、EML４−ALKの場合EML４はforward、ALKはreverseになる。次の領域で融合部位周辺１MBになる。
５’（EML4）側：chr2 41745687 42245687
３’（ALK）側：chr2 28723528 29223528
curlでNCBIから配列を取得する。
curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_000002.12&seq_start=41745687&seq_stop=42245687&strand=1&rettype=fasta&retmode=text" -o chr2_41745687_42245687_plus.fa
curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_000002.12&seq_start=28723528&seq_stop=29223528&strand=2&rettype=fasta&retmode=text" -o chr2_28723528_29223528_minus.fa
次のシェルスクリプトで連結する。
(
grep -v '^>' chr2_41745687_42245687_plus.fa | tr -d '\n'
grep -v '^>' chr2_28723528_29223528_minus.fa | tr -d '\n'
) > chr2_41745687_29223528.txt
注：AlphaGenomeの入力配列は”ACGTN"のみ可。
chr2_41745687_29223528.txtが融合遺伝子配列：これをAlphaGenomeの入力配列とする。
AlphaGenomeによるEML4-ALK転写予測（RNA-Seq予測）
パイソンスクリプト
次のスクリプトでRNA−Seqの予測を行う。組織は肺。対照のEML4, ALKも融合遺伝子配列と同様塩基配列を下記スクリプトを使って直接入力で分析した。ゲノム領域指定は使っていない。出力は１塩基位置ごとのRNA-Seqリードで、単位はreads per million (RPM)。下記スクリプトではdata.txtに出力される。
融合遺伝子配列についてはsplice siteの予測も行った（スクリプトは略）。
import os
from dotenv import load_dotenv
import matplotlib.pyplot as plt
from alphagenome.models import dna_client
from alphagenome.data import genome
from alphagenome.visualization import plot_components
import numpy as np
# --- 設定 ---
load_dotenv()
API_KEY = os.getenv("ALPHAGENOME_API_KEY")
dna_model = dna_client.create(API_KEY)
f = open('chr2_41745687_29223528.txt', 'r')
CORE_SEQ = f.read()
f.close()
# ---分析対象配列の長さの調整：１MB=1048576（２の20乗）のため ---
TARGET_LEN = 1048576
if len(CORE_SEQ) > TARGET_LEN:
raise ValueError("入力シーケンスがターゲット長を超えています。")
padded_sequence = CORE_SEQ.center(TARGET_LEN, 'N')
# --- 合成 interval ---
synthetic_interval = genome.Interval(
chromosome="chrSynthetic",
start=0,
end=len(padded_sequence),
)
region_2kb=synthetic_interval.resize(dna_client.SEQUENCE_LENGTH_2KB)
# --- 推論 ---
output = dna_model.predict_sequence(
sequence=padded_sequence,
requested_outputs=[dna_client.OutputType.RNA_SEQ],
ontology_terms=["UBERON:0002048"],
interval=synthetic_interval,
)
# --- データ保存 ---
rna_seq = output.rna_seq
print(rna_seq.values.shape)
np.savetxt('data.txt', rna_seq.values)
# --- プロット ---
plot_components.plot(
components=[plot_components.Tracks(output.rna_seq)],
interval=output.rna_seq.interval,
annotations=[plot_components.IntervalAnnotation([region_2kb], alpha=0.8)],
)
plt.show()
グラフ表示
出力のグラフ表示。融合部位は中央の縦線。パネルは上からEML4(RNA-Seq),ALK(RNA-Seq),EML4-ALK(RNA-Seq),EML4-ALK(Splice Site)。それぞれのパネルは４つのトラックから構成されている。EML4(RNA-Seq)は上から polyA RNA・forward strand, total RNA・forward, polyA RNA・reverse, total RNA・reverse。ALK (RNA-Seq)は上から polyA RNA・reverse, total RNA・reverse, polyA RNA・forward, total RNA・forward。EML4-ALK(RNA-Seq)は、上の２つのトラックはEML4部分がforwardでALK部分がreverse、下の２つのトラックはEML4部分がreverseでALK部分がforward; 上から polyA RNA, total RNA, polyA RNA, total RNA。EML4-ALK(RNA-Seq)は、すべてEML4部分がforwardでALK部分がreverse；上からpolyA RNA, total RNA, splice site donor, splice site acceptor。
２番目のパネルALK (RNA-Seq)と３番目のパネルEML4-ALK(RNA-Seq)の融合部位（中央縦線）周辺に注目する。ALK (RNA-Seq)では全くRNA発現がないが、EML4-ALK(RNA-Seq)では融合部位の下流のエクソンでRNA発現がある。
発現量算出
EML4-ALK(RNA-Seq)について発現しているエクソンの部分に絞ってデータ解析を行った。
今回は予備的な解析で、上のグラフを見ながら適当に分析した。塩基位置520000-620000の範囲で0.2以上のデータを抽出して平均と標準偏差を求めた。融合遺伝子で発現するALKエクソンはこの範囲内にあり、その発現量は0.2以上であるためだ。計算結果は、平均0.498 RPM、標準偏差0.173、該当塩基部位は2206だった。
まとめ
肺組織についてALKではmRNA発現は認められないが、EML4-ALKで発現しているので、これまでの研究で観察されている事象が再現されている。一応初期条件は達成できた状態だが、どのような応用が可能かどうか、少し考える必要がある。
補足：APIキーをパイソン仮想環境下の環境変数として設定する
上記のパイソンスクリプトでは、AlphaGenome API Keyをpython-dotenv moduleを使って環境変数として設定している。方法は次のサイトを参照。
Pythonの venv 仮想環境内で環境変数を恒常的に設定する方法 #API - Qiita
