基因组分析：VCF文件中位点提取以及R软件中的分析

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
)

标签：vcf,VCF,--,基因组,012,snp,recode,txt,位点
来源： https://blog.csdn.net/sunfishtan/article/details/120433066