VET:一个基于R语言的VCF数据提取工具,支持按基因ID、物理位置、样品名称提取指定变异信息

简介: VET:一个基于R语言的VCF数据提取工具,支持按基因ID、物理位置、样品名称提取指定变异信息

VET:Vcf Export Tools

工具简介

VET是一个基于R语言开发的变异位点信息批量提取工具,主要功能是根据VCF数据集,按照基因ID、样品ID、变异位点ID等参数,实现批量提取,同时支持变异位点结构注释,一步搞定变异数据的快速提取。

########## WelCome to VCF Export Tools ###########
>>>>>>>>>>>>>>>> Design By BioNote <<<<<<<<<<<<<<<<<
可选参数:
[1]根据基因ID提取变异数据
[2]根据物理位置提取变异数据
[3]根据样品名称提取变异数据
[4]根据SNP名称提取变异数据
--------------------------------------------------
[INFO]第一个参数填选项,第二个参数填项目备注名称
[INFO]第三个参数选择是否过滤样本,默认为“N”
[INFO]第四个参数选择是否进行结构注释,默认为“N”
[INFO]运行方法 $ Rscript ./run.R 1 test Y Y
>>>>>>>>>>>>>>>> 程序版本:V 2.0.1 <<<<<<<<<<<<<<<<
##################################################

功能与应用

基因测序后经过上游分析得到的VCF文件储存了所有样本对应的所有变异信息,通常数据量非常大,在实际使用中需要根据情况对指定信息进行提取。目前已有vcftools或bcftools等工具能够实现上述操作,但是用的时候参数比较复杂,整个过程略显繁琐。


本工具集成了R、tidyverse、Python、vcftools、bcftools、snpEff等常用工具,开发便捷式流程实现批量操作。

主要应用场景是对大规模VCF原始数据集进行提取,支持多种个性化方式

  • 按基因ID提取指定基因内变异信息
  • 按材料名称提取某些材料变异信息
  • 按变异位点名称筛选指定变异信息
  • 按照物理位置提取指定区段内数据

支持流式操作,提取后筛选指定样品并对每个变异位点进行结构注释(判断错义突变、移码突变等),最终将结果文件打包生成压缩包。

使用方法

  • 第一步:输入待提取的信息
  • 第二步:运行Run.R脚本
  • 第三步:下载结果文件

可以批量操作,无需手写代码。

原理介绍

VCF是生信研究中储存基因变异信息的重要格式,通常经过测序上游分析后得到一份具有丰富信息量的vcf或者vcf.gz文件。

以“##”开头的行表示注释信息,一般记录了字段类型和历史命令,这部分相当于一个日志信息。剩下的数据部分类似一个表格,大体上每行是一个变异位点,每列是一个材料样本。

变异位点(简称SNP)是按照染色体上不同位置进行统计的,展示不同材料中在某个位置碱基差异。

按照突变的类型可以分为3种类型:

  • 缺失:Del,某些碱基不见了
  • 插入:Ins,新出现了某些碱基
  • 替换:SNP,单核苷酸多态性变异

一般常见的VCF文件主要由上述信息组成,对于大规模测序得到的多个样品合并VCF文件,可能包含几千万行×几千列(亿级数据量)。

在实际进行分析时,可能只需要考虑某几个基因或者是一小段区间内的变异数据,因此需要对VCF文件进行提取,只取出想要的一小部分,这个过程涉及到Linux下不同软件的相互配合。

VET 源代码

首先,建立项目文件夹并生成以下结构:

Aug  9 16:30 00_scripts
Aug 17 15:24 01_INPUT_GeneID.txt
Aug 17 16:09 01_out_byGeneID
Aug  9 18:30 02_INPUT_Postion.txt
Aug 10 11:25 02_out_byPostion
Aug 10 11:02 03_INPUT_SampleName.txt
Aug  4 15:13 03_out_bySampleName
Aug  4 15:31 04_INPUT_SNP.txt
Aug  4 15:14 04_out_bySNP
Aug  6 11:34 05_INPUT_filevcf.txt
Aug  6 11:33 05_out_bySnpEff

程序初始化

下面的代码为程序初始化过程,将会加载tidyverse等软件包,并读取重要参数,完成后将获得输出提示。

#!/usr/local/bin/Rscript
# VCF Export Tools 基因型变异数据批量提取工具,快捷提取VCF文件
# 依赖软件:Python、bcftools、tidyverse、snpeff
suppressPackageStartupMessages(library("cli"))
suppressPackageStartupMessages(library("tidyverse"))
cli::cli_text("########## WelCome to VCF Export Tools ###########
 \n >>>>>>>>>>>>>>>> Design By Jewel <<<<<<<<<<<<<<<<<
 \n可选参数:
 \n\t[1]根据基因ID提取变异数据
 \n\t[2]根据物理位置提取变异数据
 \n\t[3]根据样品名称提取变异数据
 \n\t[4]根据SNP名称提取变异数据
 \n--------------------------------------------------
 \n[INFO]第一个参数填选项,第二个参数填项目备注名称
 \n[INFO]第三个参数选择是否过滤样本,Y为过滤指定样本
 \n[INFO]第四个参数为'Y'时将对vcf文件进行变异结构注释
 \n[INFO]例如 $ ./run.R 1 test Y Y
 \n>>>>>>>>>>>>>>>> 程序版本:V 2.0.1 <<<<<<<<<<<<<<<<
 \n ##################################################")
 args <- commandArgs(T)
if(length(args)!=4){stop("参数输入有误,请检查输入格式,示例“./Run.R 1 jobname Y/N Y/N")}
# CONFIG SETTING:
db_file <- "wgs_all.vcf.gz" # 设置数据库名称
db_name <- "WGS"
# 程序初始化,删除上次输出结果文件----------
OPT <- args[1] # 程序子选项
JOB <- args[2] # 项目备注信息
SAM <- args[3] # 是否过滤样本
EFF <- args[4] # 是否结构注释
print(str_c("INFO   当前选择的数据库为:",db_file))
print(str_c("INFO   当前项目名称为: ",OPT," <-> ",JOB))
print(str_c("INFO   是否对样品进行过滤(Y为过滤指定样本,否则不过滤):"),SAM)
print(str_c("INFO   是否对vcf进行结构变异注释(Y为进行注释,否则不注释):"),EFF)
system("rm -rf ./01_out_byGeneID/*")
system("rm -rf ./02_out_byPostion/*")
system("rm -rf ./03_out_bySampleName/*")
system("rm -rf ./04_out_bySNP/*")
system("rm -rf ./05_out_bySnpEff/*")
cli::cli_text("INFO   系统输出文件夹初始化完成")

根据基因ID提取变异信息

根据输入的参数进行判断,如果选项为1,则执行下面的步骤,主要调用Python程序进行信息检索,并由bcftools工具批量提取变异信息,若需要根据指定样品进行过滤,则利用view功能对样品进行筛选,最后生成结果压缩文件。

if (OPT == "1"){
  cli::cli_text("INFO   待提取的基因ID如下,将自动自取上下游3000bp内的变异数据")
  id <- read.table("./01_INPUT_GeneID.txt",header = F)
  print(id$V1)
  cli::cli_text("INFO   基因ID信息整理完毕,接下来开始检索物理区间")
  system("Rscript prefix_gene_filter.R ./01_INPUT_GeneID.txt")
  cli::cli_text("INFO   接下来执行Python脚本调用bcftools提取基因变异信息")
  system(str_c("python bcftools_view_filiter_Chr.py --input ./00_scripts/id.txt --vcf ./",
               db_file))
  cli::cli_text("INFO   提取完成,对结果进行打包压缩")
  if (SAM == "Y"){
    for (i in 1:nrow(id)){
      system(str_c("bcftools view --force-samples -S ",
                   "./03_INPUT_SampleName.txt ",
                   id$V1[i],".vcf.gz > ",
                   id$V1[i],".vcf"))
    }
    system("mv ./*vcf ./01_out_byGeneID/")
    system("rm -rf ./*vcf.gz")
    system(str_c("tar -czvf ",format(Sys.Date(), "%Y_%m_%d"),"_",JOB,"_ExportFrom_",db_name,
                 "_LOTSample_Filter_ByGeneID",".tar.gz ./01_out_byGeneID/* ./Tips.pdf"))
  }else{
    system("mv ./T* ./01_out_byGeneID/")
    system(str_c("tar -czvf ",format(Sys.Date(), "%Y_%m_%d"),"_",JOB,"_ExportFrom_",db_name,
                 "_AllSample_Filter_ByGeneID",".tar.gz ./01_out_byGeneID/* ./Tips.pdf"))
  }
  cli::cli_text("INFO   任务运行结束,请及时下载结果文件,下次运行前将清空结果文件")
}

根据物理位置提取变异信息

根据指定的物理区间判断染色体的和起止位置,并结果VCF文件筛选指定区间内的变异数据,采用bcftools的 filter功能进行实现,提取完成后进行打包压缩。

if (OPT == "2"){
  cli::cli_text("INFO   待提取物理区间如下,正在提取中......")
  region <- read.table("./02_INPUT_Postion.txt",header = F)
  for (i in 1:nrow(region)){
    print(str_c("Index: ",i,"   Region: ",region$V1[i],"   Info: ",region$V2[i]))
    system(str_c("bcftools filter ",db_file," --regions ",region$V1[i]," > ",region$V1[i],"_",region$V2[i],".vcf"))
    system("mv Chr* ./02_out_byPostion/")
    system("rename : _ ./02_out_byPostion/*")
    system("rename - _ ./02_out_byPostion/*")
  }
  cli::cli_text("INFO   提取完成,对结果进行打包压缩")
  system(str_c("tar -czvf ",format(Sys.Date(), "%Y_%m_%d"),"_",JOB,"_ExportFrom_",db_name,
               "_AllSample_Filter_ByPositin",".tar.gz ./02_out_byPostion/* ./Tips.pdf"))
  cli::cli_text("INFO   任务运行结束,请及时下载结果文件,下次运行前将清空结果文件")
}

今天分享的内容就到这里,还有两个功能正在开发中,之后有时间再分享关于提取指定位点名称和结构注释的方法,目前这项工具还未完成,如需抢先体验请后台联系,后续将开源至Github,欢迎转发分享。


参考资料:

https://mp.weixin.qq.com/s/DdXyqiW7c7lCp4103flQZQ

END

© 素材来源于网络,侵权请联系后台删除

笔记合集,点击直达

相关文章
|
3月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
4月前
|
存储 数据采集 数据处理
R语言数据变换:使用tidyr包进行高效数据整形的探索
【8月更文挑战第29天】`tidyr`包为R语言的数据整形提供了强大的工具。通过`pivot_longer()`、`pivot_wider()`、`separate()`和`unite()`等函数,我们可以轻松地将数据从一种格式转换为另一种格式,以满足不同的分析需求。掌握这些函数的使用,将大大提高我们处理和分析数据的效率。
|
2月前
|
数据挖掘 C语言 C++
R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。
【10月更文挑战第21天】时间序列分析是一种重要的数据分析方法,广泛应用于经济学、金融学、气象学、生态学等领域。R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。本文将介绍使用R语言进行时间序列分析的基本概念、方法和实例,帮助读者掌握R语言在时间序列分析中的应用。
48 3
|
2月前
|
机器学习/深度学习 并行计算 数据挖掘
R语言是一种强大的统计分析工具,广泛应用于数据分析和机器学习领域
【10月更文挑战第21天】R语言是一种强大的统计分析工具,广泛应用于数据分析和机器学习领域。本文将介绍R语言中的一些高级编程技巧,包括函数式编程、向量化运算、字符串处理、循环和条件语句、异常处理和性能优化等方面,以帮助读者更好地掌握R语言的编程技巧,提高数据分析的效率。
52 2
|
2月前
|
数据采集 数据可视化 数据挖掘
R语言与Python:比较两种数据分析工具
R语言和Python是目前最流行的两种数据分析工具。本文将对这两种工具进行比较,包括它们的历史、特点、应用场景、社区支持、学习资源、性能等方面,以帮助读者更好地了解和选择适合自己的数据分析工具。
36 2
|
3月前
R语言基于表格文件的数据绘制具有多个系列的柱状图与直方图
【9月更文挑战第9天】在R语言中,利用`ggplot2`包可绘制多系列柱状图与直方图。首先读取数据文件`data.csv`,加载`ggplot2`包后,使用`ggplot`函数指定轴与填充颜色,并通过`geom_bar`或`geom_histogram`绘图。参数如`stat`, `position`, `alpha`等可根据需要调整,实现不同系列的图表展示。
|
3月前
|
数据采集 数据可视化 数据挖掘
R语言在金融数据分析中的深度应用:探索数据背后的市场智慧
【9月更文挑战第1天】R语言在金融数据分析中展现出了强大的功能和广泛的应用前景。通过丰富的数据处理函数、强大的统计分析功能和优秀的可视化效果,R语言能够帮助金融机构深入挖掘数据价值,洞察市场动态。未来,随着金融数据的不断积累和技术的不断进步,R语言在金融数据分析中的应用将更加广泛和深入。
|
4月前
|
数据采集 机器学习/深度学习 数据挖掘
R语言数据清洗:高效处理缺失值与重复数据的策略
【8月更文挑战第29天】处理缺失值和重复数据是数据清洗中的基础而重要的步骤。在R语言中,我们拥有多种工具和方法来有效地应对这些问题。通过识别、删除或插补缺失值,以及删除重复数据,我们可以提高数据集的质量和可靠性,为后续的数据分析和建模工作打下坚实的基础。 需要注意的是,处理缺失值和重复数据时,我们应根据实际情况和数据特性选择合适的方法,并在处理过程中保持谨慎,以避免引入新的偏差或错误。
|
4月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
4月前
|
数据处理
R语言数据合并:掌握`merge`与`dplyr`中`join`的巧妙技巧
【8月更文挑战第29天】如果你已经在使用`dplyr`进行数据处理,那么推荐使用`dplyr::join`进行数据合并,因为它与`dplyr`的其他函数(如`filter()`、`select()`、`mutate()`等)无缝集成,能够提供更加流畅和一致的数据处理体验。如果你的代码中尚未使用`dplyr`,但想要尝试,那么`dplyr::join`将是一个很好的起点。