网站更换主机注意,wordpress 三栏制作,哈尔滨市建设厅网站,wordpress封装根据VCF设计Marker序列 问题描述#xff1a;如果有两个样品的测序数据#xff0c;并通过GATK等上游分析得到了变异位点信息#xff0c;现在我想要找任意两个样品的差异SNP#xff0c;并提取该位置上下游50bp序列#xff0c;用于设计引物#xff0c;应该怎么做#xff1f… 根据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(vcffix),as.data.frame(vcfgt)) 读取VCF文件并合并固定的信息和基因型信息生成一个数据框用于后续数据清洗。 4. 判断变异位点类型 df$type - NAfor (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$typediff),]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(vcffix),as.data.frame(vcfgt))# 判断变异位点类型df$type - NAfor (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$typediff),]# 提取单点突变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 - NAfor (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 多平台发布
