从FASTA文件中批量提取指定序列【Python脚本】
文章目录
- 前言
- 一:读取含特定字符的序列并输出
- 演示
- 二:读到某一个字符之前的全部输出
- 使用方法
- 三:输出前n条序列
- 使用方法
- 总结
前言
背景:学测序流程的时候,做到mapping的时牛的基因组有两个多G,因为是在个人PC上初步学习,建立index实在太慢了,而且临时也没有现成的index。于是决定只挑基因组前十条染色体拿来练习(所以需要从基因组文件里选取序列,尝试自己用python写脚本处理)。自己的python学习之前因为脑子实在不聪明也就刚学到字典部分,借鉴了网上的一些内容整理了一下,写了几种情况关于利用python脚本提取fasta基因序列。
一:读取含特定字符的序列并输出
比如现在有一个从NCBI上下载的fasta格式的文件,里面有很多条序列,我现在需要从中选取序列名中含有“XXX”的这些序列,那就可以用我在网上找到并改进的这个脚本
举个例子(用法中的三个参数):
[fasta_file]
下载的数据源fasta文件
[namelist_file]
把我们要找的名字另建一个txt文件,一行一个
[outfile_name]
输出的结果文件,自己命名
# -*- coding: utf-8 -*-import sys
import time
start=time.clock()def usage():print('Usage: python script.py [fasta_file] [namelist_file] [outfile_name]')def main():outf = open(sys.argv[3],'w')dict = {}with open(sys.argv[1], 'r') as fastaf:for line in fastaf:if line.startswith('>'):name = linedict[name] = ''else:dict[name] += line.replace('\n','') #读取整个fasta文件构成字典with open(sys.argv[2],'r') as listf:for row in listf:row = row.strip()for key in dict.keys(): #选取包含所需名称的基因名和序列if row in key:outf.write(key)outf.write(dict[key] + '\n')outf.close()try:main()
except IndexError:usage()end=time.clock()
print 'Running time: %s Seconds'%(end-start) #运行时间
这个脚本不能用来处理大文件!!!
演示
[fasta_file]
这是需要被筛选读取的源fasta文件示例
[namelist_file]
这是事先准备的包含要选出来的序列名的txt文件(每个名字独立一行),比如这里我要选取在上面那个test文件中(也就是[fasta_file])基因名包含05
/01
/04
的所有序列
然后,运行该脚本,我是在Linux下(也可以直接在Windows下的python环境,我这里Linux用的python 2.7)
被读取的[fasta_file]也就是我测试的test.fasta大小为4.5M(这个脚本不能用来处理大文件!!!),这里实测耗时0.047秒,然后会生成下面的结果文件
[outfile_name]
这是我输出的结果(即[outfile_name])
二:读到某一个字符之前的全部输出
背景:找到上面那个方法之后用了一下,发现对于这种几个G的大文件根本行不通,先要读取变成字典,再从字典里去找,太耗时了(我花了几个小时也没运行完得到结果)。所以说不能被别人的脚本限制了思维,想了想我反正只要前十条,那读到第十一条染色体是时候停止读写不就行了?于是自己重写了个
同样,这个也适用于**“输出某个字符之前的所有内容”**这一场景,而在这个例子中:
'result.fasta'
就是我最后的结果文件
'fasta.fasta'
这里就是我的基因组文件,即源数据fasta文件
上面两个参数如果和脚本不在同一个目录下则可以自己在脚本中添加地址,文件名自行替换
脚本中的'chromosome 11'
可以根据自己的需要替换
# -*- coding: utf-8 -*-import sys
import time
start=time.clock()outf=open('result.fasta','w')
inf=open('fasta.fasta', 'r')
for line in inf:if 'chromosome 11' not in line: #没遇到11号染色体就一直读并输出outf.write(line)else: break #读到11号就跳出循环
outf.close()
inf.close()end=time.clock()
print 'Running time: %s Seconds'%(end-start)
使用方法
直接把该脚本放到需要处理的fasta文件目录下运行就可以(比如我直接IDLE),也可以手动修改成 一:读取含特定字符的序列并输出 的版本,同样结果也生成在该目录下
这个脚本可以用来处理大文件,因为是读一行写一行。
三:输出前n条序列
背景:朋友给了我这个解法我才知道基因组文件是多少条染色体就有多少条序列,比如牛是30条染色体,基因组fna文件中就是30条序列,从1到30按顺序排列,那我要前十条染色体的序列的话,逐条读取输出,并在第11条处终止就行了,如下
if chr_num<=10:
这条条件语句可以根据需要序列数自行更改
fa_file = open("fasta.fasta","r")
out_file=open("result.fasta",'w')chr_num=0
for each_line in fa_file:each_line = each_line.strip()if each_line.startswith(">"):chr_num +=1 #利用序列名判断if chr_num<=10: out_file.write(each_line+"\n") #输出前10条else:breakfa_file.close()
out_file.close()
使用方法
同前第二种,不赘述了
可以处理大fasta文件
总结
上述这三种都可以根据需要自行调整那几个关键参数,第一种是我借鉴的其他博主的,但发现不适用于大的fasta序列文件,后两种自己写的,实测对将近3G的基因组也就一分钟以内。但第一种可以挑选多个条件(不同基因名)的序列,而后两种的限定条件比较窄,特定字符或者特定数目。
边学边练习(实际上我就刚刚学完了字典……所以这三个脚本基本都是一定看得懂的吧),肯定有很多不对的地方,欢迎交流,请多指教~
参考资料:
[根据ID从FASTA文件中批量提取序列【Python脚本】
从FASTA文件中批量提取指定序列【Python脚本】相关推荐
- 根据ID从FASTA文件中批量提取序列【Python】
根据ID从FASTA文件中批量提取序列[Python] 生信问题记录 我的需求 input: FASTA文件,含六千余个蛋白序列.命名为FA.fasta txt文件,经过interpro注释后,筛选出 ...
- 从勘界图批量提取宗地红线到shapefile工具,从CAD图中批量提取指定类型图形到shp数据实现方法。
在实际工作中可能会遇到需要从勘界图dwg中提红线的工作,当面对成百上千的勘界图时,人工逐一提取将非常繁琐耗时.下面介绍一个利用FME从勘界图批量提取宗地红线的方法. 关注薇信工众号:"GIS ...
- python从文件中读取数据_使用Python脚本从文件读取数据代码实例
这篇文章主要介绍了使用Python脚本从文件读取数据代码实例,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 最近自学Python的进度比较慢,工作之 ...
- (连载1.1)从招股说明书pdf文件中批量提取财务报表数据
目录标题 背景说明 阅读代码时注意点 源代码 背景说明 本文选取的是当前日期上交所科创板所有的上市公司样本. 用八爪鱼从上交所公告页面爬取公告下载链接,使用迅雷批量下载. 阅读代码时注意点 流程思路: ...
- 利用NCO或者CDO从nc文件中批量提取数据
nco提取数据方法: ncrcat -v bsf MMEAN*.nc -o bsf.nc cdo提取数据方法: cdo select,name=PRECT NBF1850_f19_tn11*.nc P ...
- WordSR 在多个Word文件中批量查找替换
这阵子在看一些技术文档,都是 word 格式的,需要在多个Word文件中批量查找指定的内容,找不到免费的合适的软件,顺手开发了这个工具软件,下载地址 WordSR v0.2,下面是一些版权信息和软件介 ...
- linux提取fasta文件的id,从大的fasta文件中提取特定的fasta序列
我想使用以下脚本从大的fasta文件中提取特定的fasta序列,但输出为空.从大的fasta文件中提取特定的fasta序列 transcripts.txt文件包含我想从assembly.fasta到s ...
- java批量提取文件夹名称_bat 批量提取指定目录下的文件名
bat 批量提取指定目录下的文件名 下面是批量获取指定目录下的文件名的核心代码 @echo off echo text input set input= set /p input=: echo %in ...
- fasta文件中序列的排序
同样的名为read_1.fa 的fasta文件,里面有若干序列,如: >@r1 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAA ...
最新文章
- 路径前面加/和不加/
- 吴恩达深度学习课程deeplearning.ai课程作业:Class 2 Week 2 Optimization methods
- MFC_选择目录对话框_选择文件对话框_指定目录遍历文件
- 分布式数据库的最新发展情况
- NGINX优化之路(一)
- observable java_RxJava之Observables类型理解
- Node.js的安装下载和运行JS代码和常用命令和按键
- 如何批量将 Excel 转换为 jpeg、png、bmp 图片
- VBA从工作表另存为工作簿
- 计蒜客——整数转换成罗马数字
- html 播放微信amr音频文件,微信amr文件打开的方法
- chrome去广告插件 去掉百度热搜
- 软件配置---重装系统---品牌电脑重启快捷键表
- 如何实现有效的项目进度控制
- 5G的速度到底能有多快
- python 爬取企业注册信息_读书笔记(十)——python简单爬取企查查网企业信息,并以excel格式存储...
- word加密文档忘记密码了如何打开
- 计算机科学数电吗,“不插电的计算机科学”, 你试过吗?
- 【相机】(2)——WebView中打开相机、文件选择器的问题和解决方法
- Apache APISIX 玩转 Tongsuo 国密插件
热门文章
- ExcludeClipRect和无闪烁图像
- 海思HI35xx平台软件开发快速入门之H264解码实例
- 【正解】LaTex插入空白页
- Android项目开发:简易计步器
- python设计一个学生类姓名年龄成绩_C# 编写学生类Student,包含学生姓名,成绩,设计一个友员函数sortDegree(),将学生成绩按大到小排序。...
- springboot 集成kafka 实现多个customer不同group
- css样式的格式是什么,css的语法格式是什么
- 关于固态硬盘的一些总结
- PCI与PCIe学习一——硬件篇
- Winodws 7 专业番茄花园版 v 1.0