]遺伝子多型解析を数か月でできるようになるために〜第2回
- 『Applied Statistical Genetics With R』の読み会シリーズの目次
- 2011/05/27
- 第１章
- 遺伝子とは
- スタディデザイン
- 機能性を標的にするか・代理を標的にするか
- 機能性多型を標的に
- 遺伝子なり領域を標的にしてＬＤで
- 標的を定めるか、定めないか
- ターゲット
- ゲノムワイド連鎖解析・ゲノムワイド関連解析等の仮説フリー
- Genome情報とその他の中間形質情報との違い
- Genome情報と形質との関連を調べることと、中間形質と形質を調べることの違いについて
- phenotypeとcovariates
- phenotypeが複数の基準からなるときのことなどをディスカッション
- 実習
- ウェブ上のテキストファイルを直接読みこみ
fmsURL <- "http://people.umass.edu/foulkes/asg/data/FMS_data.txt" fms <- read.delim(file=fmsURL,header=T,sep="\t")
- geneticsパッケージを使う
$`African Am` Number of samples typed: 43 (97.7%) Allele Frequency: (2 alleles) Count Proportion A 60 0.7 G 26 0.3 NA 2 NA Genotype Frequency: Count Proportion A/A 23 0.53 A/G 14 0.33 G/G 6 0.14 NA 1 NA Heterozygosity (Hu) = 0.4268126 Poly. Inf. Content = 0.3328711
install.packages("genetics") library(genetics)
- 色々関数
- head(),tail()
- attach()は、データベースをアタッチする
- そうすると、関数に「何か」を与えると、その名前を変数としていないかをattachedされたオブジェクトに関して探索してくれる
attach(fms) head(id) #idという変数をattachedなfmsに探しに行く head() # では、そもそも探しに行かない
- 欠損値など
> 関数<-c("is.null()","is.na()","is.nan()","is.finite()","is.infinite()","complete.cases()") > 対象<-c("NULL","NA","NaN","有限か否か","無限か否か","欠測か否か") > data.frame("関数"=関数,"対象"=対象) 関数 対象 1 is.null() NULL 2 is.na() NA 3 is.nan() NaN 4 is.finite() 有限か否か 5 is.infinite() 無限か否か 6 complete.cases() 欠測か否
#http://www.people.umass.edu/foulkes/asg.html fmsURL <- "http://people.umass.edu/foulkes/asg/data/FMS_data.txt" fms <- read.delim(file=fmsURL,header=T,sep="\t") # read.table関数 #参考 #-------------------------------------------------------------------------- #fname<-c("read.csv","read.csv2","read.delim","read.delim2","read.fwf") #header<-c(rep("TRUE",4),FALSE) #sep<-c(",",";",rep("\t",3)) #dec<-c(".",",",".",",",".") #youto<-c("区切りがカンマの場合","区切りがコロンの場合（小数点がカンマ)","区切りがタブの場合","区切りがタブの場合(小数点がカンマ)","(1行の各欄の桁数widthを指定して読み込む)") #table_variation<-data.frame("関数名"=fname,"header"=header,"sep"=sep,"dec"=dec,"用途"=youto) #table_variation #--------------------------------------------------------------------------- #read.delim関数やread.csv関数は、read.table関数からの派生（ラッパー関数）。 #これらにおいては、header sep quote decは既定なので指定しなくてよい。 #--------------------------------------------------------------------------- #table_variation[3,] #--------------------------------------------------------------------------- #つまり、以下の3つは同じ意味。 #--------------------------------------------------------------------------- fms <- read.delim(file=fmsURL,header=T,sep="\t") fms <- read.delim(file=fmsURL) fms <- read.table(file=fmsURL,header=T,sep="\t") #--------------------------------------------------------------------------- #一番上の式は、「馬から落馬」している #ちなみに、 fms <- read.table(file=fmsURL,header=T) #でも実はいけるよう。sep=""とsep="\t"は指定しなくても読んでくれる。 #sep=","は読んでくれなかった。きちんと指定するか、read.csv関数を使用する。 #read.tableを使ったらheaderをちゃんと指定すること。 #TはTRUE、FはFALSE。（headerがある(TRUE)、ない(FALSE)） #いろいろ読み込んでみてください。 # example0 タブ区切り、ヘッダーあり # example0.1 コンマ区切り、ヘッダーあり # example0.2 スペース区切り、ヘッダーなし #仕組みさえわかれば、あとは自分が使いやすい形式を。 #表からデータを取り出してくる。 #Ageだけ取り出したい場合 #--------------------------------------------------------------------------- fms$Age #--------------------------------------------------------------------------- #1~25番目だけ #--------------------------------------------------------------------------- fms$Age[1:25] #--------------------------------------------------------------------------- #データフレームで data.frame(fms$Age)[1:25] #これはエラー。 #データフレームにしたら行と列の指定が必要。 #--------------------------------------------------------------------------- data.frame(fms$Age)[1:25,] data.frame(fms$Age,fms$id)[1:25,] #--------------------------------------------------------------------------- #いちいちfms$での指定が面倒なので #--------------------------------------------------------------------------- attach(fms) data.frame(id,actn3_r577x, actn3_rs540874, actn3_rs1815739, actn3_1671064, Term, Gender, Age, Race, NDRM.CH, DRM.CH)[1:20,] #--------------------------------------------------------------------------- #集計してくれるsummary関数 #--------------------------------------------------------------------------- GenoCount<-summary(actn3_rs540874) GenoCount #--------------------------------------------------------------------------- #もちろん、genotype以外でも集計してくれる。 #数字なら、最大・最小・中央値・平均値などをだしてくれる。 #--------------------------------------------------------------------------- summary(Race) summary(Age) #--------------------------------------------------------------------------- #Allele Frequencyを計算していく。 #--------------------------------------------------------------------------- attach(fms) GenoCount<-summary(actn3_rs540874) GenoCount NumbObs<-sum(!is.na(actn3_rs540874)) NumbObs #--------------------------------------------------------------------------- #数を合計してくれるsum関数。しかし... #--------------------------------------------------------------------------- sum(actn3_rs540874) class(actn3_rs540874) #--------------------------------------------------------------------------- #これではだめ。 #factorは数えられない。 #--------------------------------------------------------------------------- length(actn3_rs540874) #--------------------------------------------------------------------------- #数えられるが、NAの数も含まれている。 #そこで、NAだけ除いたベクトルを定義し直す、という方法も考えうる。 #[1から20を]取り出すなら #--------------------------------------------------------------------------- aaa<-actn3_rs540874[1:20] aaa #--------------------------------------------------------------------------- #[NA以外を]取り出すなら #--------------------------------------------------------------------------- aaa<-actn3_rs540874[!is.na(actn3_rs540874)] aaa #--------------------------------------------------------------------------- #その大きさ #--------------------------------------------------------------------------- length(aaa) NumbObs #--------------------------------------------------------------------------- #!は論理演算子の一つ。「〜でない」 #is.na()は、NAに対してTRUEを、それ以外に対してFALSEを与える。 #--------------------------------------------------------------------------- aaa<-c(1,2,3,NA,5,6,7,NA,9) is.na(aaa) !is.na(aaa) #--------------------------------------------------------------------------- #classはlogical #--------------------------------------------------------------------------- class(is.na(aaa)) #--------------------------------------------------------------------------- #logicalのベクトルにsum関数を当てはめると #--------------------------------------------------------------------------- sum(!is.na(aaa)) #--------------------------------------------------------------------------- #TRUEのものを合計してくれる。 #その他の論理演算子 #--------------------------------------------------------------------------- 関数<-c("is.null()","is.na()","is.nan()","is.finite()","is.infinite()","complete.cases()") 対象<-c("NULL","NA","NaN","有限か否か","無限か否か","欠測か否か") data.frame("関数"=関数,"対象"=対象) #--------------------------------------------------------------------------- #fmsDataの、Ageを合計してみる。 #--------------------------------------------------------------------------- sum(Age) #--------------------------------------------------------------------------- #これではNA。 #--------------------------------------------------------------------------- sum(!is.na(Age)) class(!is.na(Age)) #--------------------------------------------------------------------------- #括弧内がlogicalなので、TUREの数をカウントしてくれたのみ。 #ちなみに、Ageはもともとinteger（整数値）。 #numericとの違いは、 #--------------------------------------------------------------------------- aaa<-c(1.2,2.3,3.4,4.5,5.6) aaa class(aaa) as.integer(aaa) #--------------------------------------------------------------------------- #NAを取り除いたinteger vectorにして、合計する。 #--------------------------------------------------------------------------- sum(Age[!is.na(Age)]) #--------------------------------------------------------------------------- #もっとスマートな方法もありました。 #--------------------------------------------------------------------------- sum(Age,na.rm=TRUE) #--------------------------------------------------------------------------- #NAをremoveする。 #しかし、length関数は2つの引数を認めてくれないので、後者の方法は使えない。 # #閑話休題。 #allele frequencyの続き。 #--------------------------------------------------------------------------- GenoCount GenoFreq<-as.vector(GenoCount/NumbObs) GenoFreq #--------------------------------------------------------------------------- #この場合、左からAA,GA,GG,NAのFrequencyが計算。 #ちなみに、as.vectorにしなくても計算には支障ないですが、AAとかいう表示がついて回#ります。 #--------------------------------------------------------------------------- aaa<-c(GenoCount/NumbObs) aaa aaA<-(2*aaa[1]+aaa[2])/2 aaA #--------------------------------------------------------------------------- #普通に計算。 #--------------------------------------------------------------------------- FreqA<-(2*GenoFreq[1]+GenoFreq[2])/2 FreqA FreqG<-(2*GenoFreq[3]+GenoFreq[2])/2 FreqG #--------------------------------------------------------------------------- #楽をしたいのでpackageを使う。 #--------------------------------------------------------------------------- install.packages("genetics") library(genetics) #--------------------------------------------------------------------------- #--------------------------------------------------------------------------- Geno<-genotype(actn3_rs540874,sep="") summary(Geno) #--------------------------------------------------------------------------- #sep=""をきちんと設定しないと動いてくれません。注意。 #試し用 #--------------------------------------------------------------------------- aaa<-c(rep(c("CC","AA","CA","CA"),10)) bbb<-c("C/A","C/C","A/A","C/A","C/A","C/C","A/A","C/A","C/A","C/C","A/A","C/A","C/A","C/C") #--------------------------------------------------------------------------- #genotype関数は、次のように変換してくれる関数です。 #--------------------------------------------------------------------------- aaa genotype(aaa) class(aaa) #--------------------------------------------------------------------------- #真似をして作っても… #--------------------------------------------------------------------------- bbb<-c("C/A","C/C","A/A","C/A","C/A","C/C","A/A","C/A","C/A","C/C","A/A","C/A","C/A","C/C") summary(bbb) #--------------------------------------------------------------------------- #classを変えてみる #--------------------------------------------------------------------------- summary(as.factor(bbb)) #--------------------------------------------------------------------------- #これならいける #--------------------------------------------------------------------------- summary(as.genotype(bbb)) #--------------------------------------------------------------------------- #でも、↓ではgenotype属性がついてくれない。 #--------------------------------------------------------------------------- aaa as.genotype(aaa) #--------------------------------------------------------------------------- #やはり、genotype関数を使うのが便利そうです。 #その他 #--------------------------------------------------------------------------- hgdpURL<-"http://people.umass.edu/foulkes/asg/data/HGDP_AKT1.txt" hgdp<-read.delim(file=hgdpURL) hgdp[1:20,] #--------------------------------------------------------------------------- #--------------------------------------------------------------------------- vircoURL<-"http://people.umass.edu/foulkes/asg/data/Virco_data.csv" virco<-read.csv(file=vircoURL) #--------------------------------------------------------------------------- #problems #1.3 #--------------------------------------------------------------------------- library(genetics) attach(fms) summary(genotype(actn3_1671064,sep="")) #--------------------------------------------------------------------------- tapply(genotype(actn3_1671064,sep=""), Race, summary) #--------------------------------------------------------------------------- # tapply(A,B,f):AをBごとにfする #1.4 #--------------------------------------------------------------------------- library(genetics) attach(hgdp) summary(genotype(AKT1.C6024T,sep="")) tapply(genotype(AKT1.C6024T,sep=""), Geographic.area, summary) #--------------------------------------------------------------------------- #1.5 #--------------------------------------------------------------------------- library(genetics) attach(virco) summary(virco[,c("P1","P10","P30","P71","P82","P90")]) #--------------------------------------------------------------------------- MutNumVirco<-colSums(virco[,c("P1","P10","P30","P71","P82","P90")]!="-") Length<-apply(virco[,c("P1","P10","P30","P71","P82","P90")],2,length) MutNumVirco/Length #---------------------------------------------------------------------------
