算法丨根据基因型VCF文件自动识别变异位点并生成序列fasta文件,基于R语言tidyverse

简介: 算法丨根据基因型VCF文件自动识别变异位点并生成序列fasta文件,基于R语言tidyverse

根据VCF文件生成序列fasta文件

首先提出一个问题:假如有一个基因型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表示某一基因型,参数r0r1rx分别表示三种不同的结果。具体的实现逻辑是,若基因型为“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

© 素材来源于网络,侵权请联系后台删除

笔记合集,点击直达

相关文章
|
12天前
|
机器学习/深度学习 存储 算法
解锁文件共享软件背后基于 Python 的二叉搜索树算法密码
文件共享软件在数字化时代扮演着连接全球用户、促进知识与数据交流的重要角色。二叉搜索树作为一种高效的数据结构,通过有序存储和快速检索文件,极大提升了文件共享平台的性能。它依据文件名或时间戳等关键属性排序,支持高效插入、删除和查找操作,显著优化用户体验。本文还展示了用Python实现的简单二叉搜索树代码,帮助理解其工作原理,并展望了该算法在分布式计算和机器学习领域的未来应用前景。
|
5天前
|
算法 数据安全/隐私保护 计算机视觉
基于FPGA的图像双线性插值算法verilog实现,包括tb测试文件和MATLAB辅助验证
本项目展示了256×256图像通过双线性插值放大至512×512的效果,无水印展示。使用Matlab 2022a和Vivado 2019.2开发,提供完整代码及详细中文注释、操作视频。核心程序实现图像缩放,并在Matlab中验证效果。双线性插值算法通过FPGA高效实现图像缩放,确保质量。
|
14天前
|
存储 算法 Java
解锁“分享文件”高效密码:探秘 Java 二叉搜索树算法
在信息爆炸的时代,文件分享至关重要。二叉搜索树(BST)以其高效的查找性能,为文件分享优化提供了新路径。本文聚焦Java环境下BST的应用,介绍其基础结构、实现示例及进阶优化。BST通过有序节点快速定位文件,结合自平衡树、多线程和权限管理,大幅提升文件分享效率与安全性。代码示例展示了文件插入与查找的基本操作,适用于大规模并发场景,确保分享过程流畅高效。掌握BST算法,助力文件分享创新发展。
|
1月前
|
存储 算法 Serverless
剖析文件共享工具背后的Python哈希表算法奥秘
在数字化时代,文件共享工具不可或缺。哈希表算法通过将文件名或哈希值映射到存储位置,实现快速检索与高效管理。Python中的哈希表可用于创建简易文件索引,支持快速插入和查找文件路径。哈希表不仅提升了文件定位速度,还优化了存储管理和多节点数据一致性,确保文件共享工具高效运行,满足多用户并发需求,推动文件共享领域向更高效、便捷的方向发展。
|
1月前
|
存储 算法 安全
基于哈希表的文件共享平台 C++ 算法实现与分析
在数字化时代,文件共享平台不可或缺。本文探讨哈希表在文件共享中的应用,包括原理、优势及C++实现。哈希表通过键值对快速访问文件元数据(如文件名、大小、位置等),查找时间复杂度为O(1),显著提升查找速度和用户体验。代码示例展示了文件上传和搜索功能,实际应用中需解决哈希冲突、动态扩容和线程安全等问题,以优化性能。
|
9月前
|
机器学习/深度学习 存储 算法
用kNN算法诊断乳腺癌--基于R语言
用kNN算法诊断乳腺癌--基于R语言
|
5月前
R语言基于表格文件的数据绘制具有多个系列的柱状图与直方图
【9月更文挑战第9天】在R语言中,利用`ggplot2`包可绘制多系列柱状图与直方图。首先读取数据文件`data.csv`,加载`ggplot2`包后,使用`ggplot`函数指定轴与填充颜色,并通过`geom_bar`或`geom_histogram`绘图。参数如`stat`, `position`, `alpha`等可根据需要调整,实现不同系列的图表展示。
|
6月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
7月前
|
算法
Ngnix02 --- Ngnix的功能特性及常见功能,Ngnix常用的功能模块,有不同算法,根据不同算法进行转发,ip_hash、url_hash、fair,核心组成 ngnix二进制可执行文件
Ngnix02 --- Ngnix的功能特性及常见功能,Ngnix常用的功能模块,有不同算法,根据不同算法进行转发,ip_hash、url_hash、fair,核心组成 ngnix二进制可执行文件
|
8月前
|
Python
R语言遍历文件夹求取其中所有栅格文件的平均值
通过NAvalue(tif_file_all) <- -10000这句代码,将值为-10000的像元作为NoData值的像元,防止后期计算平均值时对结果加以干扰。   接下来,我们通过file.path()函数配置一下输出结果的路径——其中,结果遥感影像文件的名称就可以直接以其所对应的条带号来设置,并在条带号后添加一个_mean后缀,表明这个是平均值的结果图像;但此外,这个仅仅是文件的名字,还需要将文件名与路径拼接在一起,才可以成为完整的保存路径,因此需要用到file.path()函数。最后,将结果图像通过writeRaster()函数加以保存即可,这句代码的解释大家同样参考R语言求取大量遥感
200 0

热门文章

最新文章