根据VCF设计Marker序列
问题描述:如果有两个样品的测序数据,并通过GATK等上游分析得到了变异位点信息,现在我想要找任意两个样品的差异SNP,并提取该位置上下游50bp序列,用于设计引物,应该怎么做?
本文将分享一种基于R语言从VCF文件快速获得引物设计序列的方法,用于从变异位点的vcf文件中寻找两样品的差异位点,并寻找参考基因组指定位置上下游区间,设计Marker引物序列。
1. 加载软件包
library(tidyverse)library(vcfR)
首先加载vcfR和tidyverse包,分别用于vcf文件的读取和数据操作。
2. 设置参数
out_name<-str_replace(file_name,".vcf","_marker.csv")#设置输出文件名 ref<-"xxx_assembly.fa"#参考基因组序列位置 dir_samtools<-"~/miniconda3/envs/work/bin/samtools"#samtools安装位置file_name<-"xxx.vcf"#输入文件名
以上代码定义了输入和输出文件,以及samtools的位置,用于后续与参考基因组的交互检索。其中xxx.vcf文件是至少包含两个样品的变异信息数据,默认比较前两个样品。
3. 读取并合并数据
df<-cbind(as.data.frame(vcf@fix),as.data.frame(vcf@gt))vcf<-read.vcfR(file_name)
读取VCF文件并合并固定的信息和基因型信息,生成一个数据框,用于后续数据清洗。
4. 判断变异位点类型
for(iin1:nrow(df)){ if(df[i,10]==df[i,11]){ df$type[i]<-"same" }else{df$type[i]<-"diff"} }df$type<-NA
遍历每行SNP数据,通过比较特定列来判断两样品的变异位点是相同还是不同。
5. 提取不同位点和单点突变
filter_snp<-filter[grep("^s",filter$ID,value=F),]filter<-df[which(df$type=="diff"),]
然后筛选出不同的位点,并进一步提取单点突变(SNP),通过筛选算法只提取SNP,过滤插入缺失突变。
6. 生成中间变异位点信息
filter_snp$info<-str_c("[",filter_snp$REF, "/",filter_snp$ALT,"]")
这句代码生成了变异位点的信息,格式为“[参考基因型/替代基因型]”。
7. 获取参考序列函数
cmd<-str_c(dir_samtools,"faidx",ref,"",Chr,":",pos_a,"-",pos_b) tem<-system(cmd,intern=T) return(paste(tem[2:length(tem)],collapse="")) }get_seq<-function(Chr,pos_a,pos_b){
该函数使用samtools来获取参考基因组的序列信息,只需要输入染色体名称,起始位置和终止位置,就可以自动返回这段区域的序列信息。
8. 迭代获取参考序列信息
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) }for(iin1:nrow(filter_snp)){
这部分代码迭代遍历差异SNP,并使用get_seq
函数获取每个SNP附近的序列,以便于设计引物。
9. 输出结果
write.csv(filter_snp,file=out_name,quote=F)
最后,差异SNP及其周围序列的信息被保存为CSV文件,下载后就可以用于直接设计引物了。
总结
上述R脚本提供了一种方法来识别两个材料序列之间的差异位点,并设计marker引物序列,可用于生物学研究,有助于解放双手,早点下班哈哈哈哈。
完整代码如下:
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(iin1: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(iin1: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)library(vcfR)
本文由 mdnice 多平台发布
如果觉得《根据VCF文件探测差异SNP并设计引物序列》对你有帮助,请点赞、收藏,并留下你的观点哦!