根据变异位点设计引物序列
今天碰到一个新问题：假如有一个vcf文件储存了两个样品的变异位点基因型数据，每行代表一个位点，我现在想找出两样本差异的SNP位点，再把差异位点用[REF/ALT]的形式表示，然后将其在参考基因组上下游100bp的序列找出来放在差异位点前后位置，得到一个序列文本，用于设计引物。
解决思路
- 如何判断差异SNP？
通过循环判断两个样品的基因型信息实现，相同时为same，不同时为diff - 如何提取差异位点？
通过tidyverse系列函数filter实现筛选，只选取s开头的SNP位点 - 如何获取参考基因组某段序列？
通过samtools调用faidx功能实现序列查询 - 如何生成引物设计信息？
通过合并字符串生成最终结果
下面是详细的实现方法，可以批量对vcf文件的差异位点生成引物设计数据，测试环境为linux R4.2.3，支持云端计算，如有建议或者需要欢迎联系后台交流。
1. 加载所需的库
首先，需要加载两个R包：vcfR和tidyverse。这两个包提供了处理VCF文件和数据处理的功能。
library(vcfR) library(tidyverse)
2. 设置运行参数
在代码中设置了几个运行参数：
file_name：数据文件的文件名out_name：输出文件的文件名ref：参考基因组序列的位置dir_samtools：samtools软件的安装位置
file_name <- "xx.vcf" out_name <- "xx_marker.csv" ref <- "ref.fa" dir_samtools <- "/miniconda3/envs/work/bin/samtools"
3. 读取并合并数据
使用read.vcfR函数读取VCF文件，并将其转换为数据框。然后，使用cbind函数将固定信息和基因型信息合并到一个数据框中。
vcf <- read.vcfR(file_name) df <- cbind(as.data.frame(vcf@fix), as.data.frame(vcf@gt))
4. 判断变异位点类型
为了确定每个位点的变异类型，使用一个循环来遍历数据框中的每一行。如果两个样本的基因型相同，则将类型设置为"same"，否则设置为"diff"。
df$type <- NA for (i in 1:nrow(df)) { if (df$A[i] == df$B[i]) { df$type[i] <- "same" } else { df$type[i] <- "diff" } }
5. 提取不同类型的位点
使用which函数和逻辑向量来筛选出类型为"diff"的位点，并将结果保存在一个名为filter的新数据框中。
filter <- df[which(df$type == "diff"), ]
6. 提取单点突变
为了进一步分析，只提取具有"diff"类型且以字母"s"开头的位点，这通常表示单点突变。
filter_snp <- filter[grep("^s", filter$ID, value = FALSE), ]
7. 打印结果
使用str_c函数和print函数打印出找到的变异位点总数和差异SNP的个数。
print(str_c("结果：共找到变异位点 ", nrow(df), "个！其中包括差异SNP ", nrow(filter_snp), "个！"))
8. 生成中间变异位点信息
使用str_c函数将参考基因组中的参考序列信息与变异位点的REF和ALT信息进行拼接，然后将结果保存在filter_snp$info列中。
filter_snp$info <- str_c("[", filter_snp$REF, "/", filter_snp$ALT, "]")
9. 获取参考基因组序列信息
定义了一个名为get_seq的函数，用于获取参考基因组序列信息。该函数接受染色体名称(Chr)和起始位置(pos_a和pos_b)作为参数，并使用samtools软件从参考基因组中提取相应的序列信息。
get_seq <- function(Chr, pos_a, pos_b) { cmd <- str_c(dir_samtools, " faidx ", ref, " ", Chr, ":", pos_a, "-", pos_b) tem <- system(cmd, intern = TRUE) return(paste(tem[2:length(tem)], collapse = "")) }
10. 迭代获取参考序列信息
使用一个循环来迭代处理每个变异位点，调用get_seq函数获取变异位点附近的参考序列信息，并将结果保存在filter_snp$out列中。
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) }
11. 保存结果为csv文件
最后使用write.csv函数将结果保存为一个csv文件，到了这一步，只需把文件放到polymarker网站上提交一下就能得到引物序列了。
write.csv(filter_snp, file = out_name, quote = FALSE)
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
