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

NGSデータ解析まとめ

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

MiSeqデータのMHC領域へのマッピング(10):VCF fileをRに読み込む

2015-04-13

MiSeqデータのMHC領域へのマッピング(7):samtools mpileup, bcftoolsでSNPsを抽出する - NGSデータ解析まとめ

(1) VCF fileの修正(下準備)
まず、上記エントリで作成したVCF formatのファイル(MHC_test4.var.vcf)を、テキストエディタで開いてみます。

##fileformat=VCFv4.1
##samtoolsVersion=0.1.18 (r982:295)
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=DP4,Number=4,Type=Integer,Description="# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads">
##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same">
##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">
##INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">
##INFO=<ID=G3,Number=3,Type=Float,Description="ML estimate of genotype frequencies">
##INFO=<ID=HWE,Number=1,Type=Float,Description="Chi^2 based HWE test P-value based on G3">
##INFO=<ID=CLR,Number=1,Type=Integer,Description="Log ratio of genotype likelihoods with and without the constraint">
##INFO=<ID=UGT,Number=1,Type=String,Description="The most probable unconstrained genotype configuration in the trio">
##INFO=<ID=CGT,Number=1,Type=String,Description="The most probable constrained genotype configuration in the trio">
##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=PC2,Number=2,Type=Integer,Description="Phred probability of the nonRef allele frequency in group1 samples being larger (,smaller) than in group2.">
##INFO=<ID=PCHI2,Number=1,Type=Float,Description="Posterior weighted chi^2 P-value for testing the association between group1 and group2 samples.">
##INFO=<ID=QCHI2,Number=1,Type=Integer,Description="Phred scaled PCHI2.">
##INFO=<ID=PR,Number=1,Type=Integer,Description="# permutations yielding a smaller PCHI2.">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases">
##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	./mapping_results/DRR003760/DRR003760.sort.bam	./mapping_results/DRR003761/DRR003761.sort.bam	./mapping_results/DRR003762/DRR003762.sort.bam
6	3427377	.	ATTACTTA	ATTA	71.2	.	INDEL;DP=2;VDB=0.0061;AF1=1;AC1=6;DP4=0,0,0,2;MQ=60;FQ=-35	GT:PL:GQ	1/1:110,6,0:7	0/1:0,0,0:3	0/1:0,0,0:3
6	29829653	.	G	T	35.8	.	DP=13;VDB=0.0599;AF1=0.531;AC1=3;DP4=1,3,3,6;MQ=21;FQ=-8.43;PV4=1,7e-08,1,1	GT:PL:GQ	0/1:0,0,0:3	1/1:76,27,0:24	0/0:0,12,29:10
6	29829659	.	T	C	77.4	.	DP=13;VDB=0.0459;AF1=1;AC1=6;DP4=0,0,4,8;MQ=22;FQ=-32.7	GT:PL:GQ	1/1:0,0,0:8	1/1:88,27,0:34	1/1:25,9,0:16
6	29829682	.	T	C	10.9	.	DP=6;VDB=0.0254;AF1=1;AC1=6;DP4=0,0,1,3;MQ=21;FQ=-28.5	GT:PL:GQ	1/1:0,0,0:3	1/1:43,12,0:12	1/1:0,0,0:3
6	29829850	.	T	C	97	.	DP=12;VDB=0.0001;AF1=1;AC1=6;DP4=0,0,9,3;MQ=21;FQ=-28.9	GT:PL:GQ	1/1:130,36,0:38	1/1:0,0,0:4	1/1:0,0,0:4

VCFファイルの最初の方には、上のようにヘッダ行(#で始まるいくつかの行)があり、それが終わると、各SNPサイトの情報が書かれた行が始まります。

このうち、各SNPサイトの情報が書かれた行の一つ上の行(#が1個だけのヘッダ行)は、VCF formatの各列の項目を示しています。RでVCFファイルを読み込む場合、この行を1番上にヘッダとして入れたい。そこで、まずテキストエディタで開いたVCF fileにおいて、

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	./mapping_results/DRR003760/DRR003760.sort.bam	./mapping_results/DRR003761/DRR003761.sort.bam	./mapping_results/DRR003762/DRR003762.sort.bam

の最初にある"#"を削除します。次に、最後の3つの項目における"./mapping_results/DRR00376*/DRR00376*.sort.bam"における"./mapping_results/DRR00376*/"までの部分も、邪魔なので削除します。

ファイルを上書き保存し、閉じます。

(2) VCF fileをRに取り込む
まず、Rを起動します。Rを起動すると、Rコンソールが立ち上がりますが、最初に、上の
「その他」メニューから
「作業ディレクトリの変更」
を選択し、読み込みたいVCF fileが置いてあるディレクトリを選択します。

次に、Rコンソールで以下のコマンドを実行します。

MHC_test4_var <- read.table("MHC_test4.var.vcf", header=T, sep="\t")

このコマンドは、MHC_test4.var.vcfをRに取り込んで、変数”MHC_test4_var"に格納する、ということをしています。

エラーが出なければ、うまくデータが読み込まれているかどうかを以下のコマンドで確認します。

head(MHC_test4_var,20)

うまくいけば、VCF fileのSNPsのリストが、header付きで20行分表示されると思います。もっと表示したい場合は、引数の20を50などもっと大きくします(数字は表示する行数です)。

VCFファイルには何個の多型サイト(主にSNPs)が含まれているのか? これを知るためには、VCF fileの行数を数えます。以下のコマンドを入力します。

length(MHC_test4_var[,1])

"MHC_test4_var"の行数が表示されます。

特定の列のみを表示する場合、以下のようにします。
たとえば1,2,3,4,5,9,10,11,12列目だけを表示したい場合は、

MHC_test4_var2 <- MHC_test4_var[,c(1,2,3,4,5,9,10,11,12)]

また、特定の列のみを削除したい場合、たとえば"INFO"列が長くて邪魔なので除く場合は

MHC_test4_var3 <- MHC_test4_var[,colnames(MHC_test4_var) != "INFO"]

とします。

Rを終了するときは、

q()

と入力します。Rコンソールの出力結果を保存しないときは「保存しない」を選択します。