【plink 和vcftools使用】从 VCF 文件中提取指定 SNP 的 REF/ALT 方法

从 VCF 文件中提取指定 SNP 的 REF/ALT 方法

---

方法一：基于 vcftools 直接提取

步骤说明

1. 准备位点文件 positions.txt

# 格式要求：两列，无表头，分别为染色体和物理位置（空格/TAB分隔）
# 示例内容：
1   3078680
2   3142683

2. 提取指定位点的 VCF 子集

vcftools --vcf raw.vcf --positions positions.txt --recode --out extracted

• --positions：指定目标位点文件
• 输出文件：extracted.recode.vcf

3. 提取 REF/ALT 信息

# 假设 VCF 头部有 39734 行注释（以 # 开头）
tail -n +39735 extracted.recode.vcf | cut -f 1,2,4,5 > genotype.list

• 关键点：
  • 使用 tail -n +N 跳过头部注释（N 需根据实际 VCF 文件调整）
  • cut -f 1,2,4,5 提取 CHROM, POS, REF, ALT 列

4. 结果示例

#CHROM  POS     REF     ALT
1       223     T       A
1       224     A       G

---

方法二：通过 Plink 转换后提取

步骤说明

1. 将 VCF 转为 Plink 格式

plink --vcf raw.vcf --const-fid -snps-only just-acgt --keep-allele-order \--biallelic-only strict --allow-extra-chr --make-bed --out vcf-plink \--chr-set 95

• --chr-set 95：适用于非人类物种（人类数据可移除此参数）
• --keep-allele-order：保持 VCF 中 REF/ALT 的原始顺序

2. 重命名 SNP ID（可选）

mv vcf-plink.bim vcf-plink-old.bim
awk '{OFS="\t"; $2=$1":"$4; print $0}' vcf-plink-old.bim > vcf-plink.bim

• 将 SNP ID 格式改为 染色体:位置（示例：1:3078680）

3. 准备 SNP 列表文件 snp.list

# 格式要求：一列，每行为 SNP ID（需与 .bim 中命名一致）
# 示例内容：
1:223
1:224

4. 提取目标 SNP

plink --bfile vcf-plink --extract snp.list --keep-allele-order \--make-bed --out this --chr-set 95 --allow-extra-chr

• 输出文件：this.bim 包含目标 SNP 信息

5. 结果示例（this.bim 格式）

CHROM  SNPID       CM  POS    ALT  REF
1      1:223       0   223    A    T
1      1:224       0   224    G    A

• 注意：Plink 的 .bim 文件中，REF/ALT 列顺序！

---

方法对比

特性	vcftools 方法	Plink 方法
优势	直接提取，无需格式转换	适合后续关联分析或质控
劣势	需处理 VCF 头部注释行	REF/ALT 顺序可能反转，需加参数限制keep-allele-order
适用场景	快速提取少量位点	需要结合 Plink 进行下游分析

---

注意事项

1. VCF 头部行数：
   使用 grep -n "^#" raw.vcf | tail -1 查看实际注释行数，替换 tail -n +N 中的 N。

2. REF/ALT 顺序：

   • VCF 格式：第 4 列为 REF，第 5 列为 ALT
   • Plink 的 .bim：第 5 列为 ALT，第 6 列为 REF

3. 染色体格式：
   非标准染色体（如 chr1）需在 Plink 命令中添加 --allow-extra-chr。