今回は先日書いたGWASのQCに関して、『Primer to Analysis of Genomic Data using R』の第3章のサンプルデータを使用して実施する。
はじめに
ここで使用するデータはヒトではなく、83頭の羊54,977SNPsを使用してトレーニングする。今回のSNPsはイルミナ社のSNP chipで測定されたものであり、データはGenome Studioで生成されたものである。
今回のデータセットには2つのファイルが含まれており、genotypes fileとmap fileである。一つ目のファイルには、全サンプル、全SNPsに対する遺伝子多型のデータが含まれていて、二つ目のファイルにはSNPの位置情報（染色体など）が含まれている。
今回は'RSQLite'と呼ばれるR上でSQL様なデータベースの構築ができるパッケージを使用している。詳細に関してはパッケージのPDFを参照されたい。
早速、はじめる。
1. 下準備（'RSQLite'を使用してデータセットの構築をする）
パッケージ'RSQLite'を読みこむ。
library(RSQLite)
空のDB「SNPsmall」を作成する。
con <- dbConnect(dbDriver("SQLite"),dbname="SNPsmall")
conの中に、snpmapとSNPをdbWriteTableのfunctionを使用して書き込む。
dbWriteTable(con, "snpmap", "SNPmap.txt", header=T, append=T, sep="\t") dbWriteTable(con, "SNP", "SNPsample.txt", header=T, append=T, sep="\t")
conおよびSNP、snpmapの内容を確認する。
con <- dbConnect(dbDriver("SQLite"),dbname="SNPsmall") #conの中にどんなtableが含まれているか確認する。 dbListTables(con) # [1] "SNP" "snpmap" #それぞれのtable（SNP, snpmap）内にどんなfields（列）が含まれているか確かめる。 dbListFields(con,"snpmap") # [1] "name" "chromosome" "position" dbListFields(con,"SNP") # [1] "snp" "animal" "allele1" "allele2" "x" "y" "gcscore"
sample idをSNPから取り出して確認する→ベクトルにして、ペーストする。
animids <- dbGetQuery(con, "select distinct animal from SNP") head(animids) # animal # 1 sample1 # 2 sample5 # 3 sample6 # 4 sample7 # 5 sample8 # 6 sample9 #animidsをベクトルに変換した animids<-as.vector(animids$animal) hold <- dbGetQuery(con, paste("select * from SNP where animal='", animids[1], "'", sep="")) dim(hold) # [1] 54977 7
SNP のidを取り出して、snpidsとする
snpids <- as.vector(dbGetQuery(con, "select distinct name from snpmap")[,1]) length(snpids) # [1] 54977
ここまででデータベースを正しく構築できた。ここから本題のQCに入る。（それぞれのデータフォーマットがあるのでここまでの作業はRSQLiteの使い方のトレーニングと捉えていただければ十分である。）
2. Quality Control（Genotype Calling and Signal Intensities）
今回の事例で使用しているillumina社のデータはA/Bで表記されるか、実際の塩基（A,T,G,C）で記載されることが多い。
（今回のデータはA/Bで記述されている）
con <- dbConnect(dbDriver("SQLite"),dbname="SNPsmall") snp <- dbGetQuery(con, paste("select * from SNP where snp ='", snpids[1], "'", sep="")) dim(snp) # [1] 83 7 head(snp) # snp animal allele1 allele2 x y gcscore # 1 250506CS3900065000002_1238.1 sample1 A B 0.833 0.707 0.8446 # 2 250506CS3900065000002_1238.1 sample5 A B 0.829 0.714 0.8446 # 3 250506CS3900065000002_1238.1 sample6 A B 0.816 0.730 0.8446 # 4 250506CS3900065000002_1238.1 sample7 B B 0.031 1.132 0.8446 # 5 250506CS3900065000002_1238.1 sample8 B B 0.036 1.146 0.8446 # 6 250506CS3900065000002_1238.1 sample9 B B 0.037 1.150 0.8446
通常snpデータはcharacterになっているので、factorに変換する。
snp$allele1 <- factor(snp$allele1) snp$allele2 <- factor(snp$allele2) summary(snp$allele1) # A B # 36 47 summary(snp$allele2) # A B # 6 77
snpのallele1とallele2を結合してgenotypeという列を作り、snpのXとYのシグナル強度を遺伝子多型ごとに色と形を変更してプロットする。
XとYは二つのアレルの強度であり、イルミナ社特有のデータである。本来、片方のアレルのhomozygousはXが大きく、Yが小さい。もう一方のアレルのhomozygousはXが小さく、Yが大きいという結果が期待される。すなわち、heterozygousはその間をとることが予測される。
snp <- data.frame(snp, genotype=factor(paste(snp$allele1, snp$allele2, sep=""), levels=c("AA", "AB", "BB"))) plot(snp$x, snp$y, col=snp$genotype, pch=as.numeric(snp$genotype), xlab="x", ylab="y", main=snp$snp[1], cex.main=0.9) legend("bottomleft", paste(levels(snp$genotype), "(", summary(snp$genotype), ")", sep=""), col=1:length(levels(snp$genotype)), pch=1:length(levels(snp$genotype)), cex=0.7)
3. Quality Control（MAF）
少ないアレルの頻度（minor allele frequency）についてのQCを行う。
すべてのSNPsが多型を持っているわけでなかったり、非常に低頻度のアレルを持っているものは十分に想定される。しかしながら、この場合の関連解析には注意が必要であり、統計学的な検出力は極めて乏しい。そうした潜在的な影響を取り除くためにMAF < 0.5のSNPsについては解析から除外することがある。
#allele1とallele2を含んだベクトルを作成する。 alleles <- factor(c(as.character(snp$allele1), as.character(snp$allele2)), levels=c("A", "B")) #それぞれAとBのアレル頻度を求める。 summary(alleles)/sum(summary(alleles))*100 # A B # 25.3012 74.6988
MAFはA: 25.3%（もしくは0.253）と記載する。
4. Quality Control (Hardy-Weinberg Equilibrium)
HWE平衡は観察された遺伝子型とそこから推定される遺伝子型の両者の頻度をカイ二乗検定で確かめる。
#観察された遺伝子型（obs）を記述する。 obs <- summary(factor(snp$genotype, levels=c("AA","AB", "BB"))) print(obs) # AA AB BB # 6 30 47
まず、観察データからAとBのアレル頻度を求める。
hwal <- summary(factor(c(as.character(snp$allele1), as.character(snp$allele2)), levels=c("A", "B"))) hwal <- hwal/sum(hwal) print(hwal) # A B # 0.253012 0.746988 #ここで先ほどの推定の遺伝子多型頻度と観察数をかけて、推定された遺伝子型（exp）を記述する exp <- c(hwal[1]^2, 2*hwal[1]*hwal[2], hwal[2]^2)*sum(obs) names(exp) <- c("AA", "AB", "BB")
最後にイエーツの補正（連続性の補正）を行って、カイ検定を実施する。
xtot <- sum((abs(obs-exp)-c(0.5,1,0.5))^2/exp) pval <- 1-pchisq(xtot, 1) print(pval) # [1] 0.9136328
p値が0.91であることから、HWEから逸脱していないという帰無仮説は棄却されなかった。
まとめ
今回はイルミナ社のSNPデータに基づいた形式のQC例になったが、参考になるコマンドは多くある。
その他のデータ形式でも対応できるようにしておく必要がありそう。
（Chapter3. 4. 2まで終了）
20160121
RF
