根据VCF文件探测差异SNP并设计引物序列

根据VCF设计Marker序列

问题描述：如果有两个样品的测序数据，并通过GATK等上游分析得到了变异位点信息，现在我想要找任意两个样品的差异SNP，并提取该位置上下游50bp序列，用于设计引物，应该怎么做？

本文将分享一种基于R语言从VCF文件快速获得引物设计序列的方法，用于从变异位点的vcf文件中寻找两样品的差异位点，并寻找参考基因组指定位置上下游区间，设计Marker引物序列。

---

1. 加载软件包

library(vcfR)
library(tidyverse)

首先加载vcfR和tidyverse包，分别用于vcf文件的读取和数据操作。

2. 设置参数

file_name <- "xxx.vcf" #输入文件名
out_name <- str_replace(file_name,".vcf","_marker.csv") #设置输出文件名
ref <- "xxx_assembly.fa" # 参考基因组序列位置
dir_samtools <- "~/miniconda3/envs/work/bin/samtools" # samtools安装位置

以上代码定义了输入和输出文件，以及samtools的位置，用于后续与参考基因组的交互检索。其中xxx.vcf文件是至少包含两个样品的变异信息数据，默认比较前两个样品。

3. 读取并合并数据

vcf <- read.vcfR(file_name)
df <- cbind(as.data.frame(vcf@fix),as.data.frame(vcf@gt))

读取VCF文件并合并固定的信息和基因型信息，生成一个数据框，用于后续数据清洗。

4. 判断变异位点类型

df$type <- NA
for (i in 1:nrow(df)){
      if (df[i,10] == df[i,11]){
            df$type[i] <- "same"
      }else{df$type[i] <- "diff"}
}

遍历每行SNP数据，通过比较特定列来判断两样品的变异位点是相同还是不同。

5. 提取不同位点和单点突变

filter <- df[which(df$type=="diff"),]
filter_snp <- filter[grep("^s",filter$ID,value = F),]

然后筛选出不同的位点，并进一步提取单点突变（SNP），通过筛选算法只提取SNP，过滤插入缺失突变。

6. 生成中间变异位点信息

filter_snp$info <- str_c("[",filter_snp$REF,
                     "/",filter_snp$ALT,"]")

这句代码生成了变异位点的信息，格式为“[参考基因型/替代基因型]”。

7. 获取参考序列函数

get_seq <- function(Chr,pos_a,pos_b){
      cmd <- str_c(dir_samtools," faidx ",ref," ",Chr,":",pos_a,"-",pos_b)
      tem <- system(cmd,intern = T)
      return(paste(tem[2:length(tem)],collapse = ""))
}

该函数使用samtools来获取参考基因组的序列信息，只需要输入染色体名称，起始位置和终止位置，就可以自动返回这段区域的序列信息。

8. 迭代获取参考序列信息

for (i in 1:nrow(filter_snp)){
      seq_head <- get_seq(filter_snp$CHROM[i],
                          as.numeric(filter_snp$POS[i]) - 100,
                          as.numeric(filter_snp$POS[i]) - 1)
      seq_tail <- get_seq(filter_snp$CHROM[i],
                          as.numeric(filter_snp$POS[i]) + 1,
                          as.numeric(filter_snp$POS[i]) + 100)
      filter_snp$out[i] <- str_c(seq_head,filter_snp$info[i],seq_tail)
}

这部分代码迭代遍历差异SNP，并使用get_seq函数获取每个SNP附近的序列，以便于设计引物。

9. 输出结果

write.csv(filter_snp,file = out_name,quote = F)

最后，差异SNP及其周围序列的信息被保存为CSV文件，下载后就可以用于直接设计引物了。

总结

上述R脚本提供了一种方法来识别两个材料序列之间的差异位点，并设计marker引物序列，可用于生物学研究，有助于解放双手，早点下班哈哈哈哈。

---

完整代码如下：

library(vcfR)
library(tidyverse)

# 设置运行参数
file_name <- "xxx.vcf" #输入文件名，重测序提取的数据文件
out_name <- str_replace(file_name,".vcf","_marker.csv") #设置输出文件名，csv文件
ref <- "xx_assembly.fa" # 参考基因组序列位置
dir_samtools <- "~/miniconda3/envs/work/bin/samtools" # samtools安装位置

# 读取并合并数据
vcf <- read.vcfR(file_name)
df <- cbind(as.data.frame(vcf@fix),as.data.frame(vcf@gt))

# 判断变异位点类型
df$type <- NA
for (i in 1:nrow(df)){
      if (df[i,10] == df[i,11]){
            df$type[i] <- "same"
      }else{df$type[i] <- "diff"}
}

# 提取两份材料中不同的位点
filter <- df[which(df$type=="diff"),]
# 提取单点突变
filter_snp <- filter[grep("^s",filter$ID,value = F),]

print(str_c("结果：共找到变异位点 ",nrow(df),
            "个！其中包括差异SNP ",nrow(filter_snp),"个！"))
# 生成中间变异位点信息
filter_snp$info <- str_c("[",filter_snp$REF,"/",filter_snp$ALT,"]")

# 定义一个函数，获取参考基因组序列信息
get_seq <- function(Chr,pos_a,pos_b){
      cmd <- str_c(dir_samtools," faidx ",ref," ",Chr,":",pos_a,"-",pos_b)
      tem <- system(cmd,intern = T)
      return(paste(tem[2:length(tem)],collapse = ""))
}

# 迭代获取参考序列信息
filter_snp$out <- NA
for (i in 1:nrow(filter_snp)){
      seq_head <- get_seq(filter_snp$CHROM[i],
                          as.numeric(as.numeric(filter_snp$POS[i]) - 100),
                          as.numeric(as.numeric(filter_snp$POS[i]) - 1))
      seq_tail <- get_seq(filter_snp$CHROM[i],
                          as.numeric(filter_snp$POS[i]) + 1,
                          as.numeric(filter_snp$POS[i]) + 100)
      filter_snp$out[i] <- str_c(seq_head,filter_snp$info[i],seq_tail)
}

write.csv(filter_snp,file = out_name,quote = F)

本文由 mdnice 多平台发布