2023年12月7日发(作者：)

如何从vcf文件中批量提取一系列基因的SNP位点？

需求

客户的一个简单需求：

我有一批功能基因位点，想从重测序的群体材料中找到这些位点，如何批量快速获得？

示例文件

代码实现

cat $1 |while read gene chr from to

do

#echo $chr $from $to

if echo $2 |grep -q '.*.$';then

vcftools --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$';then

vcftools --vcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out $gene.$chr.$from-$to

fi

done

运行

sh

，或

sh

生成结果：

补充说明

以上代码中利用了vcftools工具，以及shell中读取每行文件的每个字段进行赋值。

vcftools还能提取某个具体位置的SNP：

vcftools --gzvcf --positions specific_ --recode --out specific_

specific_文件格式如下：

1 842013

1 891021

1 903426

1 949654

1 1018704

除了vcftools，bcftools和plink等工具也能实现类似的功能。

bcftools filter --regions 9:470 >

但bcftools要求vcf必须是gz格式，如不是，则需要进行转化（直接用gzip不行）：

bcftools view -Oz -o

bcftools index

需要格外注意的是，vcf中的染色体名称要和提取文件中的染色体名保持一致，如Chr1或chr1或1。

或者：

bcftools view -S >sub_

可以是“染色体+具体位置”两列，也可以是“染色体+起始+终止”三列：

chr1 27639

chr1 60383

chr2 60469

chr3 60516

chr4 60534

#或者

chr1　　1　　1000

chr1　　2000　　4500

在plink中，可以指定特定的样本（keep)或SNP(extract)。

指定样本提取：

plink --bfile file --noweb --keep --recode --make-bed --out sample

第一列为提取的样本Family ID，第二列为Within-family ID（IID）。

指定位点提取：

plink --bfile file --extract --make-bed --out snp

文件中一个SNP名称一行。