zhuyi8715的个人博客分享 http://blog.sciencenet.cn/u/zhuyi8715

博文

简单的Python脚本提取对应位置基因序列(fasta文件)

已有 25283 次阅读 2014-6-8 10:09 |系统分类:科研笔记| Python, 提取, 基因序列, fasta

    最近,用Python脚本提取,在基因号已知,位置已知条件下,相对应位置的基因序列时发现,这样很简单但是很实用的脚本,在网上却比较难找。而且,能被找到的脚本,相对于具有初级编程能力的人而言,有点难。本人写了相对于初学者同样很简单脚本分享给大家。

   首先,我将fa文件处理为单行(嫌麻烦,没有写成scaffold_x一行,序列一行的样子,如图三),将下面的序列处理(图一):

(补充)经过:

import re
fr=open(r'F:desktopcorrelxxx.fa','r')
fw=open(r'F:desktopcorrelxxx_use.fa','w')
line=fr.read()
r=line.replace('n','')
s=re.sub('>','n>',r)
fw.write(s)
fr.close()
fw.close()

得到(图二):


当然你如果不嫌麻烦也可以处理成(图三):


假设我含有位置信息源文件(图四):

第一列为基因号,最后一列为基因在fa文件中的位置信息;


本人采用图二的形式,具体脚本(脚本一);

#author:Wang Binzhong

# -*- coding:utf-8 -*-
fr=open(r'F:desktopCX.txt','r')#读取含有位置信息的文件
fa=open(r'F:desktopxxxx.fa','r')#读取处理好的基因序列文件
fw_1=open(r'F:desktopfa_3.txt','w')#写入
line_cr=fr.readlines()
line_fa=fa.readlines()
for eachline in line_cr:
   sp=eachline.strip().split('t')
   title_1=eachline.find('scaffold')
   start_1=eachline.find(':',title_1)+1
   end_1=eachline.find('-',start_1)
   d_1=eachline[title_1:start_1-1].strip()#scaffold名称
   d_2=eachline[start_1:end_1].strip()#首位的位置
   d_3=eachline[end_1+1:].strip()#末尾的位置
   for each_seq in line_fa:
       if d_1 == each_seq[:int(len(d_1))+5].strip('ATGC'):#如果对应的名称在行中,就可以用以下的规则写入文本
           fw_1.write(sp[0]+'t'+each_seq[len(d_1)+int(d_2):len(d_1)+int(d_3)].strip()+'n')#改为: fw_1.write('>'+sp[0]+'n'+each_seq[len(d_1)+int(d_2):len(d_1)+int(d_3)].strip()+'n')可以省略第二步(脚本二),一步完成

           break            
fr.close()
fa.close()
fw_1.close()

表头没有'>',同时也没有换行处理,所以需要继续处理(图五):

没有写连续的脚本,重新写了一个(脚本二):

import re
fr=open(r'F:desktopfa_3.txt','r')
fw=open(r'F:desktopfa_4.fa','w')
line_fr=fr.readlines()
s_1=''
for eachline in line_fr:
   s_1=re.sub('t','n',eachline)
   fw.write(re.sub('pp','>pp',s_1))
fr.close()
fw.close()

最终得到:

程序比较简单,Python初学者都可以懂。当然,如果有错误的地方可以留言指出,希望能为需要的同学提供帮助。这个程序只是针对于正链的

注:之前写的出现了一个bug,经过修改后发布成功提取序列,希望对各位有帮助,有用的话可以引用。

鉴于某些人盗版,转载请注明网址。




 




https://blog.sciencenet.cn/blog-783116-801490.html


下一篇:Endnote插入文献末尾显示句号
收藏 IP: 122.205.72.*| 热度|

1 李志远

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 03:40

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部