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で以下のコマンドを入力します。
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の表）を作成します。以下のコマンドを実行します。
genotypes <- cbind(genotype1, genotype2)
genotypes <- cbind(genotypes, genotype3)
次に、以下のコマンドを実行することで、このtableをVCF fileのtableである"MHC_test4_annotated"に結合します。
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"に結合しています）。
alt_allele_freq <- (genotype1 + genotype2 + genotype3) / 6
MHC_test4_annotated <- cbind(MHC_test4_annotated,alt_allele_freq)
(4) 各Genotypeごとの「遺伝子型頻度」を計算する
遺伝子型頻度の計算は、ちょっとややこしいです。各SNPごとに、特定のgenotypeの個体数を数える必要があります。とりあえず以下のコマンドを実行して下さい（ALT_alleleのホモ接合体の頻度をSNPごとに計算）。
no_samples <- 3
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に結合します。
freq_genotypes <- cbind(freq_alt_homo, freq_hetero)
freq_genotypes <- cbind(freq_genotypes,freq_ref_homo)
MHC_test4_annotated <- cbind(MHC_test4_annotated,freq_genotypes)
長くなったので、続きは次のエントリで。ここまでのRコマンドはこちらの"genotype.R"にすべて書いています。
