シングルセル遺伝子発現解析のための簡便なシステムが登場してすでにもう久しく、かなり普及したといってよいのではないでしょうか。しかし、簡便になった分、システムの「クセ」に気を留めずに始める方も少なくないような気がします。ラフですが、「クセ」とその対応の仕方について下に簡単に記しました。当研究室が分担機関として進めているメダカNBRPの事業としての活動をまとめたものですが、他生物での研究にも活かしていただけると考えています。
（おことわり：下記、不備があり、何度か手直しをしております）
メダカのシングルセル遺伝子発現解析のための遺伝子モデルの改変
森上賢治、工樂樹洋（ともに分子生命史研究室）
シングルセル遺伝子発現解析（下図1に例）では、往々にして、タンパク質コード領域のさらに下流、すなわち3’非翻訳領域（UTR）の転写産物配列が得られる。こういったリードをリファレンスゲノム配列にマッピングして遺伝子発現量を定量するのだが、その際にタンパク質コード領域だけを考慮した遺伝子モデル（GFF/GTFファイルとして配布されている）を使用した場合には、得られる肝心の転写産物配列が往々にして3’UTRに位置するために、発現量のカウントに加味されにくい。NCBIやEnsembなどの公共のデータベースで配布されている遺伝子モデルにもUTRは含まれているが、不十分であることが多く、これがネックになりやすい（Univ Cambridgeのコース“Analysis of single cell RNA-seq data”の資料のセクション3.5を参照）。
図1 典型的なシングルセル遺伝子発現解析のアウトプット
この3'UTRへのリードの偏り（下図2の中段）が原因で、シングルセル遺伝子発現定量の結果に、興味ある遺伝子の発現量が正しく反映されていない、であるとか、全体のプロファイルが予想と大きく異なる、という声をよく耳にする。シングルセル遺伝子発現解析システムを販売している業者のほうでも、(cellranger mkrefのガイドにあるように、情報整備のままならない生物のための最低限の手立てが10x Genomics社によって用意されているものの)こういった問題についての周知が十分にはなされておらず、対策どころか問題自体を整理して文字にした記事も希少である。貴重な試料を用いてサンプルを調製し、シークエンスデータを得て解析を進めて初めて問題に気づいたという研究者も少なくないようだ。本記事では、どの生物でも研究者が直面しうるこの問題の周辺を事前に把握したうえで取り組まれる方が増えるよう、我々のメダカでの経験に基づいて情報提供を試みる。なお、この試行は、少なくとも、神戸の理研CLSTおよび理研BDRの分子配列比較解析チームにて、技師の西村理博士がすでに手掛けており、その一部の努力は遺伝子モデルのさほど整っていない哺乳類のシングルセル解析(Bilgic et al., eLife 2023)などに活かされている。
図2 メダカゲノム内のある遺伝子領域（3'末端）の例
mRNAの終点がどのように決まるのか、意外と知られていない。ヒトの場合には塩基配列AAUAAAが主要なポリアデニル化シグナル（polyadenylation signal, PAS; Proudfoot et al., Gen Dev 2011）とされ、ほかにAUUAAAなども使われるらしい。この配列モチーフが認識された後にRNAが切断されいわゆるポリAが付加されるが、その位置は一定ではない。開始コドンやイントロンの始点・終点ほどは厳密に同定しにくいのである。mRNAの始点（転写開始点、すなわちtranscription start site, TSS）の決定にはその目的に特化したCAGE法が用いられるように、mRNAの終点の決定を目標にしたデータ取得法として、Quant-seq（Lexogen社がキットを販売）などがある。余談だが、Quant-seqキットを試用したことのある工樂は、2019年のPlant Animal Genome会議にて、Lexogen社主催の企業ワークショップで講演した際にQuant-seqキットを用いて取得したデータを3’UTRの終点の認識に利用し、10X Chromiumデータの徹底利用に繋げるアイデアを披露した。その後、この話題からは興味が遠のいたものの、2023年度から加わったメダカのNBRPの一員となったことをきっかけに、再び視野に入ってくることとなった。
コーディング領域偏重の遺伝子モデルでは不十分との考えによる改善の努力は、一部の生物では行われてきた。ゼブラフィッシュでは、とくに3’ UTRに着目したEnsemblとNCBIの比較に基づいて両者の情報を取り込んだ改善版遺伝子モデルが作出され(Lawson et al., eLife 2020)、公開されている（Lawson研究室のページ）。メダカにおいても同様の問題が存在し、一部の遺伝子にUTRを付加した遺伝子モデルを用意して解析を進めた例もあるが（Siddique et al., Sci Rep 2021）、包括的な対応はどうやらなされていない。筆者たちが属し、NBRPメダカプロジェクトにおいて生命情報の整備を担当する国立遺伝研分子生命史研究室では、メダカにおいてこの問題に対処すべく、同NBRPプロジェクトに参加する研究者と共に、既存の遺伝子モデルの改善を試みた。
我々がまず最初に注目したのは、stringtieというツールである。ショートリードアラインメント結果をインプットすると、タンパク質コード領域の推定に基づく遺伝子モデルにUTRを付与するなど、転写産物エビデンスを加味した遺伝子モデルを出力する。次に注目したのが、2023年に入ってから公開されたpeaks2utr(Haese-Hill et al., Bioinformatics, 2023)というツールである。まさに10X Chromiumシステムを用いたシングルセル遺伝子発現解析のためのツールという触れ込みで公表され、インプットデータとしては、同システムを用いて取得したデータを想定している。つまり、それ以外のデータの取得を前提とはしていないのである。いっぽうで、前者のstringtieを利用する際には、別途バルクRNA-seqデータなどの存在が前提となる。我々は、peaks2utrが最も有望と目しながらも、インストールにてこずったため並行してstringtieの試用を進めた。stringtieのほうでもある程度有望な改善がみられていたが、途中まで進めたところでpeaks2utrのインストールの問題（下記）が解決したため、途中からはpeaks2utrに絞って進めた。なお、peaks2utrの動作機序については、開発者が出版した論文に付属のSupplementary Textも参照されたい。
改善を目指す以上、遺伝子モデルの全体的な比較と評価のための手段と判断基準を用意しておくことが重要である。我々は、手製のスクリプトを作成するとともに、AGAT(Another GFF/GTF Analysis Toolkit)という既成プログラム群のうち agat_sp_statistics.pl や agat_compare_two_annotations.pl などを利用し、GTFファイル内の情報を確認することをルーティン化している。
いざ、peaks2turをインストールしようとすると大きな壁が立ちはだかった。ピークコールのための他のプログラムMACSの利用を前提としているが、着手当初（2023年6月頃）は、MACSのインストールに必要なCythonのバージョンアップのため、peaks2utrのインストールが完了できないという問題である（参考：MACS3のCython不適合を指摘したGithubのIssue）。MACS2は新バージョンのCythonに先に適合したため(v2.2.9.1)、MACS2を利用する手立ても考えた。幸い、その後すぐにMACS3も適合したため（v3.0.0b3）、現在では問題なくインストールできることを少なくとも我々のほうでは確認できている。
peaks2utrへインプットする転写産物の情報として、開発者は10X Chromiumデータ（例、メダカ下垂体雌雄scRNA-seqデータ）を想定しているが、我々はそれに加えて、いわゆるバルクRNA-seq情報（メダカ11組織 SRP044784）も使用した。後者を用いて、我々がメダカのEnsembl由来のGTFファイルに対してpeaks2utrを使用した結果、7175個の転写産物の既存のUTRが延伸されていた（図2の下段に一例）。UTRを延伸することで起きうる弊害が、（隣の遺伝子に付随する3’末端のリードの集積を拾ってしまうことによる）隣の遺伝子が占める領域との融合である。peaks2utrはこれが起きないように設計されているという（Haese-Hill et al., 2023）。
ちなみに、RNA-seqデータをもとにUTRを推定するためのツールとしては、ExUTRや遺伝子推定プロブラムAugustus/Brakerの作者が開発したGUSHRのほかに, GETUTR, UTRme, F3UTERなどがある。
周辺情報を日本語でまとめておこうという意図で本記事を書いてみたが、今回の一連の試行プロセスについてはこちらの英語ページで細かく説明しており、我々が改変したGTFも同サイトからリンクを辿れば入手可能である。自分が扱っている生物でUTR重視の遺伝子モデルを作出したいという方にはぜひ目を通していただき、活用の可能性をご検討いただければ幸いである。この記事の草稿を書いて、公開準備をしている間にヒトやマウスでの遺伝子モデルが不完全であることの影響に注目した論文が出版された（Pool et al., Nat Methods, 2023)。今後、さらにこの問題に気を留める必要が出てくるであろう。
メダカNBRP 生命情報チームによるオリジナルページ
3' UTRを重視した遺伝子モデルの改変ワークフロー https://github.com/Squalomix/utr-modeling/
MedakaBase https://medakabase.nbrp.jp/
Medaka Omics Reference https://github.com/Squalomix/medaka-annex
謝辞
国立遺伝学研究所の丹羽大樹博士、川口茜博士、佐藤茉菜博士、矢倉勝博士、原雄一郎博士には、オミクスデータ取得とインシリコ解析全般について、ナショナルバイオリソースプロジェクト（NBRP）の成瀬清博士および安齋賢博士にはメダカのオミクスデータ整備状況とその利用について、理化学研究所生命機能科学研究センターの門田満隆博士および西村理博士、国立遺伝学研究所の久保郁博士および松田光司博士にはゼブラフィッシュなどの生物の遺伝子モデルの整備状況や改善のための手法について、それぞれ貴重な議論をさせていただいたことに関してここに感謝を申し上げます。