R语言丨根据VCF文件自动填充对其变异位点并生成序列fa文件

根据VCF文件自动填充对其变异位点并生成序列fa文件

首先提出一个问题：
假如有一个重测序结果VCF文件，里面包含了很多个样本在几百个突变位点（snp和iad）的基因型数据，现在想根据这份原始数据，得到一个fasta序列文件，包含每个样品在这些位点的各自序列信息，应该怎么做？

解决思路与方法简介

方法一：Excel手工处理

1. 将vcf文件转成Excel表格
2. 判断每个变异位点的类型是snp或者iad
3. 如果是iad，将REF和ALT补全，保证长度一致
4. 根据基因型判断填充内容
5. 按顺序合并每个样品对应的位点信息
6. 生成fa序列文件

方法二：R语言算法处理

1. 读入vcf文件并切分成基因信息、样品列表、位点信息三部分
2. 对位点信息进行处理，枚举每个碱基，判断长度之和
3. 若总长度为2，则说明为单点突变SNP
4. 若总长度大于2，则为插入缺失突变，进一步锁定单碱基
5. 填充占位符，保证插入缺失的两端对齐
6. 根据位点类型循环填充基因型信息
7. 根据样品列表迭代生成所有样品的fa序列

主要算法和步骤

加载数据

rm(list=ls())
library(tidyverse)
library(xlsx)
df <- read.xlsx("./genevcf.xlsx",sheetIndex = 1)
num_snp <- ncol(df)-12 #这里修改表头列数
num_sam <- nrow(df)-7  #这里修改表头行数
print(str_c("SNP number:",num_snp,"    Sample number:",num_sam))

这段代码的主要功能是读取一个名为“genevcf.xlsx”的Excel文件的第一个工作表，并根据读取到的数据计算SNP位点的数量和样本的数量。

其中，通过使用library(tidyverse)和library(xlsx)导入了两个R语言包以进行数据处理和读写Excel文件。

首先，read.xlsx()函数被用来读取Excel文件，其中sheetIndex = 1表示选取第一个工作表作为读取的对象。

之后，代码从数据中获取SNP位点的数量，这里假设该Excel文件的表头有12列，将除去这12列的剩余列数赋值给变量num_snp。

同理，假设该Excel文件的表头有7行，则将除去这7行的剩余行数赋值给变量num_sam。最后，通过print()函数和str_c()函数将计算得到的SNP位点数和样本数打印出来。

拆分与整理数据

df_snp <- cbind(df[1:7,1],df[1:7,13:ncol(df)]) #snp信息
df_snp <- t(df_snp)
colnames(df_snp) <- df_snp[1,]
df_snp <- df_snp[-1,-7] 
df_snp <- as.data.frame(df_snp)
df_snp$len <- nchar(str_c(df_snp$REF,df_snp$ALT))
rownames(df_snp) <- 1:nrow(df_snp)
df_snp$type_0 <- NA # 纯合0/0
df_snp$type_1 <- NA # 纯合1/1
df_snp$type_x <- NA # 杂合1/0和0/1

1. 将df_snp矩阵进行转置，即将行列互换。
2. 将矩阵的第一行的值作为列名。
3. 删除第一行第七列的值。
4. 将矩阵转化为数据框。
5. 将REF和ALT两列的字符串拼接在一起，计算字符串长度，将得到的新列命名为len。
6. 将数据框的行名重置为1到行数。

计算基因型组合方式

for (i in 1:nrow(df_snp)){
      if (df_snp$len[i]==2){ #单个snp
            df_snp$type_0[i] <- str_c(df_snp$REF[i],df_snp$REF[i])
            df_snp$type_1[i] <- str_c(df_snp$ALT[i],df_snp$ALT[i])
            df_snp$type_x[i] <- str_c(df_snp$REF[i],df_snp$ALT[i])
      }else{
            n_REF <- nchar(df_snp$REF[i])
            n_ALT <- nchar(df_snp$ALT[i])
            if (n_REF == 1){ # REF为单个碱基，ALT为多碱基
                  df_snp$type_0[i] <- str_c(df_snp$REF[i],paste(rep("-",n_ALT),collapse = ""))
                  df_snp$type_1[i] <- str_c("-",df_snp$ALT[i])
                  df_snp$type_x[i] <- str_c(df_snp$REF[i],df_snp$ALT[i])
            }else{ # ALT为单碱基
                  df_snp$type_1[i] <- str_c(df_snp$ALT[i],paste(rep("-",n_REF),collapse = ""))
                  df_snp$type_0[i] <- str_c("-",df_snp$REF[i])
                  df_snp$type_x[i] <- str_c(df_snp$REF[i],df_snp$ALT[i])
            }
      }
}

1. 对于数据框df_snp中每一行：
2. 如果该行的len等于2，说明该行代表单个SNP，则将该行的REF和ALT两列的值拼接得到三个新的列：type_0，type_1和type_x。其中type_0和type_1分别代表两个等位基因，type_x代表两个等位基因的组合。
3. 如果该行的REF为单个碱基，ALT为多碱基，则将该行的REF和ALT两列的值拼接得到三个新的列：type_0，type_1和type_x。其中type_0代表REF基因型，type_1代表ALT基因型，type_x代表两个等位基因的组合。在得到type_0时，会在REF基因前面加上与ALT基因数量相等的“-”符号。
4. 如果该行的ALT为单个碱基，则将该行的REF和ALT两列的值拼接得到三个新的列：type_0，type_1和type_x。其中type_1代表ALT基因型，type_0代表REF基因型，type_x代表两个等位基因的组合。在得到type_1时，会在ALT基因前面加上与REF基因数量相等的“-”符号。

变异类型转换函数

df_gene <- df[8:nrow(df),13:ncol(df)] 
rownames(df_gene) <- 1:nrow(df_gene)
colnames(df_gene) <- str_c("snp",1:ncol(df_gene))
function_snp_convert <- function(n,r0,r1,rx){
      if (n == "0|0"){
            return(r0)
      }else{
            if (n == "1|1"){
                  return(r1)
            }else{
                  return(rx)
            }
      }
}

首先，代码从数据框df中提取出基因型信息存储到df_gene中，同时对行列名进行了重新命名，行名从1开始递增，列名以snp 加上相应的列数编号。接下来定义function_snp_convert函数，这个函数被设计成可以计算不同基因型的组合方式，并分别返回对应的结果值。其中，参数n表示某一基因型，参数r0、r1和rx分别表示三种不同的结果。具体的实现逻辑是，若基因型为“0|0”，则返回r0；若基因型为“1|1”，则返回r1；否则如果不是以上两种情况，则返回rx。

遍历替换基因数据

for (x in 1:nrow(df_gene)){
      for (y in 1:ncol(df_gene)){
           df_gene[x,y] <- function_snp_convert(df_gene[x,y],
                                 df_snp$type_0[y],
                                 df_snp$type_1[y],
                                 df_snp$type_x[y])
      }
}

使用了两个循环，逐行逐列地读取df_gene数据框中的元素，并调用function_snp_convert函数将读取的元素代入函数相应的参数中进行计算，将计算结果存储回df_gene数据框中对应的位置，以此实现了对df_gene数据框中所有元素根据不同基因型进行计算的功能。

for循环的x和y分别表示df_gene数据框中的行下标和列下标，每次循环时调用function_snp_convert函数，并传入相应的参数，将计算结果存储回df_gene数据框的对应元素中。其中，df_gene[x, y]实际上就是df_gene数据框中第x行第y列位置的元素，而df_snptype_0[y]、df_snptype0[y]、dfsptype_1[y]和df_snp$type_x[y]分别表示df_snp数据框中第y列三个不同的基因型组合方式。

生成fasta序列文件

df_name <- df[8:nrow(df),1] 
fasta <- c()
for (index in 1:length(df_name)){
      print(df_name[index])
      fasta <- c(fasta,str_c(">",df_name[index]))
      fasta <- c(fasta,str_c(df_gene[index,],collapse = ""))
      print(str_c(df_gene[index,],collapse = ""))
}
rownames(df_gene) <- df_name
colnames(df_gene) <- df_snp$ID
write.csv(df_gene,"ID_SNP_ATCG.csv",quote = F)
write.table(fasta,"all.fa",quote = F,row.names = F,col.names = F)

在这段代码中，首先提取了数据框 df 第8到最后一行的第1列数据，作为样品名称列表 df_name。

然后初始化了一个空向量 fasta。接下来的 for 循环会迭代样品名称列表 df_name 的所有元素。每次循环，会打印该元素名称，并将 “>” 和该元素名称拼接成的字符串添加至 fasta 中。

然后将 df_gene 的第 index 行拼接为字符串后加入 fasta 中。最后，将 df_name 赋值为 df_gene 数据框的行名，df_snp 的 ID 列作为 df_gene 的列名，该代码的作用是生成一个 fasta 格式的序列文件。

本文由mdnice多平台发布