コロナウイルス変異株分析をしたイギリス発の論文の分析手法について

概要

注意：この記事は、2021-01-24に初稿を起こしていますが、急激に変わりゆく現状をもはや反映できておらず、最新の状況から乖離があることに留意してください。過去を振り返る参考程度に参照してください。この記事のタイトルにおける「変異株」とはB.1.1.7のことを指しています。

イギリスで現在も流行中のコロナウイルス変異株は、いよいよ日本にも拡大が始まる兆候が見えつつある今、このウイルスに関する定量的な分析を実施した論文を紹介します。私はよくあるデータ分析屋さん一人で、医療や疫学などとは全く接点がありませんが、この論文をお正月に読んで、不謹慎にも分析の鮮やかさに感動しました。この記事では、データ分析手法を考えるという観点から、イギリスでの変異株の流行をどのように把握したのかについて記載されている分析手法の一端を説明することを試みます。なお、当方は医学・疫学・分子生物系の知識に疎いので、その観点で間違い等あればご指摘いただければ幸いです。また、この数週間でも状況が大きく変わっているので、これを書いた時点(2021-01-24)での情報に基づく点をご理解ください。

対象論文¹ は大晦日の日付で、Imperial Collage London のサイトに掲載されました（ただし未査読論文であることに注意）。年末年始返上でこの文献を上梓された研究者の皆さんに頭が下がります。

Voltz+(2020-12-31) Report 42 - Transmission of SARS-CoV-2 Lineage B.1.1.7 in England: insights from linking epidemiological and genetic data

問題の背景

まず、イギリスで流行中のコロナウイルス"変異株"について説明します。このウイルスはB.1.1.7あるいは501Y.V1、VOC-202012/01などと呼ばれます(以下、論文の記載に準じ、B.1.1.7で統一します)。通常、コロナウイルスのようなRNAウイルスの突然変異率は高い傾向にありますが、この変異株B.1.1.7は23箇所の変異があり、特にスパイク状の形をしたタンパク質の部位に多数の変異があるのが特徴です²。このうちスパイクの501番のアミノ酸がNからYになることから、その変異を指してN501Yのように呼ばれることがあります。同様にN501Yの変異を持つ変異株として、南アフリカ型の501Y.V2、ブラジルの501Y.V3などが見つかっています。特にこのN501Yの変異があることにより、ヒトの細胞でスパイクタンパク質が結合するACE2と呼ばれる受容体への親和性を高めることにより、感染力が高まるのではと言われています（が、まだ確定的な証拠がありません）

このウイルスB.1.1.7の拡大が最初に顕著になったのはイギリスですが、イギリスにおいてこの事実が大きく取り上げられたのは、昨年2020年12月19日、イギリスのジョンソン首相によるロックダウンの再発動の決断です³。ジョンソン首相は、専門家を伴って会見を実施し「従来よりも最大70％感染しやすいようだ」という発言をし、その事実をもとにロックダウンに踏み切る決断をしました。その後のイギリスでの変異株拡大は周知のとおりです。

この発表は定量的な裏付けがよくわからないというのが私の印象でしたが、その後続々と「感染力が強い」ことの裏付けとなる詳細な報告が出てきました。中でも12/31にImperial Collage London から出た報告¹は、非常に短期間で出てきた内容としては、データ分析の観点から鮮やかであるという印象を受けました。以下ではこの論文の内容の一部を紹介します。

イギリスの変異ウイルスの検出経緯

（新型）コロナウイルスは、RNAと呼ばれる塩基が多数連なった構造をしています。この塩基の配列に準じてタンパク質が生成されることでウイルスが機能しますが、この塩基配列がコピーされる際にコピーが失敗するなどして、一部の塩基が置き換わったりなくなったりすることにより、タンパクのアミノ酸配列にも変異が生じます。これが変異株です。ウイルスの変異を正確にガチで見つけるには、本来塩基配列をすべて読み解く解析（以下、ゲノム解析）を実施する必要がありますが、現在の技術をもってしてもゲノム解析は時間がかかり、かつ高コストで、簡易検査における陽性検体のすべての配列を解析する体制を組めている国は現在ないといってよく、イギリスにおいても解析割合は全体の5-10%にとどまるようです²。

ゲノム解析以外に、PCR検査で新型コロナウイルスの存在を突き止める方法があります。これは、新型コロナウイルスRNAの特徴的な塩基配列をもつ部位を数か所特定し、RNAにその部位がすべて存在することをもって当該RNAが新型コロナウイルスであることを確定させます。この特定部位を見つけるのにつかわれるセンサーの役割を果たすものがPCRのプライマーと呼ばれます⁴。PCR検査はゲノム解析と比較すると低コストで所要時間も短くて済むという利点があります。ただし、PCR検査で変異株を見つけるには、変異した部位を検出できるようにプライマーが設計されている必要があり、そうでない場合にはコロナウイルスの存在は分かっても変異株であるかどうか一般にはわかりません。

ここでイギリスにとって幸か不幸か、偶然がありました。イギリスでコロナウイルスの検出用に使っているPCRのプライマーのうち1つが、対象となる変異ウイルスB.1.1.7の1つの変異のうち、スパイクタンパク質の69-70番のアミノ酸に対応する部分の欠失(69-70 deletionなどと文献に記載される)に関与していました。そのため、変異ウイルスB.1.1.7は新型コロナウイルスであるにもかかわらず、イギリスのPCR検査の手順でこの部位を正しく検出できず、陰性（偽陰性）となっていたのです。本来機能するはずのプローブが機能しないということで、この偽陰性が発生する状況は Spike Gene Target Failure 略して SGTF と呼ばれます⁵。

本来であればこの現象はプライマー設計の失敗で終わるのでしょうが、このSGTFの発生数が変異ウイルスの増加を顕著にとらえていました。12月に出された報告書に掲載された以下のグラフ⁶は、PCR検査における12月初旬までのSGTFの測定割合(黒実線, S dropoutと示されている時系列)と、同じ時期にサンプリングされた検体に対して実施されたゲノム解析によるB.1.1.7の検出割合(赤点線、B1.1.7と示されている時系列)です。両者は期間の末ごろにかけて同じ増加の傾向を示しています。PCR検査はゲノム解析と比較してサンプル数が多いことに加え迅速性にも優れるため、グラフにおけるSGTF検出割合はその先の期間で増加している傾向を示しており、ゲノム解析をしなくてもPCR検査によるSGTF発生数を先行指標として使うことで、簡易かつ迅速に流行を判定できることを示唆していました。

ただし、ここで注意が必要なのは、過去期間においてB.1.1.7が検出されていない時期にもSGTFが発生していた点です。これはSGTFの発生の原因となる69-70 deletion変異が、B.1.1.7ではないウイルスにも存在し、
SGTFが日常的に発生していたということです。 実際に、ゲノム解析によって69-70deletionを持つ他の変異ウイルス(B.1.1.7でないもの)を検出した割合(青点線、Other del6970と書かれているもの)が増えた時点(グラフのAug-Sepの間の小さなピーク)で、SGTFも増えてしまっています。つまりSGTFをそのままB.1.1.7の発生数と考えてしまうと、流行初期時点では誤検知を含んでしまうという問題があります。

データ分析による対象変異株の数の推定

おさらいすると、変異ウイルスとSGTFについては次の関係になっていました。

• 変異ウイルスB.1.1.7は、ゲノム解析によりほぼ確実に判定できるが、全量を判定するのはとても時間がかかり現実的ではない。
• 一方、このウイルスはイギリスにおける標準的な手順に従ってPCR検査を実施すると確実にSGTFを起こす。
• ところが、都合の悪いことに、同じようにSGTFを起こす別のウイルスもいるため、変異ウイルスB.1.1.7流行初期ではSGTFを呈する別ウイルスと誤検知をする可能性がある。

なので、検査の機動性の観点からは、PCR検査で発生するSGTFの数をカウントし、この数が変異株の数であると断定したいところなのですが、誤検知という厄介な問題があります。にもかかわらず大量に集まった過去のPCR検査の結果（そこにはSGTFの結果も含まれているもの）をすぐに捨てずに、なんとか変異株B.1.1.7の流行をとらえるマーカーとして使いたい！と研究者たちがあきらめなかったところが偉いところかもしれません。

彼らは次のように2段階で考えることにしました。

1. ゲノム解析のデータから、地域/日ごとに測定されたSGTFの検出が本当に（今確認したい）B.1.1.7をどれだけ含むかの検知率(True Positive Rate ; TPR)を算出する回帰モデルを立てる。つまりゲノム解析において、69-70deletionが発生しているケースで、どれだけの割合が変異ウイルスB.1.1.7となっているかの割合についての実績データから、真のTPRを算出する回帰モデルを作る。
2. このTPRを、注目している地域/日におけるSGTFの測定数に対して掛け算することにより、その地域/日における変異ウイルスB.1.1.7の検出数であると推定する。

ここで、1.についてもう少し補足します。ゲノム解析は全量の解析ではなくサンプル調査であるものの、ゲノム解析によるB.1.1.7の確定数を、ゲノム解析調査で69-70deletionが発生しているサンプル数で割り算すると、その時点でのSGTFの変異ウイルスマーカーとしてのTPRを簡易に求めることができます。ただしこの値はサンプル数も小さく揺らぎが大きい傾向にあります。

文献においては、次の2つの工夫により、測定誤差を吸収する工夫をしてます。

• TPRは時間的に緩やかに上昇することを考慮し、一般化加法モデルを使って、
  地域ごとに経過日数に対する3次のスプライン関数をフィッティングする。
• ただしTPRは隣接する地域間で連動する変動揺らぎがあることをモデルに取り込む⁷

モデリングの詳細はこれ以上は述べられていませんが、TPRの算出をただの割り算とせず、統計モデルの形に落とし込んで精緻に推定するステップを踏んだという点が重要です。実際にこの手順で算出したTPRを各週ごとのマップの形で図示した結果(対象論文のFigure S1からの抜粋)を以下に示します。流行初期では、TPRはほぼ全国的に0に近い値であったところ、流行が拡大するにつれて1の値に近づきます。

実際の推定ではこの値にSGTFの観測数をかけることで変異ウイルスの推定数とします。たとえばTPRが0の場合、SGTFの観測数があったとしても、変異ウイルスB.1.1.7は「1つも含まれていなかった」として推定することになります。このようにして得られた変異ウイルス数の推定結果の例を示します。次の図(対象論文の Figure 1-C,1-D) は、3地域におけるB.1.1.7の検出数・推定数を対数オッズにしたもので、左がゲノム解析の検出数、右がSGTFに対してTPRの算出地を掛け算して推定した値です。SGTFから推定した値の方が粒度が細かく、しかもゲノム解析よりも直近日での増加傾向をつかんでおり、流行初期における伝播状況をとらえやすいことが分かります。

これらの結果を用いて、従来の新型コロナウイルスに対する実効再生産数 $R_t$ の増加量を推定したところ（下図、対象論文のTable2）、ゲノム解析のサンプルに対して決定された値(Data of Variants=Genomic)では、中央値で0.48～0.68となり、各手法の信頼区間もかなり広かったところ、SGTFにTPRを掛けて推定した値(Data of Variants=TPR-adjusted SGTF)を用いると、$R_t$ の増加は（モデル種別や対象地域区分で違いはあるものの）信頼区間はやや狭くなって推定値の確実性が増し、中央値ベースではゲノム解析の値よりも小さいものの、依然として0よりは十分に高い0.36-0.48程度という推定値を得ています。現在「イギリスの変異ウイルスの実効再生産数は0.4程度上がる」とされている値は、TPRで補正されたSGTFから推定された結果を用いて求められた値を根拠にしているところが大きいと思われます。

考察

これらの分析には重要な示唆があります。

本来、変異ウイルスの検出に必要なゲノム解析は、全量検査をおこなうには高コストかつ時間がかかるものです。イギリスででもゲノム解析が行われていたものの、サンプル調査であり網羅性・速報性が低いという問題点がありました。一方で、偶然であるとはいえ、PCR検査の副産物としてのSGTFがイギリスでの変異株の増加を（不正確さが付きまとうものの）捉えられてるということを積極的に利用し、この値をサンプル調査から求めたTPRによって補正することにより時系列的・地域的な流行拡大をとらえることができました。当該研究者は数として勝るSGTFの結果と、少数のサンプルに対するゲノム解析結果を組み合わせることで、ある種の厳密さを捨てつつも、PCR法による測定結果の速報性と数量の豊富さを重視することで、むしろ感染の性質に関してより詳細な結果を得ることができています。

イギリスにおける流行初期で、SGTFを積極的に生かし、かつこれらの解析を短期で実施することで感染拡大の評価に用いて政治判断にまで結び付ける材料の一翼を担ったという点は、非常に注目すべき点です。これは、医学的見地に詳しい専門家に加え、統計解析の専門家たちが自分たちの知見を加えてより具体的な結果を出した好事例と私は見ます。この分析によって国境に対するポリシーを素早く変えて、世界が変異ウイルスに対する準備期間を得ることができたと言えるでしょう。

私は医療や疫学の素人なので、それらの観点でのこれ以上の考察は憚られますが、医療に限らず様々な領域（例えばマーケティングなど）の領域で、非常に高コストな調査・検査が要求されるシチュエーションにおいても、そのような調査などをせずとも、より機動性の高い代替調査との関連性を見ておくことで、高コストな調査の代替/先行指標/マーカーとして補完しながら利用することで、むしろ精緻な結果を出すことができることもあり得るという事例としてみることもできます（データフュージョンと言われたりします）。ただしこのような方法は適用を間違うとノイズを多く拾ってしまう場合もあって、慎重に行うべき分析でもあります。対象論文は未査読ではあるものの、ノイズを避けつつより詳細な結果を出すことに成功した一事例として、個人的には非常に鮮やかな結果だなと感じました。

1. Voltz+(2020-12-31) Report 42 - Transmission of SARS-CoV-2 Lineage B.1.1.7 in England: insights from linking epidemiological and genetic data ↩ ↩²

2. 国立感染症研究所(2021-01-02) 感染性の増加が懸念されるSARS-CoV-2新規変異株について（第4報） ... 以下2021-05-12追記。現時点でのイギリスにおけるゲノム解析割合は、50%-60%に達するというレポートがありますので「5%-10%に留まります」とした記述は現時点での実態とはもはや乖離しています。 ↩ ↩²

3. NHKニュース(2020-12-20) コロナ ロンドンなど外出制限へ “変異ウイルスで感染急拡大” ↩

4. 例えば次を参照。 PCRとRT-PCRの基本原理 ↩

5. 本稿執筆時点で日本において感染研が推奨しているプロトコルでは、PCR検査のプライマーはB.1.1.7の変異領域が対象となっておらず「これまでと同様に使⽤可能である」とされています。そのため、同プロトコルに準じたPCR検査を主にしている日本ではSGTFは発生しないことになります。ただし、日本において現在23箇所の変異の1つであるN501Y変異を直接ターゲット部位とするPCR検査を開発中であるという報告もあり、進捗が期待されます。（例えば東京都健康安全研究センターの報告：[変異株を追え！～新たな迅速判定法の開発～] (https://note.com/tokyo_icdc/n/n82950294c56a) などを参照。2021-01-24閲覧） ↩

6. この図は以下からの引用。Investigation of novel SARS-COV-2 variant / Variant of Concern 202012/01 ↩

7. このあたり、論文においてはGaussian Markov Random Field を適用したとなっています。具体的なモデリングの定式化やソースコードが見当たらないので詳細は不明ですが、おそらくスプラインからの誤差項について、隣接する地域間で正の共分散を持つようにモデリングしたものと思われます。 ↩