fasta文件与fastq文件相互转化Python脚本

fa文件与fq文件互相转换

今天分享的内容是fasta文件与fastq文件的基本知识,以及通过Python实现两者互相转换的方法。

测序数据公司给的格式通常是fq.gz,也就是fastq文件,计算机的角度来说,生物的序列属于一种字符串,就是一堆字母,这些字母就蕴含了遗传信息。

通过序列拼接将fastq转换为fasta,通过短序列比对将fastq与fasta合并为bam,通过变异检测将bam中突变位点提取出来转换为vcf,这就是上游分析的套路。

fastq文件基本格式
alt

可以看出fq文件包含了更多的信息,比如测序质量,碱基信息等,这些是通过测序仪产生的数据。

fasta文件基本格式

alt 对比一下可以看出,fa文件主要是两部分,大于号开头的是序列的ID,下一行是序列,相比于fq文件,少了质量信息。

将fasta文件转换为fastq文件

分享一个Python脚本实现这个操作:

import sys

fa_f = sys.argv[1]
fq_f = sys.argv[2]
len_max = int(sys.argv[3])  # fq len max: 150bp

with open(fa_f, 'r'as f, open(fq_f, 'w'as pf:
    while True:
        name = f.readline().strip()
        if not name:
            break
        if name.startswith('>'):
            tmp_name = name.strip('>')
            read_id = '@'+tmp_name
        seq = f.readline().strip()
        if len(seq) <= len_max:
            read_info = '{rd_id}\n{seq}\n+\n{qual}\n'.format(rd_id=read_id, seq=seq, qual='F'*len(seq))
        else:
            read_info = ''
            cnt = int(len(seq) / len_max)
            for i in range(cnt):
                tmp_id = read_id + '_' + str(i)
                tmp_seq = seq[i*len_max:(i+1)*len_max]
                read_info += '{rd_id}\n{seq}\n+\n{qual}\n'.format(rd_id=tmp_id, seq=tmp_seq, qual='F'*len_max)
            lseq = len(seq) % len_max
            if lseq != 0:
                tmp_id = read_id + '_' + str(cnt)
                tmp_seq =seq[-lseq:]
                read_info += '{rd_id}\n{seq}\n+\n{qual}\n'.format(rd_id=tmp_id, seq=tmp_seq, qual='F'*lseq)
        pf.write(read_info)

使用的方法也很简单,把这个脚本保存为xx.py,然后运行并添加三个参数,第一个是原始fasta文件名,第二个是输出文件名,第三个参数是数字,表示每条序列的最大长度,超过该长度的序列将会被切分成多条。

原理解释

刚刚这段Python脚本的功能是将fasta格式的序列文件转换为fastq格式的序列文件,并且可以对序列进行分割,使得每条序列的长度不超过指定的最大长度。

功能:

读取输入的fasta格式的序列文件。 将fasta序列文件中的序列转换为fastq格式。 如果序列长度超过指定的最大长度(len_max),则将长序列分割成多个子序列,每个子序列长度不超过len_max。 将转换后的fastq格式的序列写入输出文件中。

原理:

通过命令行参数传入fasta格式的序列文件路径(fa_f)、要生成的fastq序列文件路径(fq_f)和最大序列长度(len_max)。 使用with open()语句打开fasta序列文件和要生成的fastq序列文件。

逐行读取fasta序列文件,每次读取两行:第一行为序列ID,第二行为序列信息。 对于每条序列,如果序列长度不超过指定的最大长度,则直接转换为fastq格式;否则,将长序列分割成多个子序列,每个子序列长度不超过len_max。

将fastq文件转换为fasta文件

同样,我们也可以使用Python将fq文件转换为fa文件:

import sys
import gzip 

fq_in = sys.argv[1]
fa_out = sys.argv[2]
reads_count = sys.argv[3]  # if set [-1], means output all reads

with gzip.open(fq_in, 'r') as f, open(fa_out, 'w') as pf:
    cnt = 0
    while True:
        rd_id = f.readline()
        if not rd_id or cnt==int(reads_count):
            break
        seq = f.readline()
        tmp = f.readline()
        qual = f.readline()
        pf.write('>'+rd_id+seq)
        cnt+=1

这段Python代码是一个简单的脚本,用于将gzip压缩的Fastq文件(.fq.gz文件)转换为普通的Fasta文件(.fa文件), 下面是代码的原理和作用:

首先,导入了sys和gzip模块,sys用于接收命令行参数,gzip用于解压缩.fq.gz文件。 从命令行参数中获取输入Fastq文件路径(fq_in)、输出Fasta文件路径(fa_out)和要输出的reads数量(reads_count)。

使用gzip.open函数打开输入的Fastq文件,以只读模式打开。使用open函数打开输出的Fasta文件,以写入模式打开。 设置一个计数器cnt,用于记录已经处理的reads数量。

进入一个无限循环,循环中读取Fastq文件中的每个reads信息: 读取reads的ID行(以'@'开头的行)作为rd_id。 读取reads的序列行作为seq。 读取reads的空行(通常为'+')作为tmp。 读取reads的质量信息行作为qual。

将reads的ID和序列信息写入输出的Fasta文件中,格式为>rd_idseq。 计数器cnt加一。 如果读取的reads数量达到指定的reads_count值,则退出循环。 循环结束后,关闭输入和输出文件。

总的来说,将压缩的Fastq文件解压缩并转换为Fasta格式,同时可以根据指定的reads数量控制输出的reads数量。代码中使用了gzip模块解压缩文件,以及文件读取和写入操作,最终实现了Fastq到Fasta的转换。

以上就是今天分享的全部内容,感谢您的阅读,如果感觉有用,欢迎收藏或者转发,您的支持是我更新的最大动力。

参考资料:
https://blog.csdn.net/sinat_32872729/article/details/117353884
https://blog.csdn.net/weixin_46128755/article/details/127947650
https://zhuanlan.zhihu.com/p/77874271

本文由 mdnice 多平台发布

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/729882.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

CVHub | 万字长文带你全面解读视觉大模型(建议收藏!)

本文来源公众号“CVHub”&#xff0c;仅用于学术分享&#xff0c;侵权删&#xff0c;干货满满。 原文链接&#xff1a;万字长文带你全面解读视觉大模型 0 导读 众所周知&#xff0c;视觉系统对于理解和推理视觉场景的组成特性至关重要。这个领域的挑战在于对象之间的复杂关系…

形参化类 ‘Result‘ 的原始使用

在编程中&#xff0c;特别是在面向对象编程&#xff08;OOP&#xff09;中&#xff0c;当我们说“形参化类”或“参数化类”&#xff0c;我们实际上是指泛型&#xff08;Generics&#xff09;的概念。泛型允许在定义类、接口和方法时使用类型参数。这样&#xff0c;你可以创建可…

统计某个字段有哪些值

统计某个字段有哪些值 POST mgb_test_alias/_search {"size": 0,"aggs": {"gameid_list": {"terms": {"field": "tag.abc","size": 100000}}},"track_total_hits":true }

前端实现点击复制内容到粘贴板功能

Document.execCommand() Document.execCommand()是操作剪贴板的传统方法&#xff0c;各种浏览器都支持。 复制时&#xff0c;先选中文本&#xff0c;然后调用document.execCommand(‘copy’)&#xff0c;选中的文本就会进入剪贴板。 当需要选中的内容是文本输入标签(input、te…

算法练习第十六天| 104.二叉树的最大深度、559. N 叉树的最大深度、111.二叉树的最小深度、222.完全二叉树的节点个数

优先掌握递归的方式 104.二叉树的最大深度 /*** Definition for a binary tree node.* public class TreeNode {* int val;* TreeNode left;* TreeNode right;* TreeNode() {}* TreeNode(int val) { this.val val; }* TreeNode(int val, TreeNode le…

大数据开发(Spark面试真题-卷一)

大数据开发&#xff08;Spark面试真题&#xff09; 1、什么是Spark Streaming&#xff1f;简要描述其工作原理。2、什么是Spark内存管理机制&#xff1f;请解释其中的主要概念&#xff0c;并说明其作用。3、请解释一下Spark中的shuffle是什么&#xff0c;以及为什么shuffle操作…

京东获得JD商品详情 API 返回值说明

jd.item_get 获取API测试 公共参数 名称类型必须描述keyString是调用key&#xff08;必须以GET方式拼接在URL中&#xff09;secretString是调用密钥api_nameString是API接口名称&#xff08;包括在请求地址中&#xff09;[item_search,item_get,item_search_shop等]cacheStrin…

99.qt qml-单例程序实现

在之前讲过: 58.qt quick-qml系统托盘实现https://nuoqian.blog.csdn.net/article/details/121855993 由于,该示例只是简单讲解了系统托盘实现,并没有实现单例程序,所以多次打开后就会出现多个exe出现的可能,本章出一章QML单例程序实现, 多次打开始终只显示出第一个打开…

DiT:Scalable Diffusion Models with Transformers

TOC 1 前言2 方法和代码 1 前言 该论文发表之前&#xff0c;市面上几乎都是用卷积网络作为实际意义上的&#xff08;de-facto&#xff09;backbone。于是一个想法就来了&#xff1a;为啥不用transformer作为backbone呢&#xff1f; 文章说本论文的意义就在于揭示模型选择对于…

AI为什么需要GPU

人工智能&#xff08;AI&#xff09;需要GPU&#xff08;图形处理单元&#xff09;主要是因为深度学习模型的训练和推理过程需要大量的计算资源。GPU相比于传统的中央处理单元&#xff08;CPU&#xff09;在并行计算方面具有明显的优势&#xff0c;能够更有效地处理大规模的数据…

用python写一个自动进程守护,带UI

功能是指定程序关闭后自动重启&#xff0c;并点击1作为启动 原来的想法是群成员说的某软件打包后&#xff0c;软件进程被杀后&#xff0c;界面白屏。所以写了个计算器重启demo进行进程守护 import subprocess import time import pyautogui import psutil #用计算器做演示。 d…

WiFi模块助力敏捷办公:现代办公室的关键角色

随着信息技术的飞速发展&#xff0c;现代办公室正经历着一场数字化和智能化的变革。在这一变革过程中&#xff0c;WiFi模块作为无线通信技术的核心组成部分&#xff0c;扮演着关键的角色&#xff0c;为敏捷办公提供了强大的支持。本文将深入探讨WiFi模块在现代办公室中的关键角…

Spring Boot工作原理

Spring Boot Spring Boot 基于 Spring 开发&#xff0c;Spirng Boot 本身并不提供 Spring 框架的核心特性以及扩展功能&#xff0c;只是用于快速、敏捷地开发新一代基于 Spring 框架的应用程序。也就是说&#xff0c;它并不是用来替代 Spring 的解决方案&#xff0c;而是和 Spr…

安康杯安全知识竞赛上的讲话稿

各位领导、同志们&#xff1a; 经过近半个月时间的准备&#xff0c;南五十家子镇平泉首届安康杯安全生产知识竞赛初赛在今天圆满落下帏幕&#xff0c;经过紧张激烈的角逐&#xff0c; 代表队、 代表队和 代表队分别获得本次竞赛的第一、二、三名让我们以热烈的掌声表示祝…

指针的进阶小tips

前情提要 指针是变量&#xff0c;存地址&#xff08;唯一&#xff09;指针4/8个字节&#xff08;32/64位平台&#xff09;指针有类型&#xff0c;其类型决定指针整数的步长&#xff0c;指针解引用操作的时候的权限。指针的运算 字符指针 int main() {char am;char *pc&c…

使用插件vue-seamless-scroll 完成内容持续动态

1、安装插件 npm install vue-seamless-scroll --save 2、项目中引入 //单独引入import vueSeamlessScroll from vue-seamless-scrollexport default {components: { vueSeamlessScroll},}//或者在main.js引入import scroll from vue-seamless-scrollVue.use(scroll)3、页面使…

SRS服务器ffmpeg 推流rtmp超时中断

ffmpeg错误显示 failed to update header with correct duration failed to update header with correct filesize. Error writing trailer of rtmp://----- broken pipe SRS日志错误显示 serve error code2056 kickoffforidle : service cycle : rtmp stream service: timeou…

【 React 】super()和super(props)有什么区别

相关文章 react 中的 super super(props) 1. ES6类 在ES6中&#xff0c;通过extends关键字实现类的继承&#xff0c;方式如下&#xff1a; class sup{constructor(name){this.namename;}printName(){console.log(this.name)} } class sub extends sup{constructor(name,age)…

基于Pytorch搭建分布式训练环境

Pytorch系列 文章目录 Pytorch系列前言一、DDP是什么二、DPP原理terms、nodes 和 ranks等相关术语解读DDP 的局限性为什么要选择 DDP 而不是 DP代码演示1. 在一个单 GPU 的 Node 上进行训练&#xff08;baseline&#xff09;2. 在一个多 GPU 的 Node 上进行训练临门一脚&#x…

【深度学习笔记】稠密连接网络(DenseNet)

注&#xff1a;本文为《动手学深度学习》开源内容&#xff0c;部分标注了个人理解&#xff0c;仅为个人学习记录&#xff0c;无抄袭搬运意图 5.12 稠密连接网络&#xff08;DenseNet&#xff09; ResNet中的跨层连接设计引申出了数个后续工作。本节我们介绍其中的一个&#xf…