基因组分析: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
)
上一篇:012—JAVA8新特性(lambda表达式)


下一篇:Office-012 xlsx快捷操作