利用R语言实现vcf转txt算法,基于vcfR和tidyverse

news/2024/10/30 11:30:14/

算法:vcf转txt并自动规范化

引言

vcf文件是存放基因变异信息的一种方式,本文提供一种算法,用于读取vcf文件并转换等位基因展示方法、替换染色体展示格式、以及自动识别非唯一变异并进行修改,用于对变异信息进行整理。


主要步骤与设计思路

  • 读取VCF文件并分为三部分储存
  • 提取变异信息并批量替换
  • 修改染色体格式
  • SNP位点的判断与校正
  • 单点碱基差异唯一化

项目运行环境

  • centos7 linux
  • R4.2.3

具体操作步骤

加载R包与数据

library(tidyverse)
library(vcfR)
library(do)
library(R.utils)
df <- read.table(paste0("02_ordata/",job,".filter.vcf"),header = F)
vcf <- read.vcfR(paste0("02_ordata/",job,".filter.vcf.gz"))
chr_ref <- read.table("01_scripts/chr_num2str.txt",header = T)

读取VCF文件信息

fix <- vcf@fix
gt <- vcf@gt
meta <- vcf@meta

利用vcfR包读取入VCF文件后,分别提取出不同部分存放于临时变量中,以供后续使用。

批量替换变异信息

### 批量替换“|”为“/” ==================================================================
df[df == "0|0"] = "0/0"
df[df == "1|0"] = "1/0"
df[df == "0|1"] = "0/1"
df[df == "1|1"] = "1/1"
colnames(df) <- c(colnames(fix),colnames(gt))

该步骤的目的是为了将|修改为/,这是后面转hmp格式所需的条件。

替换染色体编号

###  替换染色体 =====================================================================
for (i in 1:nrow(df)){
  old_chr <- df$CHROM[i]
  for (k in 1:nrow(chr_ref)){
    if (chr_ref$chr_str[k] == old_chr){
      new_chr <- chr_ref$chr_num[k]
      df$CHROM[i] <- new_chr
    }
  }
}

利用for循环查找逐一取出染色体元素值,然后从参考信息中查找对应的正确格式,然后赋值给染色体信息,这一步中使用的chr_ref是染色体不同格式的对应信息。

参数识别与矫正

因为有插入缺失的存在,所以参考位置和实际位置的碱基并非完全唯一且差异,这将导致后面运行出错。这里提供一个算法,批量实现对SNP位点的检测与矫正。

  • snp_reverse函数
snp_reverse <- function(one,more){
  # 输入俩参,一为单二为多,返回存在于多但不与单同之值
  list_snp <- str_split(more,"")
  for (i in 1:str_length(more)){
    snp_now <- list_snp[[1]][i]
    ifelse(one==snp_now,next,return(snp_now))
  }
}

该函数输入两个参数,如“A,CATG”,首先将第二个参数分割成单个字母,然后迭代判断第一个字母是否与第二个一致,一旦出现与第一个参数不相同的值则返回该值。目的是为了让两个值长度为1且不相同。

批量处理ALT和REF位点

# 对每行的REF和ALT进行处理,将其变成不同值
for (i in 1:nrow(df)){
  ref <- df$REF[i]
  alt <- df$ALT[i]
  # 情况有三,均为单或其一为多
  if (str_length(ref) == 1){
    if (str_length(alt) == 1){
    }else{
      df$ALT[i] <- snp_reverse(ref,alt)
    }
  }else{
    if (str_length(alt) == 1){
      df$REF[i] <- snp_reverse(alt,ref)
    }else{
      print(paste0("ERROR:",df$ID[i]," this snp has more REF、ALT !"))
    }
  }
}

结果保存与输出

colnames(df)[1] <- "#CHROM"
write.table(df,paste0("03_vcf2txt/","gene_",job,".txt"),
            sep = "\t",row.names = F,col.names = T,quote = F)
print(paste0(job," Step ordata gene vcf to txt finished!"))

通过该算法能够对vcf文件进行转换,并得到规范化的txt文件,用于后续的分析。

本文由 mdnice 多平台发布


http://www.ppmy.cn/news/47344.html

相关文章

测试5年,从纯手工测试到测试开发,我是怎么拿到腾讯25koffer的?

什么都做了&#xff0c;和什么都没做其实是一样的&#xff0c;走出“瞎忙活”的安乐窝&#xff0c;才是避开弯路的最佳路径。希望我的经历能帮助到有需要的朋友。 在测试行业已经混了5个年头了&#xff0c;以前经常听到开发对我说&#xff0c;天天的点点点有意思没&#xff1f…

Java工程行业管理系统源码-专业的工程管理软件-提供一站式服务

Java版工程项目管理系统 Spring CloudSpring BootMybatisVueElementUI前后端分离 功能清单如下&#xff1a; 首页 工作台&#xff1a;待办工作、消息通知、预警信息&#xff0c;点击可进入相应的列表 项目进度图表&#xff1a;选择&#xff08;总体或单个&#xff09;项目显示1…

使用Softing edgePlug软件扩展数控机床的连接性

那些使用SINUMERIK 840D控制器来运行数控机床的制造商正面临着一个挑战——从车间提取机床性能和过程数据来进行分析。这些数据对于优化流程至关重要&#xff0c;但它们却无法通过传统方式来被获取。对此&#xff0c;制造商的应对方法是通过自定义代码来实现数据访问&#xff0…

BufferedOutputStream,BufferedInputStream是字节流,对象处理流,序列化,输入输出流,转换流

BufferedInputStream字节输入流 意思就是InputStream类及其子类都能以参数的形式放到BufferedInputStream构造器的参数 package com.hspedu.outputstream_;import java.io.*;/*** author 韩顺平* version 1.0* 演示使用BufferedOutputStream 和 BufferedInputStream使用* 使用他…

带你浅谈下Quartz的简单使用

Scheduler 每次执行&#xff0c;都会根据JobDetail创建一个新的Job实例&#xff0c;这样就可以规避并发访问的问题&#xff08;jobDetail的实例也是新的&#xff09; Quzrtz 定时任务默认都是并发执行&#xff0c;不会等待上一次任务执行完毕&#xff0c;只要间隔时间到就会执…

电子台账:生成的数据和图表导出到一个excel表中

目录 1 数据选择 1.1 选择1行数据 1.2 选择1列数据 2 图表设置 3 数据导出 为了便于进行数据分析和数据展示&#xff0c;可以把生成的汇总数据生成图表&#xff0c;然后对图表进行定制修改&#xff0c;最后把数据和图表一起导出到一个excel表中。 程序目前支持两种数据作…

开放式蓝牙耳机哪个好,分享几款舒适性高的开放式蓝牙耳机

开放式耳机的兴起是近几年来才出现的新概念&#xff0c;开放式耳机也是近几年来才开始流行起来&#xff0c;在我看来开放式耳机的兴起是科技进步的产物。随着蓝牙耳机技术和设备的发展&#xff0c;蓝牙耳机也越来越普及&#xff0c;但是也给用户带来了很多困扰。而开放式耳机就…

mulesoft MCIA 破釜沉舟备考 2023.04.18.16

mulesoft MCIA 破釜沉舟备考 2023.04.18.16 1. A mule application is required to periodically process large data set from a back-end database to Salesforce CRM using batch job scope configured properly process the higher rate of records.2. A company is design…