根据VCF文件生成序列fasta文件
首先提出一个问题：假如有一个基因型VCF文件，里面包含了很多个样本在多个突变位点（snp和iad）的基因型数据，现在想根据这份原始数据，得到一个fasta序列文件，包含每个样品在这些位点的各自对应的序列信息，应该怎么做？
解决思路与方法简介
方法一：Excel手工处理
- 将vcf文件转成Excel表格
- 判断每个变异位点的类型是snp或者iad
- 如果是iad，将REF和ALT补全，保证长度一致
- 根据基因型判断填充内容
- 按顺序合并每个样品对应的位点信息
- 生成fa序列文件
方法二：R语言算法处理
- 读入vcf文件并切分成基因信息、样品列表、位点信息三部分
- 对位点信息进行处理，枚举每个碱基，判断长度之和
- 若总长度为2，则说明为单点突变SNP
- 若总长度大于2，则为插入缺失突变，进一步锁定单碱基
- 填充占位符，保证插入缺失的两端对齐
- 根据位点类型循环填充基因型信息
- 根据样品列表迭代生成所有样品的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
- 将df_snp矩阵进行转置，即将行列互换。
- 将矩阵的第一行的值作为列名。
- 删除第一行第七列的值。
- 将矩阵转化为数据框。
- 将REF和ALT两列的字符串拼接在一起，计算字符串长度，将得到的新列命名为len。
- 将数据框的行名重置为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]) } } }
- 对于数据框df_snp中每一行：
- 如果该行的len等于2，说明该行代表单个SNP，则将该行的REF和ALT两列的值拼接得到三个新的列：type_0，type_1和type_x。其中type_0和type_1分别代表两个等位基因，type_x代表两个等位基因的组合。
- 如果该行的REF为单个碱基，ALT为多碱基，则将该行的REF和ALT两列的值拼接得到三个新的列：type_0，type_1和type_x。其中type_0代表REF基因型，type_1代表ALT基因型，type_x代表两个等位基因的组合。在得到type_0时，会在REF基因前面加上与ALT基因数量相等的“-”符号。
- 如果该行的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 格式的序列文件。
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
