R data frameの初めの100行分だけデータを抽出

R data frameの初めの100行分だけデータを抽出

df.LineList 数万行のデータで、これを使って集計したい。スクリプトを作成中はトライ&エラーのところがあって、ちょっとスクリプトを書いては試しを繰り返す。その、試しのスクリプトが機能するかを実行するたびに待ち時間が大きい。そこで、一部だけテスト用に抜き出したい、という場面です。

抜き出したdata frame をdf.testDATA へ代入します

#テスト用にはじめ100行のみのデータ
df.testDATA <- df.LineList[1:100,]

 

 

Rのdata frameに、条件によって異なる値をとるカラムの追加

下図のようなdata frame (df.LineList)があります。【ir_AE】というカラムは、 -1 の場合irAE、0の場合非irAEという情報です。

このir_AEカラムを見ながらirAEかどうかの情報を”IRAE”か”nonIRAE”かという文字列を持つカラム【IFirAE】へ書き込みたいという場合のスクリプト。dplyrパッケージのmutate if_else を使います。

library(dplyr)
df.LineList <- mutate(df.LineList, IFirAE = if_else(ir_AE == -1, true = "IRAE", false = "nonIRAE"))

実行すると次のようになります。(右端のIFirAE列が追加されている)

 

一般化したら次のようになります。(マニュアル

if_else(condition, true, false, missing = NULL, ..., ptype = NULL, size = NULL)

condition

論理ベクトル

true, false

条件の TRUE および FALSE 値に使用するベクトル。 true と false の両方が条件のサイズにリサイクルされます。 true、false、および missing (使用されている場合) は、共通の型にキャストされます。

missing

NULL でない場合は、条件の NA 値の値として使用されます。 true と false と同じサイズと型の規則に従います。

これらのドットは将来の拡張用であり、空にする必要があります。

ptype

必要な出力タイプを宣言するオプションのプロトタイプ。 指定した場合、これは true、false、および missing の共通タイプをオーバーライドします。

size

希望の出力サイズを宣言するオプションのサイズ。 指定した場合、これは条件のサイズをオーバーライドします。

Value

条件と同じサイズ、および true、false、および missing の共通型と同じ型を持つベクトル。
条件が TRUE の場合は true の一致値、FALSE の場合は false の一致値、NA の場合は欠落値の一致値 (指定されている場合)、それ以外の場合は欠落値が使用されます。

感想

SQLならチョイとできるデータの扱いもRだと不慣れで調べながらになるのがもどかしい。もちろんSQLとRを行ったり来たりすればいいのでしょうが、チョイと手間がかかるのと、Rのデータをデータベースへ読み込ませ、それをSQLで加工して、データベースからデータを吐き出して、Rへ読み込ませるというのは、ゴミが入る危険性をはらんでいるのでできれば避けたい作業なのです。

library(sqldf)
df_LineList <- df.LineList
df.new <- sqldf('SELECT *, case
                  when ir_AE == -1 then "IRAE" 
                  when ir_AE == 0 then "nonIRAE" 
                  else "" 
                  end as IFirAE FROM df_LineList')

(でも、Rのsqldfでやろうとすると、case文の文法が私の馴染んでいるやつと結構違う!

Rでマルチバイト文字の読み込みのエラー

マルチバイト文字の読み込みエラー

Rで集計しようとcsvのファイルを読み込ませようとしたところマルチバイト文字の読み込みのエラーになりました。’invalid multibyte string at…’の個所。

> JADERdata <- read.csv("0020 PhViDdata.csv", head=T)
Error in type.convert.default(data[[i]], as.is = as.is[i], dec = dec,  : 
  invalid multibyte string at '<8a><e9><90>}'こ

明示的に文字コードを教えてあげると良いと思い 、確かUTF-8だったはずとやったらはずれで

> JADERdata <- read.csv("0020 PhViDdata.csv", fileEncoding="utf-8", head=T)
Warning messages:
1: In read.table(file = file, header = header, sep = sep, quote = quote,  :
  invalid input found on input connection '0020 PhViDdata.csv'
2: In read.table(file = file, header = header, sep = sep, quote = quote,  :
  incomplete final line found by readTableHeader on '0020 PhViDdata.csv'

SJIS(シフトJIS)で読み込めました

> JADERdata <- read.csv("0020 PhViDdata.csv", fileEncoding="sjis", head=T)
>

 

 

 

京都出張 2024年02月

はじめに

出張は体力を消耗し、安宿に泊まる(会社の規定のホテル代金が、まともな所に留まるのを許してくれない)せいで不快なことが多いのですが、今回の京都出張は楽しく過ごすことが出来ました。以前お世話になった方々に学会でご挨拶できたというのもありますが、今回は宿泊を楽しむことができたというのが大きいです。
学会会場(京都みやこめっせ)までバスで行きました。写真はバス停からみやこめっせに行く途中の本山妙傳寺。
初日の教育プログラムが終わって。宿へ向かう途中。「ななじょうどおり」と「しちじょうどおり」があると知り合いから教えてもらいました。
京都駅前の地下街で、ヒラメと香味の「ぶぶ漬け」を美味しく頂きました。京都出身の知人から「ぶぶ漬け」ってホントに有るんだ!って旨のツッコミがありました。

調べたらぶぶ漬けは祇園の言葉で京都全域で使う表現ではないとのこと。

上方落語で、ぶぶ漬けを食って帰ろうとする大阪人と、ぶぶ漬け出さずに返そうとする京都人の駆け引きが出てくるお話で有名になった「ぶぶ漬け」の知名度を利用した商売かもしれない。

会社の規定で1泊あたりの宿泊料金の上限が決まっています。これは大学も他社も大体似たりよったり。そのためこれまで学会で京都に泊まる際はアパホテルだったり、カプセルホテルだったりを使っていました。

今回、私史上最高クラスのホテルに宿泊できました。

ホテルにはレコード(LPとか)を再生できる装置が共用のスペースに設置されていました。先に来てたお客さんが私の好みとは違う曲をかけていたので。写真を撮って。会釈だけして出てきました。
ホテルのコーヒーのセルフサービス(自分で淹れる)場所。部屋ごとにパックのコーヒーが置いてあるホテルは多いのですが、このホテルは各部屋にはコーヒーはなく、共用のスペースで。ミルで豆を挽いて、ドリップする様になってました。

 

淹れ終わって気がついたのですが、カップを計量計に載せて。豆毎に決められた重量になるまでお湯を注ぐように指示が貼ってありました。
豆の種類は3種類。私は薫りの良さげな左端の豆を使いました。お代わりを入れようと2度目やってきたところ、豆がありません。夜9時でクローズだそうです。

ボードゲームをする場所の様です。この近くを通ると、時々二人組で(刺しで)ゲームをしている人たちがいらっしゃいました。

大浴場から宿泊棟へもどる途中の共用スペース。卓球台が置いてありました。ひょっとしてこれがいわゆる温泉卓球?

ホテル、いや、ポテルの入り口付近

右側が大浴場。日中は宿泊客以外の方々も銭湯のように使用することができるようでした。左側は酒類など販売している売店。以下、周辺の様子の写真。

ビール他のアルコールやおつまみ類も無料で楽しむことができます。(セブンのチキンはセブンイレブンで購入してきましたが、それ以外はホテルで提供された無料アルコール・おつまみ)

月の石・ドライフルーツ・オレンジジュース

朝食の炊き込みご飯、各種おかず

焼きおにぎりのお茶漬け、温野菜・漬物・黒豆納豆・豆腐梅肉

和菓子・杏仁ヨーグルト

おまけの学会弁当

きっと、今回の出張では太ったな

電動アシスト(YAMAHA PAS) リコールの部品交換

はじめに

5年ほど子供が使っている電動アシスト自転車、タイヤが減ってきたということでタイヤ交換を街の自転車屋さんに持ち込みました。そこで指摘されたのが、ホイールの線状のサビでした。私の自転車のサビは軽度でしたが、そのまま放置して進行した場合には最悪走行中に車輪が分解して、外れるほどの問題が潜んでいるもののようでした。メーカーのリコールが発表されているのでまずはリコールの手当をしてもらうようにということでした。タイヤ交換を依頼した自転車屋さんは、取引のないメーカー(別の自転車屋さんで購入したものでした)ということで、購入した自転車店へ持ち込むことになりました。その手続きの過程で、バッテリーもリコールが発表されていることが判明し、それぞれ対応することになりました。

 

各部を撮影

自転車の号機番号です。まずはこれの番号を元に、YAMAHAのホームページでリコール対象であることを確認しました。メーカーに登録していたはずなのに、なぜ連絡が来なかったのかは謎ですが、購入店は当時「ダイエー」に入っていた自転車屋でしたが、その後経営主体がかわり「イオン」に入っているイオンバイクと言う店に変わっています。販売店自体の経営が変わったかどうかはわかりませんが、そういう理由で顧客リストの面倒な部分は引き継がなかったため経由での連絡が来なかったのか?とも思ったりしています。また、ホームページで確認したところ私の名前で別の型番の自転車が登録されていました。当時の販売店が登録用のハガキを他の自転車と取り違えて渡していたのかもしれません。

レギュレーションで厳しく各種の記録を求められる仕事をしている身からすると、いい加減な仕事をしても問題にされないような業種なんだなぁと感じてしまいます。国土交通省もメーカーに目を光らせてるけれども、自転車レベルでは流通の先までリコールが届かないような状態で放置しているという事かもしれません。

5年も屋外で直射日光や風雨にさらされてきたため、印刷はほとんど読み取れません。

この写真にあるX0Tから始まる(写真ではほとんど判別できない)文字列を、YAMAHAホームページに入力して交換対象かどうかを確認するようになっていました。

交換対象かどうかは、上記のX0Tで確認しますが、交換を申し込む際にはこのVK300…の番号も入力を求められました。

リムの状態です。線状のサビが前後ともに出ています。リムを溶接した箇所がサビに弱かったということです。リム全体の素材は錆びないものだけれども、溶接箇所に使用した具材がサビを誘発するものなのでしょうかね。

車輪のリム

リムの交換は、販売店で行いましたが、申込みはホームページから行いました。しばらくすると、YAMAHAからお電話が有りまして、近くで交換してくれる販売店をいくつか紹介してもらいました。候補の店舗をお伝えしたら、YAMAHAの担当の方が販売店と話を通してくださって、その後販売店から連絡がありました。

数日後、交換用のリムが販売店に届いたということで販売店から電話がありました。

バッテリー

バッテリーはホームページから申し込みました。2週間ほどで自宅まで、新しいバッテリーと、バッテリー放電機が届きました。

開梱すると、説明のチラシが入っていました。発火することもあるということで、使用を中止して充電もしないようにということです。(もう5年も使ってきてしまったのだが、(;_;))

こんな感じで梱包されて届きました。

はじめに見たときは、「充電器」もついてきたのかと思いましたが、普通は見ることがない「放電器」です。輸送中に出荷しないようにこれで放電してから、梱包するようにということでした。

バッテリーを保護するように、バッテリーをエアクッションでつつみ、段ボールの中にさらに段ボールで保護するような形になっていました。当該メーカーがバッテリーからの火災を深刻なそして、実際に起こりうるリスクとして認識していることがうかがい知れます。

放電器の使用方法。バッテリーを充電器と同じ要領で差し込んで使用します。数時間「ブーン」という音がしていました。放熱用のファンが回っているようです。そして5年間、体感的にはほとんど劣化することなく働いてくれたバッテリーは残量ボタンを押しても何も点灯しない状態になりました。

返送方法です。購入したときの箱の、蓋の部分の閉じる順番を逆にすると、返送用の伝票がすでに貼られていました。この仕組みはありがたいです。

送られてきたラベルが外側についていましたが、内側の蓋にはすでに返送用の伝票が貼られていました。

集荷依頼をするようにということですが、近くのやまとに持ち込みました。

Rでサンプルサイズあるいはパワーを計算

はじめに

臨床研究をするためのプロトコールに書くような、厳密性のある見積もりは統計の専門家に相談して進める必要があるのですが、専門家に依頼するとそれなりに前提を整理してお膳立てしてからお願いすることになります。また、専門家にきちんと検討してもらうには時間も必要です。しかしそこに至る前の段階で、段取りを組み立てて行くために、厳密性を犠牲にして大まかな見積もりが欲しいという場面も少なからずあります。(細かい前提をすっ飛ばすので雑です、すみません)

コードです

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%")

 

結果です

i) コントロール群で50%, 実薬群で25% (と薬効やイベントを想定), 検出力 power は80%, p=0.005 (これは意図的に汎用されているP=0.05や0.01より小さくしました)

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とだいたい一致しました。

ii) コントロール群で50%, 実薬群で 40%  (同上)

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とだいたい一致しました。

 

ゲノムワイド関連解析(GWAS)のマンハッタンプロットManhattan plot

ゲノムワイド関連解析(GWAS)のマンハッタンプロット

アイキャッチ画像のような図を論文で見ていて、どうやって描いたもので何を意味するのか自分でもとデータから図を描きつつ見方を学ぼうという事で書籍を開きながら試してみました。その時のチョコっとした学びを書き留めた備忘録です。

#初心者によるやってみた系の記事

参考にした書籍

参考にしたのは「ゼロから実践する 遺伝統計学セミナー〜疾患とゲノムを結びつける」という本です。参考というか、このホームページに書いていることは、この本に書かれていることをたどっただけです。書籍にはCGYWINのインストール、plinkの使い方から、マンハッタンプロットまで丁寧に書かれていますので、この書籍を読めば私のホームページは見る価値があまりないかもしれません。

Rでのプロット

ここで、今回使用した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’側に不均衡が集中している様なのでそのあたりを拡大したりしました。

 

 

Translate »