如何从vcf文件中批量提取一系列基因的SNP位点?
需求
客户的一个简单需求:
我有一批功能基因位点,想从重测序的群体材料中找到这些位点,如何批量快速获得?
示例文件
gene.txt
image.png
test.vcf
image.png
代码实现
run.sh
cat $1 |while read gene chr from todo#echo $chr $from $toif echo $2 |grep -q '.*.vcf.gz$';thenvcftools --gzvcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out $gene.$chr.$from-$to elif echo $2 |grep -q '.*.vcf$';thenvcftools --vcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out $gene.$chr.$from-$tofidone
运行sh run.sh gene.txt test.vcf
,或sh run.sh gene.txt test.vcf.gz
生成结果:
image.png
补充说明
以上代码中利用了vcftools工具,以及shell中读取每行文件的每个字段进行赋值。
vcftools还能提取某个具体位置的SNP:
vcftools --gzvcf test.vcf.gz --positions specific_position.txt --recode --out specific_position.vcf
specific_position.txt文件格式如下:
1 8420131 8910211 9034261 9496541 1018704
除了vcftools,bcftools和plink等工具也能实现类似的功能。
bcftools filter test.vcf.gz --regions 9:4700000-4800000 > out.vcf
但bcftools要求vcf必须是gz格式,如不是,则需要进行转化(直接用gzip不行):
bcftools view test.vcf -Oz -o test.vcf.gzbcftools index test.vcf.gz
需要格外注意的是,vcf中的染色体名称要和提取文件中的染色体名保持一致,如Chr1或chr1或1。
或者:
bcftools view -S keep.list test.vcf >sub_indv.vcf
keep.list可以是“染色体+具体位置”两列,也可以是“染色体+起始+终止”三列:
chr1 27639chr1 60383chr2 60469chr3 60516chr4 60534#或者chr1 1 1000chr1 2000 4500
在plink中,可以指定特定的样本(keep)或SNP(extract)。
指定样本提取:
plink --bfile file --noweb --keep sampleID.txt --recode --make-bed --out sample
sampleID.txt第一列为提取的样本Family ID,第二列为Within-family ID(IID)。
指定位点提取:
plink --bfile file --extract snp.txt --make-bed --out snp
snp.txt文件中一个SNP名称一行。
赞 (0)