算法：vcf轉(zhuǎn)txt并自動(dòng)規(guī)范化

引言

vcf文件是存放基因變異信息的一種方式，本文提供一種算法，用于讀取vcf文件并轉(zhuǎn)換等位基因展示方法、替換染色體展示格式、以及自動(dòng)識(shí)別非唯一變異并進(jìn)行修改，用于對(duì)變異信息進(jìn)行整理。

主要步驟與設(shè)計(jì)思路

• 讀取VCF文件并分為三部分儲(chǔ)存
• 提取變異信息并批量替換
• 修改染色體格式
• SNP位點(diǎn)的判斷與校正
• 單點(diǎn)堿基差異唯一化

項(xiàng)目運(yùn)行環(huán)境

• centos7 linux
• R4.2.3

具體操作步驟

加載R包與數(shù)據(jù)

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文件后，分別提取出不同部分存放于臨時(shí)變量中，以供后續(xù)使用。

批量替換變異信息

### 批量替換“|”為“/” ==================================================================
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))

該步驟的目的是為了將|修改為/,這是后面轉(zhuǎn)hmp格式所需的條件。

替換染色體編號(hào)

###  替換染色體 =====================================================================
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循環(huán)查找逐一取出染色體元素值，然后從參考信息中查找對(duì)應(yīng)的正確格式，然后賦值給染色體信息，這一步中使用的chr_ref是染色體不同格式的對(duì)應(yīng)信息。

參數(shù)識(shí)別與矯正

因?yàn)橛胁迦肴笔У拇嬖冢詤⒖嘉恢煤蛯?shí)際位置的堿基并非完全唯一且差異，這將導(dǎo)致后面運(yùn)行出錯(cuò)。這里提供一個(gè)算法，批量實(shí)現(xiàn)對(duì)SNP位點(diǎn)的檢測(cè)與矯正。

• snp_reverse函數(shù)

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))
  }
}

該函數(shù)輸入兩個(gè)參數(shù)，如“A,CATG”，首先將第二個(gè)參數(shù)分割成單個(gè)字母，然后迭代判斷第一個(gè)字母是否與第二個(gè)一致，一旦出現(xiàn)與第一個(gè)參數(shù)不相同的值則返回該值。目的是為了讓兩個(gè)值長(zhǎng)度為1且不相同。

批量處理ALT和REF位點(diǎn)

# 對(duì)每行的REF和ALT進(jìn)行處理，將其變成不同值
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 !"))
    }
  }
}

結(jié)果保存與輸出

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!"))

通過該算法能夠?qū)cf文件進(jìn)行轉(zhuǎn)換，并得到規(guī)范化的txt文件，用于后續(xù)的分析。

本文由mdnice多平臺(tái)發(fā)布