はじめに
このPythonコードは、内閣官房下水サーベイランス実証事業の結果公開されたデータをもとに、プログラミングを使った簡単なデータ取得、データクリーニング、 QC/QA、データ解析の例を示す教育的目的で書かれている。本情報は、下水サーベイランスデータそのものを理解することが目的ではないため、利用者はデータの特性を詳細に理解した上で、データ分析の結果そのものは注意して取り扱うべきである。また、企業・アカデミア・個人問わず、活用できるコードがあればぜひ活用していただきたい。質問等あれば京都大学招聘研究員、遠藤礼子 (endo.noriko.3p@kyoto-u.ac.jp)まで。
データソース:https://corona.go.jp/surveillance/
Code developed by: Dr. Noriko Endo
Code last modified: 2023-05-15
パッケージ(ライブラリ)をあらかじめ読み込んでおけば、いちいち難しいコードを書かなくとも「ピアソン統計を計算して」、「見栄えの良いラインプロットを作って」、「日本語やアラビア語を使いたい」などが簡単にできる。よく使うパッケージは毎回デフォルトで読み込みようにすると良い。
# 必要なパッケージのダウンロード
import gspread as gs
import pandas as pd
import numpy as np
import scipy.stats
import datetime
import random
import math
import os, sys, glob
import io
import fnmatch
from dateutil.relativedelta import relativedelta
from scipy.optimize import curve_fit
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
fp = FontProperties(fname=r'/Library/Fonts/Osaka.ttf', size=12)
import japanize_matplotlib
import plotly.io as pio
pio.renderers.default = 'notebook'
!pip install Pyppeteer
!pyppeteer-install
import seaborn as sns
sns.set(font='IPAexGothic')
%matplotlib inline
Requirement already satisfied: Pyppeteer in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (1.0.2) Requirement already satisfied: urllib3<2.0.0,>=1.25.8 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (1.26.11) Requirement already satisfied: tqdm<5.0.0,>=4.42.1 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (4.64.1) Requirement already satisfied: appdirs<2.0.0,>=1.4.3 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (1.4.4) Requirement already satisfied: certifi>=2021 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (2022.9.24) Requirement already satisfied: websockets<11.0,>=10.0 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (10.4) Requirement already satisfied: pyee<9.0.0,>=8.1.0 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (8.2.2) Requirement already satisfied: importlib-metadata>=1.4 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from Pyppeteer) (4.11.3) Requirement already satisfied: zipp>=0.5 in /Users/norikoendo/opt/anaconda3/lib/python3.9/site-packages (from importlib-metadata>=1.4->Pyppeteer) (3.8.0) chromium is already installed.
データを読み込み、後にデータ解析で使いやすいようにデータの記述用法を統一したりデータの型を揃えるデータクリーニングを行う。
用いたデータは内閣官房ウェブサイト(https://corona.go.jp/surveillance/) からダウンロードし、ローカルのコンピュータに”../data/内閣官房実証事業_final/”というフォルダを作り、ダウンロードしたすべてのデータを格納した。
# ファイルの読み込み
f_name = '../data/内閣官房実証事業_final/プロファイル.csv'
df_meta = pd.read_csv(f_name,index_col=None, header=1)
# データフレーム(表みたいなもの)のカラムを設定
df_meta.columns = df_meta.columns.str[1:]
# データクリーニング:表示名を統一
df_meta.replace('沈殿物抽出法(仮称)','沈殿物抽出法',inplace=True)
df_meta.replace('EPISENS-S法(パッシブ)','EPISENS-S法',inplace=True)
# ファイルの読み込み
## ファイル名に’ウイルス’がついている複数のファイルを一気に読み込む
rootdir = '../data/内閣官房実証事業_final'
df=[]
for root, subdirs, files in os.walk(rootdir):
for f in files:
if 'ウイルス' in f:
# if '仙台市_検査結果(ウイルス).csv' in f: ## 全部一気にデータを読む前に、ファイル1つだけ読み込んでうまくいくかチェック
# read raw data
filename = root + '/'+f
# print(filename) # ファイル名のチェック
df_tmp = pd.read_csv(filename, index_col=None, header=0)
# Convert wide data to long/tidy data
## このへんは元のデータ形式がワイドデータだったことやカラムがたくさんあったため、ちょっと乱暴な力技
df_tmp.iloc[:,0].fillna(method='ffill',inplace=True)
df_tmp.iloc[:,0] = df_tmp.iloc[:,0]+ ';' +df_tmp.iloc[:,1]
df_tmp2 = df_tmp.transpose().fillna(method='ffill')
df_tmp3 = df_tmp2.set_index([0,1,2,3,4,5,6])
df_tmp3.columns = df_tmp3.iloc[0]
df_tmp3 = df_tmp3.iloc[2:]
df_tmp4 = df_tmp3.melt(var_name='datetime',value_name='results',ignore_index=False).reset_index()
df_tmp4.drop(columns=1,inplace=True)
# データフレーム(表みたいなもの)のカラムを設定
df_tmp4.columns = ['city','loc_num','loc','protocol_num','target','junk','datetime','results']
df_tmp4.drop(index=[0,1],inplace=True)
#宇部市はあらかじめデータフォーマットを少し変更した’宇部市_検査結果(ウイルス)_NEmod.csv’を用いる。
if not f == '宇部市_検査結果(ウイルス).csv':
df.append(df_tmp4)
# convert to dataframe
df_v = pd.concat(df, axis=0, ignore_index=True)
df_v[['date','time']] = df_v['datetime'].str.split(";", expand = True)
# データクリーニング
"""
データ解析をする際に、表記名が揃っていないと扱いにくい。
例えば'SARS-CoV-2(N1)'が'SARS-Cov-2(N1)'("v"が全角小文字);'SARS-Cov-2(N1)'(”(N1)”が全角);'CDC_N1';'CDC N1'
とバラバラであると別なものと認識され、「'SARS-CoV-2(N1)'での平均を計算して」などのオペレーションが簡単にはできない。
"""
df_v['target'].unique()
df_v.replace('SARS-Cov-2(N1)','SARS-CoV-2(N1)',inplace=True)
df_v.replace('SARS-Cov-2(N1)','SARS-CoV-2(N1)',inplace=True)
df_v.replace('SARS-CoV-2(N1)(Copies/l)','SARS-CoV-2(N1)',inplace=True)
df_v.replace('SARS-CoV-2 RNAウイルス(N1)(Copies/l)','SARS-CoV-2(N1)',inplace=True)
df_v.replace('CDC_N1','SARS-CoV-2(N1)',inplace=True)
df_v.replace('CDC N1','SARS-CoV-2(N1)',inplace=True)
df_v.replace('SARS-Cov-2(N2)','SARS-CoV-2(N2)',inplace=True)
df_v.replace('CDC_N2','SARS-CoV-2(N2)',inplace=True)
df_v.replace('PMMoV(Copies/l)','PMMoV',inplace=True)
df_v.replace('phi6','Phi6',inplace=True)
# the following 3 locations don't indicate which primer was used in the file:'小松市', '埼玉県', '札幌市'
# '小松市' --> N1 only
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='小松市')] =\
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='小松市')].replace('SARS-CoV-2','SARS-CoV-2(N1)')
# '埼玉県'-> N1 only
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='埼玉県')] =\
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='埼玉県')].replace('SARS-CoV-2','SARS-CoV-2(N1)')
# '札幌市'-> N1 only
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='札幌市')] =\
df_v[(df_v['target']=='SARS-CoV-2')&(df_v['city']=='札幌市')].replace('SARS-CoV-2','SARS-CoV-2(N1)')
# Check if the parameter names are standardized
df_v['target'].unique()
array(['SARS-CoV-2(CDC N1+N2)', 'PMMoV', 'Phi6', 'SARS-CoV-2(N1)', 'BOD5', 'SARS-CoV-2 RNAウイルス(N1)(Copies/lCotton)', 'PMMoV(Copies/lCotton)', 'SARS-CoV-2(N2)', 'PMMoV回収率', 'CDC_N1 [copies/sample]', 'CDC_N2 [copies/sample]'], dtype=object)
# remove Ct values -> '甲府市' has Ct values, instead of concentrations, in results
df_v[df_v['junk'].str.contains('Ct')]['city'].unique()
df_v = df_v[~df_v['junk'].str.contains('Ct')]
# データの型(object, integer, float, dateなど)を揃える
df_v.replace('1008-1','1',inplace=True)
df_v['loc_num'] = df_v['loc_num'].astype('int64')
df_v['protocol_num'] = df_v['protocol_num'].astype('int64')
df_v.loc[df_v['date'].str.contains('~',na=False),'date'] =\
df_v[df_v['date'].str.contains('~',na=False)]['date'].str.split('~',expand=True)[0]
df_v.loc[df_v['date'].str.contains(' ',na=False),'date'] =\
df_v[df_v['date'].str.contains(' ',na=False)]['date'].str.split(' ',expand=True)[0]
df_v['date'] = df_v['date'].astype('datetime64[ns]')
"""
数量データに'定量下限値未満'や'ND'というテキストデータがあるとデータの型が揃えられないので、
np.nan(数量データとして扱われるNaN)に置き換える
"""
df_v = df_v[~df_v['results'].str.contains('ー')] # I don't know what 'ー' means...
df_v['results'].replace('定量下限値未満',np.nan,inplace=True)
df_v['results'].replace('定量下限未満',np.nan,inplace=True)
df_v['results'].replace('ND',np.nan,inplace=True)
df_v['results'].replace('N.D.',np.nan,inplace=True)
df_v['results'].replace('N.D',np.nan,inplace=True)
df_v['results'].replace('UD',np.nan,inplace=True)
df_v['results'].replace('BLQ',np.nan,inplace=True)
df_v = df_v[~df_v['junk'].str.contains('採水時刻')]
df_v = df_v[~df_v['junk'].str.contains('採水開始時刻')]
df_v = df_v[~df_v['junk'].str.contains('採水開始日時')]
df_v = df_v[~df_v['results'].str.contains('\*',na=False)] # due to inhibition
df_v['results'] = df_v['results'].astype('float')
# 他に(効率的にデータ解析をするために)除くべきデータがないかチェック
df_v.groupby(by=['target','junk']).count()
df_v = df_v[~df_v['junk'].str.contains('回収率(%)')]
df_v = df_v[~df_v['junk'].str.contains('定量')]
# ファイルの読み込み
## ファイル名に’外部環境’がついている複数のファイルを一気に読み込む
rootdir = '../data/内閣官房実証事業_final'
df=[]
for root, subdirs, files in os.walk(rootdir):
# print(subdirs)
for f in files:
if '外部環境' in f:
# read raw data
filename = root + '/'+f
# print(filename)
df_tmp = pd.read_csv(filename, index_col=None, header=0)
df_tmp.drop(index=[1,6],inplace=True)
df_tmp.iloc[:,0] = df_tmp.iloc[:,0]+ ';' +df_tmp.iloc[:,1]
df_tmp2 = df_tmp.transpose().fillna(method='ffill')
df_tmp3 = df_tmp2.set_index([0,2,3,4,5])
df_tmp3.columns = df_tmp3.iloc[0]
df_tmp3 = df_tmp3.iloc[2:]
df_tmp4 = df_tmp3.melt(var_name='datetime',value_name='results' ,ignore_index=False).reset_index()
df_tmp4.columns = ['city','loc_num','loc','env_var','unit','datetime','results']
df.append(df_tmp4)
#convert to dataframe
df_env = pd.concat(df, axis=0, ignore_index=True)
df_env[['date','time']] = df_env['datetime'].str.split(";", expand = True)
# データクリーニングをするために、どんなデータがどんなユニットで格納されているかチェック
## …ちょっと難しそうなので、データクリーニングは諦め、今回環境データは扱わないことにする
df_env['unit'].unique()
array(['℃', '-', '度', 'm3/h', 'CFU/mL', 'm3/日', 'mm', 'm^3', 'mg/L', 'm3', nan], dtype=object)
# ファイルの読み込み
"""
ファイル名に’疫学データ’がついている複数のファイルを一気に読み込む。
レポーティングの単位(市ごと、処理区行くごと、など)が違ったり、レポーティングの方法が途中で変わり統一的なデータ処理が難しかったため、
提供されたすべてのデータを使うことを諦め、各事業主体ごとに最もデータ数の多いものを疫学データの代表値として扱った。
単位の違うデータをそのまま実数値で扱うことは難しいため、この後の疫学データを用いたデータ解析はすべて増加率を用いる。
"""
rootdir = '../data/内閣官房実証事業_final'
df=[]
for root, subdirs, files in os.walk(rootdir):
# print(subdirs)
for f in files:
if '疫学データ' in f:
# if '沖縄県_疫学データ.csv' in f:
# read raw data
filename = root + '/'+f
# print(filename)
df_tmp = pd.read_csv(filename, index_col=None, header=0)
df_tmp.replace('ー',np.nan,inplace=True)
df_tmp.iloc[0,2:] = df_tmp.iloc[0,1]
df_tmp2 = df_tmp.iloc[:,1:].count()
usecol = pd.DataFrame(df_tmp2).idxmax(axis=0)[0]
df_tmp.drop(index=[1,5,6,7],inplace=True)
df_tmp3 = df_tmp[['管理番号',usecol]].transpose()
df_tmp3.columns = df_tmp3.iloc[0]
df_tmp3 = df_tmp3.iloc[1:]
df_tmp4 = df_tmp3.set_index(['実証主体名', '処理場名', '地域の単位', 'データソース'])\
.melt(var_name='date',value_name='cases' ,ignore_index=False).reset_index()
df_tmp4['date'] = df_tmp4['date'].astype('datetime64[ns]')
df_tmp4 = df_tmp4[~df_tmp4['cases'].str.contains('\*',na=False)]
df_tmp4['cases'] = df_tmp4['cases'].astype('float')
df_tmp4['cases_7d_ave'] = df_tmp4['cases'].rolling(window=7).mean()
df.append(df_tmp4)
#convert to dataframe
df_epi2 = pd.concat(df, axis=0, ignore_index=True)
df_meta.replace('高知県','高知県・高知市',inplace=True)
df_v.replace('高知市','高知県・高知市',inplace=True)
df_meta.replace('埼玉県・さいたま市・秩父市・日高市・飯能市','埼玉県+4市',inplace=True)
df_v.replace('埼玉県、さいたま市、秩父市、飯能市、日高市','埼玉県+4市',inplace=True)
df = df_v.merge(df_meta,left_on=['city','loc','protocol_num'],right_on=['共同体名','採水地点名','プロトコル名'],how='left')
調査結果を見る前に、まずはこの調査にどのようなデータが含まれているのかみてみる。用いた検査手法、採水方法、採水地点の種類に対し、どのくらいの自治体数とサンプリング箇所があったか計算した。データクリーニングさえできていれば、こういった解析がコード1行でできる。
df.groupby(by=['検査手法名','採水方法','採水地点種別'])[['city','loc']].nunique().rename(
columns={'city':'自治体数','loc':'サンプリング箇所数'})
自治体数 | サンプリング箇所数 | |||
---|---|---|---|---|
検査手法名 | 採水方法 | 採水地点種別 | ||
Direct Capture法 | グラブ | 処理場 | 2 | 9 |
EPISENS-S法 | グラブ | その他(返流水) | 1 | 2 |
ポンプ場 | 3 | 6 | ||
処理場 | 10 | 27 | ||
コンポジット | ポンプ場 | 2 | 3 | |
処理場 | 6 | 20 | ||
パッシブ | マンホール | 2 | 4 | |
PEG沈殿法 | グラブ | 処理場 | 2 | 8 |
コンポジット | ポンプ場 | 1 | 2 | |
処理場 | 2 | 5 | ||
沈殿物抽出法 | グラブ | マンホール | 1 | 2 |
処理場 | 1 | 4 | ||
パッシブ | マンホール | 1 | 1 | |
脱脂綿-Direct capture法(C+DC法) | パッシブ | 処理場 | 1 | 4 |
遠心沈殿法 | グラブ | 処理場 | 1 | 1 |
データを解釈しようとする前に、まずはデータのバラツキ・信頼区間を確かめなければならない。
バラツキのソースとしては以下のようなものがある。
実証事業ではRNA抽出液を何回かPCRにかけている自治体があるため、1によるバラツキが評価できる。
まずは、どの自治体、検査手法が各ターゲットに対し何回(何ウェル)平均でPCRをかけてたのかみてみた。
仙台市、大分市、宇部市、小松市、甲府市、高知県・高知市以外、大多数が1ウェルだけでの調査を行なっていたようである。
注1:データシートに載っていないだけで、複数のウェルでの調査を行なっていた可能性はある。\ 注2:データクリーニングの過程でCt値の情報を除いたため、ここで表示される甲府市の結果は正しい統計を反映していない。元データによると、SARS-CoV-2はN1N2同時検知を4ウェルで、PMMoV及びPhi6はそれぞれ2ウェルで行なっていたようである。
# 手法・自治体ごとに平均何ウェルでの測定をしていたか?
df_tmp = df[df['city']!='甲府市'] # Because Kohu uses Ct values instead of concentration in the data sheet
df_rep_per_sample = df_tmp.fillna(-1).groupby(by=['city','loc','datetime','target','プロトコル名','検査手法名']).count()\
[['results','junk']].reset_index()
df_ave_rep_per_city = df_rep_per_sample.groupby(by=['検査手法名','city']).mean()[['results']]
df_ave_rep_per_city.rename(columns={'results':'Number of PCR replicates'})
Number of PCR replicates | ||
---|---|---|
検査手法名 | city | |
Direct Capture法 | 宇部市 | 3.000000 |
滋賀県 | 2.986111 | |
EPISENS-S法 | 久留米市 | 1.000000 |
京都府 | 1.000000 | |
仙台市 | 2.991667 | |
四日市市 | 1.000000 | |
埼玉県 | 1.000000 | |
埼玉県+4市 | 1.000000 | |
札幌市 | 1.000000 | |
水戸市 | 1.000000 | |
沖縄県 | 1.000000 | |
浜松市 | 1.418719 | |
神奈川県 | 1.000000 | |
秋田県 | 1.000000 | |
養父市 | 1.000000 | |
香川県 | 1.000000 | |
PEG沈殿法 | 大分市 | 2.495690 |
宇部市 | 3.000000 | |
小松市 | 2.052373 | |
水戸市 | 1.000000 | |
沈殿物抽出法 | 高知県・高知市 | 2.996435 |
脱脂綿-Direct capture法(C+DC法) | 宇部市 | 3.000000 |
では、それぞれの手法及びターゲットに対し、PCRにおけるばらつきはどれくらいあるのだろうか?変動係数(Coefficient of variation = standard deviation / mean)をサンプルごとに計算し、その自治体ごとの平均、そしてその平均を計算した。
PCR測定における変動係数は0.2-0.4程度。つまり、もし1ウェルでしか測定をしていない場合、2-4割程度の濃度のばらつきは信頼区間の範囲内と捉えることができる。1ウェルでなく、複数ウェルで測定をする場合、この信頼区間を縮めることができる。
また、バラツキはPCR測定以外からも生じることに注意する。つまり、下水サーベイランスにおける1つのデータポイントの信頼区間はもっと大きい。これも、測定を何度か繰り返したり、いくつかのデータポイントの平均(時間平均、空間平均)を取ることで信頼区間をもっと縮めることができる。
df_cv_per_sample = df.groupby(by=['city','loc','datetime','target','プロトコル名','検査手法名']).std()\
/df.groupby(by=['city','loc','datetime','target','プロトコル名']).mean()
df_ave_cv_per_city = df_cv_per_sample.groupby(by=['city','検査手法名','target']).mean()[['results']]
df_ave_cv_per_city = df_ave_cv_per_city[~df_ave_cv_per_city['results'].isna()].rename(
columns={'results':'Coefficient of variation in replicates'})
df_ave_cv_per_city.reset_index().groupby(by=['検査手法名','target']).mean()
Coefficient of variation in replicates | ||
---|---|---|
検査手法名 | target | |
Direct Capture法 | Phi6 | 0.195461 |
SARS-CoV-2(N1) | 0.266040 | |
SARS-CoV-2(N2) | 0.232406 | |
EPISENS-S法 | PMMoV | 0.057236 |
SARS-CoV-2(N1) | 0.085380 | |
PEG沈殿法 | PMMoV | 0.144703 |
Phi6 | 0.141204 | |
SARS-CoV-2(N1) | 0.396853 | |
SARS-CoV-2(N2) | 0.379619 | |
沈殿物抽出法 | SARS-CoV-2(N1) | 0.375906 |
SARS-CoV-2(N2) | 0.357839 | |
脱脂綿-Direct capture法(C+DC法) | CDC_N1 [copies/sample] | 0.420688 |
CDC_N2 [copies/sample] | 0.312676 | |
Phi6 | 0.332197 | |
SARS-CoV-2(N1) | 0.274482 | |
SARS-CoV-2(N2) | 0.320254 |
本実証事業では、SARS-CoV-2に加えて、20自治体のうち18自治体がPMMoVを、2自治体(甲府市、宇部市)がPhi6を下水中のバイオマーカーとして測定した。PMMoV,Phi6は共に糞便や人口のバイオマーカーとして下水サーベイランスはよく使われており、コロナの流行度や都市によらず、ある程度一定範囲の濃度であると考えられる。そこで、早速SARS-CoV-2の結果をみたいところだが、なにかおかしなデータがないかをチェックするために、まずは多くの自治体がバイオマーカーとして使ったPMMoVの測定データを見てみる。
手法が同じであればPMMoV濃度は自治体を超えてある程度の範囲に収まることがわかる。また、PEG沈殿法は沈殿物抽出法やEPISENS-S方よりPMMoV濃度が高く出やすいことがわかる。
なお、この結果は既に報告されている独立した研究による結果と同様である。
参照:
df[df['target'].str.contains('Phi6')]['city'].unique()
array(['甲府市', '宇部市'], dtype=object)
fig, ax = plt.subplots(figsize=(8, 4))
sns.boxplot(x="city", y="results",
hue="検査手法名", # hueで"検査手法名"ごとにまとめた統計にして、と言っている。採水方法"、"採水地点種別"ごとの統計にもこのパラメータを変えるだけで簡単にできる。
data=df[df['target'].str.contains('PMMoV')],ax=ax)
ax.set(yscale ='log')
ax.set_xlabel('')
ax.set_ylabel('[copies/L]')
ax.set_title('各自治体におけるPMMoV測定値の箱髭図')
ax.set_xticklabels(ax.get_xticklabels(),rotation=90);
fig, ax = plt.subplots(figsize=(8, 4))
sns.histplot(data=df[df['target'].str.contains('PMMoV')], x="results", hue="検査手法名", log_scale=True,kde=True,ax=ax)
ax.set_xlabel('[copies/L]')
ax.set_title('各手法におけるPMMoV測定値のヒストグラム')
Text(0.5, 1.0, '各手法におけるPMMoV測定値のヒストグラム')
データクリーニング、 データのQA/QCが終わったので、いよいよSARS-CoV-2のデータを見てみる。
下水サーベイランスデータにはばらつきがあるも、(週ごとの)時間平均・(自治体ごとの)空間平均を加えると長期的なトレンドが追いやすくなった。
# WW concentration w.r.t. that in August
# 手法を超えて濃度を同じスケールで比較することはできないので、8月平均と比較した変化量という新しいデータを'results_change'という名前で作る。
df_tmp = df[(df['date']>='2022-08-01') & (df['date']<'2022-09-01')]
df_ref = df_tmp.groupby(by=['city','loc','protocol_num','target']).median()[['results']].reset_index().rename(
columns={'results':'results_ref'})
df = df.merge(df_ref,on=['city','loc','protocol_num','target'],how='left')
df['results_change'] = df['results']/df['results_ref']
# 自治体(自治体内の複数サンプリングポイントを空間平均)及び週(同じ週で取った複数日時の時間平均)及び自治体ごとの平均を計算
df_ave = df.copy()
df_ave['week_start_date']=df_ave['date']-df_ave['date'].dt.weekday.astype('timedelta64[D]')
df_ave = df_ave.groupby(by=['city','loc','week_start_date','採水地点種別','target']).mean().reset_index()
# すべての自治体('処理場','ポンプ場'のデータのみ)の自治体・週平均のデータを見てみる。プライマーごとのデータを色(hue)を変えて表示。
#sns.set_style("white")
sns.set(font='IPAexGothic',style='white')
loc_type=['処理場', 'ポンプ場']
df_plot = df_ave[df_ave['target'].str.contains('SARS')].query('採水地点種別==@loc_type')
g = sns.relplot(
data=df_plot,
kind="line", # ラインプロットをつくる(このファンクションでは、95%CIがデフォルトで網掛けで表示される)
x="week_start_date", # 横軸は日付(厳密には週平均データをつくつために週の初めの日付けに統一したもの)
y="results", # 縦軸はSARS-CoV-2濃度
hue="target", # ターゲット(プライマーの種類)ごとに線の色を変える
col='city', # 市ごとに(各カラムに)プロット作る
col_wrap=4, # カラムは4つごとに行を変える
facet_kws={'sharey': False, 'sharex': True}) # 縦軸はプロットする範囲を統一しない、横軸はプロットする範囲を統一する
#plt.legend(loc='lower center',bbox_to_anchor=(1,0.5), borderaxespad=0)
for axes in g.axes.flat:
_ = axes.set_xticklabels(axes.get_xticklabels(), rotation=30)
g.fig.subplots_adjust(top=0.95) # adjust the Figure in rp
g.fig.suptitle('自治体ごとのSARS-CoV-2濃度(週及び自治体ないで時間、空間平均)',fontsize=14);
/var/folders/9b/3rq5w84x00j6qf03q8d3wngr0000gn/T/ipykernel_1591/3790666873.py:18: UserWarning: FixedFormatter should only be used together with FixedLocator
用いたデータは必ずしもスケールが揃っていない(例:新規感染者数の報告が市全体のこともあれば、処理区域内のみのものもある;ウイルス濃度が手法Aに基づくものもあれば手法Bに基づくものもある)ため、データを直接比較することはできない。相互に比較できるようにするため、下水中ウイルス濃度および新規感染者数のデータそれぞれに対し、”8月平均と比較した変化量”という新しいデータ作る。8月を選んだのは、i)8月には多くの自治体が下水サーベイランスを開始していた、ii)8月は全数検査が続いていたため、データのノイズが少ないと思ったから。
(各自治体のデータ特性を無視して一色単にデータ解析してみても)ばらつきはあるが、相対的感染者数と相対的下水ウイルス濃度に1:1の相関がありそうだ。さらに詳しい分析や結果の解釈には各自治体のデータ特性を理解する必要がある。
# Cases w.r.t. that in August
# スケールの異なるデータを直接比較することはできないので、8月平均と比較した変化量という新しいデータを'cases_7d_ave_ref_change'という名前で作る。8月を選んだのは、8月は全数検査が続いていたため、データのノイズが少ないと思ったから。
df_tmp = df_epi2[(df_epi2['date']>='2022-08-01') & (df_epi2['date']<'2022-09-01')]
df_epi_ref = df_tmp.groupby(by=['実証主体名']).median()[['cases_7d_ave']].reset_index().rename(
columns={'cases_7d_ave':'cases_7d_ave_ref'})
df_epi2 = df_epi2.merge(df_epi_ref,on=['実証主体名'],how='left')
df_epi2['cases_7d_ave_ref_change'] = df_epi2['cases_7d_ave']/df_epi2['cases_7d_ave_ref']
df_ave2 = df_ave.groupby(by=['city','week_start_date','採水地点種別','target']).mean().reset_index()
df_ave_merge = df_ave2.merge(df_epi2[['実証主体名','date','cases_7d_ave','cases_7d_ave_ref_change']],
left_on=['city','week_start_date'],right_on=['実証主体名','date'],how='left')
# 8月平均に比べた相対的感染者数と相対的下水ウイルス濃度の比較(自治体ごと)
# 自治体によって相関やばらつきが出てきそうか確認
# ばらつきはあるが、まあまあ1:1の相関が見える。
loc_type=['処理場', 'ポンプ場']
df_plot = df_ave_merge[df_ave_merge['target'].str.contains('SARS')].query('採水地点種別==@loc_type')
df_plot.dropna(subset=["cases_7d_ave_ref_change", "results_change"],inplace=True)
g = sns.scatterplot(data=df_plot, x="cases_7d_ave_ref_change", y="results_change",
hue="city", # 異なる色で異なる自治体を表現
style="target", # 異なるマーカーで異なるターゲットを表現
alpha=0.6) # 透明度を上げて見やすくする
g.set(xscale ='log',yscale ='log')
sns.move_legend(g, "upper left", bbox_to_anchor=(1, 1))
g.plot([0.03,5],[0.03,5], 'k--', alpha=0.3, zorder=0) # 1:1のラインを追加
plt.figtext(np.log10(5.5),np.log10(3.3), '1:1ライン')
plt.xlabel('相対的感染者数(8月と比較)')
plt.ylabel('相対的下水ウイルス濃度 \n (8月と比較;週及び事業主体ごとに平均)')
plt.title('感染者数と下水中ウイルス濃度の相関');
# 8月平均に比べた相対的感染者数と相対的下水ウイルス濃度の比較(月ごと)
# 月によって相関やばらつきが出てきそうか確認
loc_type=['処理場', 'ポンプ場']
df_plot = df_ave_merge[df_ave_merge['target'].str.contains('SARS')].query('採水地点種別==@loc_type')
df_plot.dropna(subset=["cases_7d_ave_ref_change", "results_change"],inplace=True)
df_plot['month_start_date']=df_plot['date'].map(lambda x: x - relativedelta(day=1)).dt.date
g = sns.scatterplot(data=df_plot.sort_values('date'), x="cases_7d_ave_ref_change", y="results_change",
hue="month_start_date", # 異なる色で月を表現
alpha=0.6)
g.set(xscale ='log',yscale ='log')
sns.move_legend(g, "upper left", bbox_to_anchor=(1, 1))
g.plot([0.03,5],[0.03,5], 'k--', alpha=0.3, zorder=0)
plt.figtext(np.log10(5.5),np.log10(3.3), '1:1ライン')
plt.xlabel('相対的感染者数(8月と比較)')
plt.ylabel('相対的下水ウイルス濃度 \n (8月と比較;週及び事業主体ごとに平均)')
plt.title('感染者数と下水中ウイルス濃度の相関');
終わりに
コードの筆者遠藤はデータサイエンティストでもソフトエンジニアでもなく、Python歴は5年もないので、このコードはとてもエレガントとは言えないかもしれないが、この程度のプログラミングでなかなかパワフルなデータ解析が簡単にできてしまうことをわかっていただけたらありがたい。
また、データ解析を行う上で、以下のようなデータ解析の前段階も非常に重要であることを強調したい。
下水サーベイランスデータに関する具体的な考察はこのドキュメントの目的ではないが、このように多くの自治体が半年以上かけて取得した貴重なデータが活用され、下水サーベイランスの特性の理解が進み、今後ますます下水サーベイランスの導入が進むことを期待したい。