R语言中操作vcf文件的一个实例—vcfR包

参考资料链接

• https://grunwaldlab.github.io/Population_Genetics_in_R/analysis_of_genome.html
• https://github.com/grunwaldlab/Population_Genetics_in_R

安装vcfR

直接通过install.packages()函数就可以安装

install.packages("vcfR")

vcfR简介

vcfR is a package intended to help visualize, manipulate and quality filter data in VCF file.

vcf格式文件

利用各种测序数据检测变异的结果存储文件，简单划分vcf格式文件里的内容

• ‘#’号开头行
• 非#号开头行

vcfR包存储vcf文件分为三个部分

• ‘#’号开头行——meta
• 非#号开头行分为fix和gt两个部分

fix部分存储vcf文件中非#号开头行的前7列，分别是

• 染色体编号
• 碱基位置
• ID
• 参考碱基
• 变异碱基
• 质量值
• 是否过滤

gt 部分存储两部分内容

• format
• 样本基因型

本文的数据使用 pinf_sc50.vcf.gz 来自R包 pinfsc50

读入数据

vcf<-read.vcfR("pinf_sc50.vcf.gz")

可以通过@符号获取meta、fix、gt三部分数据

vcf@meta
vcf@fix
vcf@gt

另外一些操作数据的操作

queryMETA(vcf,element="AD")
getFIX(vcf)[1:2,1:7]
vcf@gt[1:5,1:7]
gt<-extract.gt(vcf,element="GT")
gt[1:2,1:7]

结合基因组注释文件对数据进行展示

install.packages("ape")
library(ape)
vcf<-read.vcfR("pinf_sc50.vcf.gz",verbose=FALSE)
dan<-ape::read.dna("pinf_sc50.fasta",format="fasta")
gff<-read.table("pinf_sc50.gff",sep="\t",quote="")
chrom<-create.chromR(name="Supercontig",vcf=vcf,
                     seq=dan,ann=gff,verbose=TRUE)
chrom
plot(chrom)

image.png

options(scipen = 999)
chromoqc(chrom)

image.png

chrom<-proc.chromR(chrom,verbose = TRUE)
plot(chrom)
chromoqc(chrom)

image.png
image.png

探索一下以上作图用到的数据

> head(chrom@var.info)
             CHROM POS    MQ  DP mask  n Allele_counts        He       Ne
1 Supercontig_1.50  41 51.30 174 TRUE 16          0,32 0.0000000 1.000000
2 Supercontig_1.50 136 52.83 390 TRUE 17          32,2 0.1107266 1.124514
3 Supercontig_1.50 254 56.79 514 TRUE 17          31,3 0.1608997 1.191753
4 Supercontig_1.50 275 57.07 514 TRUE 17          31,3 0.1608997 1.191753
5 Supercontig_1.50 386 57.40 509 TRUE 16          29,3 0.1699219 1.204706
6 Supercontig_1.50 462 58.89 508 TRUE 17          31,3 0.1608997 1.191753
> head(chrom@seq)
1 DNA sequence in binary format stored in a vector.

Sequence length: 6 

Base composition:
    a     c     g     t 
0.000 0.333 0.167 0.500 
(Total: 6 bases)
> head(chrom@ann)
                V1 V2   V3     V4     V5 V6 V7 V8
1 Supercontig_1.50  . exon  38234  39706  .  +  .
2 Supercontig_1.50  . exon  94775  98998  .  +  .
3 Supercontig_1.50  . exon 102033 105840  .  -  .
4 Supercontig_1.50  . exon 102033 105840  .  -  .
5 Supercontig_1.50  . exon 102033 105840  .  -  .
6 Supercontig_1.50  . exon 102033 105840  .  -  .
                                                     V9
1               NAME=hypothetical protein;ID=PITG_17064
2     NAME=conserved hypothetical protein;ID=PITG_17068
3 NAME=coatomer protein complex, putative;ID=PITG_17069
4 NAME=coatomer protein complex, putative;ID=PITG_17069
5 NAME=coatomer protein complex, putative;ID=PITG_17069
6 NAME=coatomer protein complex, putative;ID=PITG_17069
> head(chrom@win.info)
             CHROM window start  end length   A   C   G   T    N other genic variants
1 Supercontig_1.50      1     1 1000   1000 267 213 293 227    0     0     0       16
2 Supercontig_1.50      2  1001 2000   1000 283 206 309 202    0     0     0       15
3 Supercontig_1.50      3  2001 3000   1000 229 213 235 177  146     0     0       10
4 Supercontig_1.50      4  3001 4000   1000   0   0   0   0 1000     0     0        0
5 Supercontig_1.50      5  4001 5000   1000   0   0   0   0 1000     0     0        0
6 Supercontig_1.50      6  5001 6000   1000   0   0   0   0 1000     0     0        0

欢迎大家关注我的公众好
小明的数据分析笔记本