2015-06-22
ここでは、Rを使ってVCF fileに含まれる情報を整理する方法を考えます。具体的には、
- VCF fileに含まれる各sampleのgenotype情報から、SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)を求める
- VCF fileから必要な情報だけを抽出し、SNPsについての新しい表を作成する
- SNP positionとallele frequencyとの間の関係(散布図)を作成する
以上3つのことを行います。
VCF fileのRへの読み込み、および基本的な操作については、このエントリを参考して下さい。
(復習)VCF fileのRへの読み込み
読み込むVCF fileですが、前回Variant Effect PredictorでアノテーションをしたVCF file(MHC_test4.var.annotated.vcf)を読み込みます。読み込む方法はこのエントリを参考して下さい。
もし"MHC_test4.var.annotated.vcf"がない場合はこちらへ。
1. SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)の計算
VCF fileでは、sampleごとのgenotype情報は、Genotype fieldの列に書かれています。そこでは、
1/1: ALT alleleのホモ接合体
0/1: REFとALT allelesのヘテロ接合体
0/0: REF alleleのホモ接合体
を示しています(このエントリも参照)。これはこれでOKなのですが、この表示のままでは、対立遺伝頻度や遺伝子型頻度の計算には不便なので、これを以下のように変換します。
1/1 -> 2
0/1 -> 1
0/0 -> 0
また、genotype fieldの各列における、各alleleの統計的な評価(参照ページ)の数字やコンマも邪魔なので除きます。
(1) Genotyping fieldの変換(正規表現による)
まず、Rで以下のコマンドを入力します。
# V50 genotype1 genotype1 <- sub("1/1:[0-9]+,[0-9]+,[0-9]+:[0-9]+","2", MHC_test4_annotated[,50]) genotype1 <- sub("0/1:[0-9]+,[0-9]+,[0-9]+:[0-9]+","1", genotype1) genotype1 <- sub("0/0:[0-9]+,[0-9]+,[0-9]+:[0-9]+","0", genotype1) genotype1 <- as.numeric(genotype1)
ここでは、Rのコマンド"sub"と正規表現を使って、genotype fieldの列(50列目)の値を変換したベクトル(複数の数値からなる変数)"genotype1"を作成しています。最後のコマンドの"as.numeric"は、変換した値のデータ型を「文字列」から「数値」に変換する命令です。コマンドは1行ずつ入力しても良いですし、Rの「ファイル」-> 「新規文書」でエディタを開き、そこにペーストしてからコマンド全体を選択し、「command」を押しながら「return」を押すと、すべての行を上から順に実行します。
参考URL:
R.4.51. データ型変換 | R Financial & Marketing Library
同じことを51列目、52列目(別のsampleのgenotype fields)にも行い、それぞれgenotype2, genotype3という名前のベクトルを作ります。
(2) 変換したgenotypes(数値)をtable("MHC_test4_annotated")に結合する
まず、genotype1からgenotype3を、Rの"cbind"コマンドを使って結合し、一つのtable(m x nの表)を作成します。以下のコマンドを実行します。
# genotype1からgenotype3を表として結合 genotypes <- cbind(genotype1, genotype2) genotypes <- cbind(genotypes, genotype3)
次に、以下のコマンドを実行することで、このtableをVCF fileのtableである"MHC_test4_annotated"に結合します。
# table genotypesをもとのVCF fileに結合 MHC_test4_annotated <- cbind(MHC_test4_annotated,genotypes)
結合されていることを、headコマンドで確認します。
# 確認 head(MHC_test4_annotated,10)
参考URL:
R言語で統計解析入門: サンプルの結合(行追加)と変数の結合(列追加) 梶山 喜一郎
(3) ALT_alleleの「対立遺伝子頻度」を計算する
ここではALT alleleの対立遺伝子頻度を計算します。ALT alleleの対立遺伝子頻度は、変換したgenotypeの数値を使って
(ALT_allele_freq) = (genotype1 + genotype2 + … + genotypeN) / (2 x number_of_individuals)
で計算できます。ここではN = 3なので、以下のコマンドを入力します(ついでに結果を"MHC_test4_annotated"に結合しています)。
# allele頻度 // referenceと異なるSNP alleleの頻度を求める alt_allele_freq <- (genotype1 + genotype2 + genotype3) / 6 #allele freqを計算 # VCF fileに結果を結合 MHC_test4_annotated <- cbind(MHC_test4_annotated,alt_allele_freq)
(4) 各Genotypeごとの「遺伝子型頻度」を計算する
遺伝子型頻度の計算は、ちょっとややこしいです。各SNPごとに、特定のgenotypeの個体数を数える必要があります。とりあえず以下のコマンドを実行して下さい(ALT_alleleのホモ接合体の頻度をSNPごとに計算)。
# Tableで行ごとに特定のgenotype個数を数える # 結果はベクトルで # number of samples no_samples <- 3 # genotype 1/1 -> alt_homo freq_alt_homo <- c() x <- 0 for (i in 1:length(genotype1)) { x <- length(which(genotypes[i,]==2))/no_samples freq_alt_homo <- c(freq_alt_homo,x) }
ヘテロ接合体、REF alleleのホモ接合体も同様にして計算できます。
結果を結合してtableを作成、VCF fileに結合します。
# 結果を結合 // tableにする freq_genotypes <- cbind(freq_alt_homo, freq_hetero) freq_genotypes <- cbind(freq_genotypes,freq_ref_homo) # 結果をVCF fileに結合 MHC_test4_annotated <- cbind(MHC_test4_annotated,freq_genotypes)
長くなったので、続きは次のエントリで。ここまでのRコマンドはこちらの"genotype.R"にすべて書いています。