オンラインセミナーで行った赤坂インタシティーの模様(ほぼ写真のみ)
オンラインセミナーで行った赤坂インタシティーの模様(ほぼ写真のみ)
5年ほど子供が使っている電動アシスト自転車、タイヤが減ってきたということでタイヤ交換を街の自転車屋さんに持ち込みました。そこで指摘されたのが、ホイールの線状のサビでした。私の自転車のサビは軽度でしたが、そのまま放置して進行した場合には最悪走行中に車輪が分解して、外れるほどの問題が潜んでいるもののようでした。メーカーのリコールが発表されているのでまずはリコールの手当をしてもらうようにということでした。タイヤ交換を依頼した自転車屋さんは、取引のないメーカー(別の自転車屋さんで購入したものでした)ということで、購入した自転車店へ持ち込むことになりました。その手続きの過程で、バッテリーもリコールが発表されていることが判明し、それぞれ対応することになりました。
自転車の号機番号です。まずはこれの番号を元に、YAMAHAのホームページでリコール対象であることを確認しました。メーカーに登録していたはずなのに、なぜ連絡が来なかったのかは謎ですが、購入店は当時「ダイエー」に入っていた自転車屋でしたが、その後経営主体がかわり「イオン」に入っているイオンバイクと言う店に変わっています。販売店自体の経営が変わったかどうかはわかりませんが、そういう理由で顧客リストの面倒な部分は引き継がなかったため経由での連絡が来なかったのか?とも思ったりしています。また、ホームページで確認したところ私の名前で別の型番の自転車が登録されていました。当時の販売店が登録用のハガキを他の自転車と取り違えて渡していたのかもしれません。
レギュレーションで厳しく各種の記録を求められる仕事をしている身からすると、いい加減な仕事をしても問題にされないような業種なんだなぁと感じてしまいます。国土交通省もメーカーに目を光らせてるけれども、自転車レベルでは流通の先までリコールが届かないような状態で放置しているという事かもしれません。
5年も屋外で直射日光や風雨にさらされてきたため、印刷はほとんど読み取れません。
この写真にあるX0Tから始まる(写真ではほとんど判別できない)文字列を、YAMAHAホームページに入力して交換対象かどうかを確認するようになっていました。
交換対象かどうかは、上記のX0Tで確認しますが、交換を申し込む際にはこのVK300…の番号も入力を求められました。
リムの状態です。線状のサビが前後ともに出ています。リムを溶接した箇所がサビに弱かったということです。リム全体の素材は錆びないものだけれども、溶接箇所に使用した具材がサビを誘発するものなのでしょうかね。
リムの交換は、販売店で行いましたが、申込みはホームページから行いました。しばらくすると、YAMAHAからお電話が有りまして、近くで交換してくれる販売店をいくつか紹介してもらいました。候補の店舗をお伝えしたら、YAMAHAの担当の方が販売店と話を通してくださって、その後販売店から連絡がありました。
数日後、交換用のリムが販売店に届いたということで販売店から電話がありました。
バッテリーはホームページから申し込みました。2週間ほどで自宅まで、新しいバッテリーと、バッテリー放電機が届きました。
開梱すると、説明のチラシが入っていました。発火することもあるということで、使用を中止して充電もしないようにということです。(もう5年も使ってきてしまったのだが、(;_;))
こんな感じで梱包されて届きました。
はじめに見たときは、「充電器」もついてきたのかと思いましたが、普通は見ることがない「放電器」です。輸送中に出荷しないようにこれで放電してから、梱包するようにということでした。
バッテリーを保護するように、バッテリーをエアクッションでつつみ、段ボールの中にさらに段ボールで保護するような形になっていました。当該メーカーがバッテリーからの火災を深刻なそして、実際に起こりうるリスクとして認識していることがうかがい知れます。
放電器の使用方法。バッテリーを充電器と同じ要領で差し込んで使用します。数時間「ブーン」という音がしていました。放熱用のファンが回っているようです。そして5年間、体感的にはほとんど劣化することなく働いてくれたバッテリーは残量ボタンを押しても何も点灯しない状態になりました。
返送方法です。購入したときの箱の、蓋の部分の閉じる順番を逆にすると、返送用の伝票がすでに貼られていました。この仕組みはありがたいです。
送られてきたラベルが外側についていましたが、内側の蓋にはすでに返送用の伝票が貼られていました。
集荷依頼をするようにということですが、近くのやまとに持ち込みました。
臨床研究をするためのプロトコールに書くような、厳密性のある見積もりは統計の専門家に相談して進める必要があるのですが、専門家に依頼するとそれなりに前提を整理してお膳立てしてからお願いすることになります。また、専門家にきちんと検討してもらうには時間も必要です。しかしそこに至る前の段階で、段取りを組み立てて行くために、厳密性を犠牲にして大まかな見積もりが欲しいという場面も少なからずあります。(細かい前提をすっ飛ばすので雑です、すみません)
setwd("c://Documents/R/hogehoge/") library(pwr) # i) コントロール群で50%, 実薬群で25% (と薬効やイベントを想定) # 検出力 power は80%, p=0.005 (これは意図的に汎用されているP=0.05や0.01より小さくしました) out_i <- pwr.2p.test(h = ES.h(p1 = 0.50, p2 = 0.25), sig.level = 0.005, power = 0.80, alternative = "two.sided") plot(out_i, main = "コントロール群で50%, 実薬群で25%") # ii) コントロール群で50%, 実薬群で 40% (同上) # 検出力 power は80%, p=0.005 (これは意図的に汎用されているP=0.05や0.01より小さくしました) out_ii <- pwr.2p.test(h = ES.h(p1 = 0.50, p2 = 0.40), sig.level = 0.005, power = 0.80, alternative = "two.sided") plot(out_ii, main = "コントロール群で50%, 実薬群で 40%")
pwrライブラリーを使わずに、組み込み関数を用いると、図は出ませんが次のようになります。
power.prop.test(p1=0.50, p2=0.25, sig.level = 0.005, power = 0.80)
結果です
Two-sample comparison of proportions power calculation
n = 98.28936
p1 = 0.5
p2 = 0.25
sig.level = 0.005
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
図中に示されたn=98とだいたい一致しました。
pwrライブラリーを使わずに、組み込み関数を用いると、図は出ませんが次のようになります。
power.prop.test(p1=0.50, p2=0.40, sig.level = 0.005, power = 0.80)
この結果です
Two-sample comparison of proportions power calculation n = 657.4394 p1 = 0.5 p2 = 0.4 sig.level = 0.005 power = 0.8 alternative = two.sided NOTE: n is number in *each* group
こちらも図中の657とだいたい一致しました。
アイキャッチ画像のような図を論文で見ていて、どうやって描いたもので何を意味するのか自分でもとデータから図を描きつつ見方を学ぼうという事で書籍を開きながら試してみました。その時のチョコっとした学びを書き留めた備忘録です。
#初心者によるやってみた系の記事
参考にしたのは「ゼロから実践する 遺伝統計学セミナー〜疾患とゲノムを結びつける」という本です。参考というか、このホームページに書いていることは、この本に書かれていることをたどっただけです。書籍にはCGYWINのインストール、plinkの使い方から、マンハッタンプロットまで丁寧に書かれていますので、この書籍を読めば私のホームページは見る価値があまりないかもしれません。
ここで、今回使用したRのスクリプトを掲載しようと思っていたのですが、そのスクリプトというのが上記書籍を購入したらダウンロードできるというデータの中にあるスクリプトで、ごくわずか変更(ファイル格納ディレクトリやデータファイル名を変更)した程度なので、さすがに掲載したらまずそうだと思いとどまり、別の例で示します。
CRANから、Rのパッケージqqmanをインストールしてあるとして次になります。
library(qqman) manhattan(gwasResults, chr="CHR", bp="BP", snp="SNP", p="P" )
gwasResultsはパッケージqqman付属のデータで、次のようにSNPの名称・染色体・位置・表現型と当該SNPとの関連性がないと仮定した場合の確率P値の4項目のデータが16470行続きます。
で、このスクリプトの結果は次の図のようになります。
色が寂しいので、同じデータでウクライナを思わせるカラーで表現すると
# Load the library library(qqman) # Make the Manhattan plot on the gwasResults dataset manhattan(gwasResults, chr="CHR", bp="BP", snp="SNP", p="P", col = c("#0057B7", "#FFDD00"))
次のようになります
図の中のSNPを表示するには , annotatePval = 0.01 を追加します
# Load the library library(qqman) # Make the Manhattan plot on the gwasResults dataset manhattan(gwasResults, chr="CHR", bp="BP", snp="SNP", p="P", col = c("#0057B7", "#FFDD00"), annotatePval = 0.01)
3番染色体の集中していくつものSNPが高い-log10(p)値を示すのに頂上の変異しか名前が表示されないので annotateTop = FALSE というオプションをつけてみました。
# Load the library library(qqman) # Make the Manhattan plot on the gwasResults dataset manhattan(subset(gwasResults, CHR == 3), chr="CHR", bp="BP", snp="SNP", p="P" , col = c("#0057B7", "#FFDD00"), annotatePval = 0.001 , annotateTop = FALSE, main = "Chr 3")
3番目の染色体に関連のありそうなSNPが集中しているようなので3番目のsubsetを表示してみます。肝は subset(gwasResults, CHR ==3) の個所です。
# Load the library library(qqman) # Make the Manhattan plot on the gwasResults dataset manhattan(subset(gwasResults, CHR == 3), chr="CHR", bp="BP", snp="SNP", p="P" , col = c("#0057B7", "#FFDD00"), annotatePval = 0.001 , annotateTop = FALSE, main = "Chr 3")
気になるあたりをもう少し拡大します
# Load the library library(qqman) # Make the Manhattan plot on the gwasResults dataset manhattan(subset(gwasResults, CHR == 3), chr="CHR", bp="BP", snp="SNP", p="P" , col = c("#0057B7", "#FFDD00"), annotatePval = 0.001 , annotateTop = FALSE, main = "Chr 3", xlim = c(200,500))
上記の作図は単純で、各Chromosome上のそれぞれの位置でp値をプロットしただけのものなので、結局はp値を計算して大きな外れ値(-> association)を示す位置を目視で確認しやすくしたものだという事がわかりました。
肝となるのはp値の計算になります。それでは、p値をどうやって計算していたかというと、複数人のSNPのデータと、それぞれのサンプル(の人)の興味のある表現型のデータを用いて表現型と無関係なSNPか表現型とassociateしているSNPかを計算しています。
plink.exe (windows用)がインストールされていて、冒頭の書籍の付録のファイル 1KG_EUR_QC.bed, 1KG_EUR_QC.bim, 1KG_EUR_QC.fam とかのデータがあったとして、コマンドプロンプトで次を実施しました。
plink.exe --bfile 1KG_EUR_QC --out 1KG_QC_Pheno1 --pheno phenotype1.txt --logistic --ci 0.95
上記で作成されるファイルから、p値等のデータを抽出するのにawkを使いました。awkを使うのに、MacOSだとawkはデフォルトで入っているのですが、今回はwindowsを使いましたのでCYGWINで実行しました。
awk '{print $2"\t"($1*100000000000+$3)"\t"$12}' 1KG_QC_Pheno1.assoc.logistic > 1KG_QC_Pheno1.assoc.logistic.P.txt
この出力で生成される 1KG_QC_Pheno1.assoc.logistic.P.txt ファイルを上の方で記したRスクリプトのインプットとすることで、マンハッタンプロットを描くことができました。
さらなる練習用という事で「心筋梗塞1666症例および対照健常者3198名のGWAS」データが公開されており、使用についての制限がないためこちらのデータで試してみました。
(制限がない=NBDC非制限公開データとは、アクセスに制限を設けることなく、誰でも利用することが可能な公開データです。)
次のサイトからダウンロードしました。
https://humandbs.biosciencedbc.jp/
「ゼロから実践する 遺伝統計学セミナー〜疾患とゲノムを結びつける」に従った例でawkで行っっていた、行の抽出の部分はをsqlで行いました。
好みの問題かもしれませんが、windowsでlinux emulatorみたいなのを稼働させてawk走らせるより、Microsoft AccessのSQLで抽出するのが素直に開発しやすいと感じました。
SELECT Release_L12vsUniv06_AllSample.[#SNPID] AS SNP, Release_L12vsUniv06_AllSample.chr AS chr, Release_L12vsUniv06_AllSample.chrloc AS location, Release_L12vsUniv06_AllSample.ptrend AS p INTO [0010 GWAS_results] FROM Release_L12vsUniv06_AllSample;
で作図は次のようにしました。
library(readxl) library(qqman) GWAS_results <- read_excel("0010 GWAS_results.xlsx") manhattan(GWAS_results, chr="chr", bp="location", snp="SNP", p="p" , col = c("#0057B7", "#FFDD00"), annotatePval = 0.00001 , annotateTop = FALSE) manhattan(subset(GWAS_results, chr == 12), chr="chr", bp="location", snp="SNP", p="p" , col = c("#0057B7", "#FFDD00"), annotatePval = 0.00001 , annotateTop = FALSE, xlim = c(1.1e+8, 1.12e+08))
12番染色体の3’側に不均衡が集中している様なのでそのあたりを拡大したりしました。
解説した動画もありますので併せてご確認いただけると分かりやすいかと思います!https://t.co/D3QKwvNULS
— Hirosato Gamo (@hiro_gamo) April 27, 2023