読者です 読者をやめる 読者になる 読者になる

NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

MiSeqデータのMHC領域へのマッピング(15):RによるVCF fileの操作 (1): SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)の計算

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"にすべて書いています。