Note of Pediatric Surgery

腸内細菌、R、ときどき小児外科

菌叢解析結果の比較 ( 5 ) 検定を一括で行い結果をcsvファイルにまとめる

というわけで、以前の記事をまとめるて、if構文で繋げると以下の様になります。

[splus] BTable=NULL ATable=NULL FTable=NULL for ( i in 3:ncol(X)) { BP=bartlett.test(X[,i],X$Test)$p.value if(BP<0.01){ SD<-"" }else if(BP<0.05){ SD<-"*" }else{ SD<-"n.s." } BResult=c(colnames(X)[i],BP,SD) BTable=rbind(BTable,BResult) if(BP>0.05){ ASummary = summary(aov(X[,i]~X$Test+X$Subject)) AResult = c(colnames(X)[i],ASummary1["Pr(>F)"][1,1],ASummary1["Pr(>F)"][2,1]) ATable = rbind(ATable,AResult) }else{ FP=friedman.test (X[,i],X[,1],X[,2])$p.value if(FP<0.01){ SD<-"" }else if(FP<0.05){ SD<-"*" }else{ SD<-"n.s." } FResult=c(colnames(X)[i],FP,SD) FTable=rbind(FTable,FResult) } } colnames(BTable)=c("Phylum","P.value","SD") write.csv(BTable,"BartlettTest_Result.csv",row.names=FALSE) colnames(ATable)=c("Phylum","Test$P.value","Subject$P.value") write.csv(ATable,"AOV_Result.csv",row.names=FALSE) colnames(FTable)=c("Phylum","P.value","SD") write.csv(FTable,"FriedmanTest_Result.csv",row.names=FALSE) [/splus]

できました!出力したファイルを添付しておくと、

BartlettTest_Result.csv

"Phylum","P.value","SD" "Bacteroidetes","0.532577365726786","n.s." "Firmicutes","0.0803489400883996","n.s." "TM7","0.637276102471321","n.s."

AOV_Result.csv

"Phylum","Test$P.value","Subject$P.value" "Bacteroidetes","0.540929440422521","0.00276660508206454" "Firmicutes","0.703002431390229","0.333099501465002" "TM7","0.687474809867674","0.44333185016966"

FriedmanTest_Result.csv

""

"FriedmanTest_Result.csv"に何も入っていないのは、バートレット検定ですべて有意差がなく、分散が等しいという結果になったので問題ないです。取り敢えずこれで完成です!これで被検者が何万人いても、何百種類といったspecies levelでも解析できるはずです。初めて書いた長いスクリプトなので、感慨深いです。

注意点としては、friedman.test()は欠損値は欠損値としてNAで入力しないとエラーが出て動きません。例えばBの試験だけ被検者1のサンプルがかけていた場合は、relative abundanceをの項目を各phylum毎に"NA"で埋めておいて下さい。僕はこれではまりました。

本当は二元配置分散分析とフリードマン検定の結果をまとめて1つのcsvファイルにしたかったのですが、二元配置分散分析は各項目ごとに有意差が出てしまうので、一応別々にしました。今回は検討実験で多群2要因でやったのでややこしい検定になりましたが、普通の介入試験の結果比較であれば、1つのcsvファイルにまとめて出力できると思います。

さて、これだけでは終われません。ここまでの解析ではBacteroidetesが被検者間で有意に変化しているということがわかりそうですが、多群間の比較では、P値をそのまま適応してはいけません。そこで、次回は多重検定について書いていきたいと思います。