算法：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 多平台发布
