この記事は、頭の中味を整理するために、下書きと生成AIによる修正のサイクルを繰り返して執筆されたハイブリッド記事です。AIの出力を使用していますが「バイブ記事」ではなく、人間の手間をそれなりにかけて作成しています。

ソフトウェアの外で

生命科学のデータは、観測の対象を必要とします。

実験的な研究では、細胞株、オルガノイド、マウス、患者由来 xenograft（PDX）などが使われます。臨床研究では、医療機関で収集された血液、正常組織、腫瘍組織、病理標本、治療前後の検体などが対象になります。これらは、採血、生検、手術などによって確保されます。

サンプルは、そのままデータになるわけではありません。DNA抽出、RNA抽出、タンパク質抽出、ライブラリ調製、染色、細胞分取、質量分析、顕微鏡撮像などの操作を通じて、測定可能な対象に変換されます。

観測対象を設計する作業もあります。実験では、薬剤処理、遺伝子改変、培養条件、動物実験などを通じて、比較可能な条件を作ります。臨床研究では、患者群、治療前後、奏効例と非奏効例、原発巣と転移巣など、現実の制約のなかで比較可能な条件を設定します。

臨床研究では、現病歴、既往歴、薬剤情報といった電子カルテ情報や、レントゲン、エコー、CT、MRI、内視鏡といった医用画像も使われます。ただし、これらの情報はゲノムデータと同じ形式で扱えるわけではありません。研究データとして統合するには、別の作業が必要になります。患者本人の同意や倫理審査も必要です。

リサーチクエスチョンに基づいて研究内容を決定すべきだと言われます。一方で、研究資源には制約があります。利用できる検体、測定技術、研究費、人員、施設、倫理的条件によって、立てられる問いは影響を受けます。理想的な問いから研究が始まる場合だけでなく、観測できる内容から問いが選ばれる場合もあります。

病理医や検査技師は、腫瘍含有率、壊死、正常組織の混入、FFPEの品質などを確認し、病理診断や検体の品質を判断します。バイオバンクのように、検体を管理・保存する部門もあります。

その後、DNAシーケンシング、RNAシーケンシング、シングルセルシーケンシング、ロングリードシーケンシング、フローサイトメトリー、免疫染色、空間トランスクリプトーム、プロテオミクス、病理画像、質量分析などが行われます。

次世代シーケンサーのゲノムデータ

全エクソームシーケンスや全ゲノムシーケンスでは、DNA上の変異、コピー数変化、構造変異を調べます。対象遺伝子を絞った遺伝子パネルシーケンスもあります。RNA-Seq は遺伝子発現を測定し、scRNA-seq は細胞ごとの状態や細胞集団の違いを調べます。ChIP-Seq は転写因子やヒストン修飾がゲノム上のどこに結合しているかを調べ、ATAC-Seq はクロマチンの開きやすさを測ります。Hi-C はゲノムの三次元構造や、離れた領域同士の接触を調べます。ロングリードシーケンシングは、構造変異、反復配列、ハプロタイプ、メチル化などを調べるために使われます。微生物群集の塩基配列を読む16S rRNA解析もあります。

バルク RNA-Seq や WES は、外注も普及しています。一方で、上記のシーケンス手法のほとんどは、実施できる施設と人材が限られます。

次世代シーケンサーを使わない実験手技も併用されます。PCR やリアルタイムPCRは、特定の配列や遺伝子発現を確認するために使われます。ウェスタンブロッティングは特定のタンパク質の存在量やサイズを確認し、ELISAはタンパク質やサイトカインなどを定量します。フローサイトメトリーは、細胞表面マーカーや細胞内タンパク質を使って細胞集団を分類します。

免疫染色や病理画像は、組織の形態、タンパク質発現、細胞の配置を見るために使われます。タンパク質量や翻訳後修飾を測るプロテオーム解析や、代謝産物を定量するメタボローム解析もあります。またタンパク質の折りたたまれた立体構造を予想したり機能を解明する研究する分野もあります。

一次解析の世界

日本語では、シーケンスデータにおける一次解析は、FASTQ から SAM/BAM/CRAM、VCF/BCF、遺伝子発現カウント行列などを作る工程を指すことが多いと思います。

英語圏では、シーケンサーを作る企業や装置側の文脈もあるため、電気信号から FASTQ を生成する段階を一次解析と呼ぶ場合があります。一方、日本の研究現場では、その段階は機種やメーカーの標準パイプラインとして扱われ、あまり意識されないこともあります。Nanopore のpod5のように、生データをユーザーが触れる形式で保存するものもありますが、すべてのシーケンサーで同じように扱えるわけではありません。

DNAシーケンシングでは、bwa、minimap2 などを用いてリードをリファレンスゲノムにマッピングします。RNA-Seq では STAR、HISAT2、Salmon、kallisto、featureCounts などが使われます。Cell Rangerのようなプロプライエタリなツールの内部でも、STAR 系の処理が使われています。がんゲノム解析では、マッピング後に GATK Mutect2 などで腫瘍・正常ペアの比較を行い、体細胞変異を VCF/BCF として出力します。生殖細胞系列の変異コールでは、GATK HaplotypeCaller が使われることが多いです。

この工程には計算資源が必要になります。htop などでマシンの様子を見ると、CPU の複数コアが長時間使われていることがわかります。複数検体を処理する場合、小規模ならワークステーション、自作PC、解析サーバーが使われます。規模が大きくなると、HPC やクラウドコンピューティングが使われます。

一次解析では、測定装置から出てきた生データを、標準的な座標系とファイル形式に変換します。シーケンスであれば、リードをリファレンスゲノム座標に対応づけます。データは専用のバイナリ形式で扱われることがあります。これらのファイルは、ランダムアクセスが可能なようにBGZF圧縮され、インデックスを持ち、専用のツールで処理されます。

ここでは Unix のコマンドラインツール、パイプライン処理、ジョブスケジューラ、ファイル形式、I/O、並列処理、ログの読解が重要になります。最近では、ログの読解やエラー原因の切り分けにAIを使える場面もあります。ただし、ファイル名、サンプル対応、リファレンス、バージョン、実行条件は人間が確認する必要があります。

一次解析を運用する人は、Linux、シェルスクリプト、HPC、Docker や Apptainer、Nextflow や Snakemake、リファレンスゲノム、アノテーション、QC、ログ確認、ファイル管理に習熟していることが多いです。重要なのは、ツールを安全に動かし、サンプルの取り違え、低いマッピング率、異常なカバレッジ、バッチ効果、リファレンスの不整合などに気づくことです。最近ではnf-coreのようなワークフローを使うこともあります。

ファイルフォーマットを策定し、一次解析ツールを作る人もいます。DNA-Seq では bwa、minimap2、samtools、bcftools、htslibなどがあります。ほかにも、アセンブラ hifiasm、リード分類器 Kraken2、トランスクリプト定量ツール Salmon、ゲノムグラフを扱う odgi などがあります。

なぜ SQLite ではなく BAM など業界専用の形式が採用されているのかという問いに対しては、単純な答えはありません。ゲノムデータは巨大であり、ゲノム座標に沿ったランダムアクセス、圧縮、インデックス、ストリーミング処理が必要になります。そのため、汎用データベースではなく、用途に合わせたファイル形式が発達してきたと考えられます。

配列の検索やアラインメントでは、次世代シーケンサーより前から使われている BLAST や MAFFT などもあります。教科書的には、Smith–Watermanアルゴリズムは局所アライメントを厳密に解く方法です。ただし実行時間が大きいため、大規模データではそのまま使いにくいです。そのため、ヒューリスティックな検索や、BWT / FM-index を利用する BWA のようなマッピングツールが使われます。

この領域では、C、C++、Rust などの静的プログラミング言語が使われます。データ構造、アルゴリズム、圧縮形式、インデックス、マルチスレッド、SIMD、メモリ管理、キャッシュ効率、I/O 最適化などの知識が必要になります。計算機科学やシステムプログラミングに近い領域です。

BAM/BCF/GTF/GFF を可視化するときは、IGVやJBrowseなどの GUI ゲノムブラウザを使います。IGV は Java Swing で作られており、JBrowse は TypeScript で作られています。前者はローカルでの閲覧に向き、後者はWebベースの表示や拡張性に強みがあります。ほかに UGENE も、さまざまなファイルを表示できます。

変異解析では、GATK HaplotypeCaller や Mutect2 で VCF/BCF ファイルが作成されます。変異コールという言葉は日本語にするとやや馴染みにくいですが、「変異であると判定する」「その候補として出力する」という意味に近いと思います。コールされた変異は、ただちに真実を意味するのではなく、データ、モデル、しきい値、フィルタリングに基づく判定結果です。

Broad Institute の GATK 系ツールは Java で実装されており、htslib ではなく htsjdk が使われています。この点は、C/C++ 系のツールが多い一次解析の世界のなかでは、異なる系譜に見えます。

一次解析の出力は、そのまま結論になるわけではありません。VCF、BAM、発現カウント行列などは、次の解析の入力になります。そこから、ドライバー変異、変異シグネチャー、TMB、コピー数変化、臨床情報との関連などを調べます。RNA-Seq では edgeR や DESeq2 を用いて発現変動解析を行い、その後 GSEA やパスウェイ解析に進みます。clusterProfilerなどが使われます。

scRNA-seq では Seurat や Scanpy を用いて、正規化、次元削減、クラスタリング、細胞型注釈、trajectory 解析、細胞間相互作用解析（CellChat）などが行われます。フローサイトメトリーでは、細胞集団の割合やマーカー発現を群間で比較します。FlowJoのような専用ソフトも使われます。

免疫染色では、陽性細胞率、細胞密度、腫瘍領域と間質領域の差、空間分布を解析します。プロテオミクスでは、差次的タンパク質発現、リン酸化変化、経路解析が行われます。

二次解析の世界

一次解析では、測定装置から出てきたデータを FASTQ、BAM、CRAM、VCF、BCF、発現カウント行列などの標準的な形式へ変換します。二次解析では、それらの出力を使って、生物学的・臨床的な意味を持つ比較を行います。

ここで重要になるのは、データを文脈に戻すことです。サンプルID、治療内容、採取時点、病理診断、奏効・非奏効、予後、検体品質、バッチ情報などを、一次解析の出力と対応づけます。どの患者の、どの検体から、どの条件で得られたデータなのかを確認しながら解析します。

がんゲノム解析では、VCF や MAF からドライバー変異、TMB、変異シグネチャー、コピー数変化、LOH、クローン構造、治療反応との関連などを調べます。RNA-Seq では、発現カウント行列から正規化、群間比較、発現変動解析、GSEA、パスウェイ解析へ進みます。scRNA-seq では、細胞 × 遺伝子の行列から、正規化、次元削減、クラスタリング、細胞型注釈、trajectory 解析、細胞間相互作用解析などを行います。

この段階で、データは表や行列として扱われることが増えます。VCF から抽出した変異表、MAF、RNA-Seq の count matrix、scRNA-seq の細胞メタデータ、臨床情報テーブルなどは、DataFrame として扱いやすいです。Python では pandas や polars、R では data.frame、tidyverse の dplyr、Bioconductor がよく使われます。Jupyter Notebook、RStudio、R Markdown などの環境では、表を確認し、列を追加し、フィルタリングし、可視化しながら解析を進めることができます。

変異の MAFファイルでは、各行に一つの変異を置き、染色体、位置、参照アリル、変異アリル、遺伝子名、Variant Classification、VAF、depth、サンプルID、患者ID、ドライバー遺伝子かどうか、COSMIC や ClinVar での注釈などを列として追加していきます。

RNA-Seq では、遺伝子ごとの発現量、log fold change、p 値、補正 p 値、pathway annotation などが列として追加されます。scRNA-seq では、細胞ごとに cluster ID、細胞型、UMAP 座標、サンプルID、患者ID、QC 指標などが付与されます。

二次解析では、新しい列を作り、データを絞り込み、別の表と join し、可視化しながら解釈を更新していきます。tidy data の形に整え、しきい値を使ってデータを選別し、lookup や join のような操作を繰り返しながら、解析対象を研究上の問いに近づけていきます。

Microsoft Excel や LibreOffice Calc のような汎用表計算ソフトも、実務上は使われます。CSV や TSV を目で確認するためには、SmoothCSV のようなビューアが使われることもあります。

Julia も数値計算やシミュレーションで使われます。ただし、生命科学の二次解析では、現時点では Python と R の利用が中心です。Bioconductor の蓄積があるため、統計解析やオミクス解析の多くは R で行われます。一方で、機械学習、可視化、ファイル処理、パイプライン周辺では Python も使われます。

CLI ツールの世界では、データ構造、アルゴリズム、圧縮形式、インデックス、I/O、並列処理が前面に出ます。Notebook と DataFrame の世界では、データの意味づけが前面に出ます。どの行を残すのか。どの列を作るのか。どの群とどの群を比較するのか。外れ値をどう扱うのか。臨床情報とどう結びつけるのか。どの可視化が問いに対応しているのか。こうした判断が解析の中心になります。

ただし、DataFrame には限界もあります。すべてのデータ構造が表の表現に向いているわけではありません。木構造やグラフは、表や行列として表現できます。たとえば、エッジリストを DataFrame として持つことはできますし、隣接行列として表すこともできます。ただし、表として表現できることと、その表現が解析に適していることは同じではありません。

C や Rust などで実装された低レベルなライブラリでは、グラフは DataFrame としてではなく、整数化されたノードID、エッジ配列、隣接リスト、あるいは CSR/CSC のような疎行列表現として保持されることが多いです。探索、近傍取得、最短経路、連結成分、クラスタリングなどの操作では、表形式よりもグラフ用の表現のほうが効率的だからです。

空間トランスクリプトームや病理画像では、位置情報も重要になります。細胞系譜、腫瘍内不均一性、治療抵抗性、時系列変化などを扱う場合にも、単純な行と列だけでは不十分になります。データは、表だけでなく、画像、多次元配列、座標付きの点集合、細胞間グラフ、時間変化を含む構造として扱われます。

新しい発見を行うには、既存の集計だけでは見えにくい構造や関係を示す必要があります。しかし、それを自分で一からモデル化し、アルゴリズムとして実装するのは簡単ではありません。そのため実際には、目的に近い解析ツールを探し、手元のデータに適用できるかを検討することが多くなります。

二次解析の後半では、DataFrame を中心とした操作だけでは扱いにくい対象が増えます。グラフ、木、ネットワーク、潜在変数モデル、確率モデル、系統樹、空間モデルなどが必要になります。これらは、データを集計するためだけでなく、観測されていない構造を推定するために使われます。

たとえば、がんのクローン進化解析では、観測できるのは変異アリル頻度、コピー数、サンプル間の変異の有無、細胞集団の割合などです。一方で知りたいのは、どの変異が先に起きたのか、どの細胞集団がどの細胞集団から分岐したのか、治療後に残った細胞がどの系統に属するのか、といった直接観測できない構造です。

このように、問題を設定して、観測値から隠れた構造を推定する必要が出てきます。

そのため、変分推論、MCMC、EMアルゴリズム、混合整数計画、グラフ探索、系統樹推定など推定や探索のための手法が使われます。ゲノム空間は広く、ノイズも多く、観測できる情報は限られています。そのため、常に厳密な最適解を求めるというより、仮定を置き、近似し、複数の可能性を比較しながら、妥当な説明を探します。

一方で、常に複雑な数理モデルが必要になるわけではありません。医学研究では、生存解析
、治療群間の有効性・安全性の比較、回帰モデル、傾向スコアを用いた観察研究の調整など、標準的な統計解析として扱う問題も多くあります。また、予測や分類を目的として、ランダムフォレストなどの機械学習モデルが使われる場合もあります。

おわりに

生命科学のデータ解析には複数の層があります。検体を集める層、測定する層、ファイルを作る層、表として整理する層、統計的に比較する層、観測されていない構造を推定する層です。

それぞれの層で必要な知識は異なります。実験、臨床、病理、シーケンス技術、Linux、ファイル形式、統計、可視化、アルゴリズム、数理モデル…。そのため、生命科学データ解析をひとくくりに捉えることは難しいです。

研究を遂行するためには、異なる背景を持つ人が必要です。実験や臨床を知る人には、データがどこから来たのかを理解する強みがあります。統計や数理を知る人には、観測値から構造を推定する強みがあります。プログラミングやシステムを知る人には、大量のデータを処理し、再現可能な道具を作る強みがあります。一方で、さまざまな制約から、協力関係を構築することは容易ではないという現実もあります。

生命科学の問いを計算機科学や統計学の人が扱える形に翻訳すること。逆に、計算機科学や統計学の道具を生命科学の現場で使える形に翻訳すること。その往復が、もう少し増えればよいと思います。

この記事では非常に限られた視点から、検体や実験から、一次解析、二次解析、モデル化に近い解析までを整理しました。

この記事は以上です