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コンソールの出力結果を保存しないときは「保存しない」を選択します。
