VCF文件时基因组分析中最常见的文件类型,有时需要从中提取部分信息进行后续分析,在《vcf、plink文件格式互转》中我们已经提及了SNPs的提取方法:
#在file.txt中, snp名字作为一列,无header,输出格式为vcf
vcftools --gzvcf my.vcf --snps snps.txt --recode --recode-INFO-all --out filter.snp
在实际应用中,还有许多变化,这里列出两个例子:
1.根据染色体(CHR)和物理位置(PS)筛选SNP
有些情况下,文件中rs这一列时空的,因此需要结合chr和ps筛选SNP,那么这时候需要一个两列信息的.txt文件,代码如下
#Filter based on chr and position
vcftools --gzvcf sample.vcf --positions snps.txt --recode --recode-INFO-all --out sample.filter.snp
2.将基因型转为012矩阵
筛选出来的snp矩阵是以基因型的形式储存的,为了方便统计,有时候需要将其转换为012矩阵
#把genotype转为012矩阵
vcftools --vcf snp.recode.vcf --012 --out snp_matrix
3.在R软件中分析vcf文件
这里需要用到一个软件包:vcfR
library(vcfR)
a<-read.vcfR( "recode.vcf")
a中会包含三部分,
meta,metadata
fix:vcf文件的前七列;
gt:genotype
a1<-data.frame(a@fix)
a2<-data.frame(a@gt)
a2<-a2[,-1]
rownames(a2)<-a1$ID
提取基因型
a3<-extract.gt(
a,
element = "GT",
mask = FALSE,
as.numeric = T,#是否转换为数值型
return.alleles = F,#是否输出allel的基因型信息(例如,0/1或A/T)
IDtoRowNames = TRUE,
extract = TRUE,
convertNA = TRUE#是否将"." 转为 NA
)