最近，用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，经过修改后发布成功提取序列，希望对各位有帮助，有用的话可以引用。

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

转载本文请联系原作者获取授权，同时请注明本文来自王彬忠科学网博客。

链接地址：http://blog.sciencenet.cn/blog-783116-801490.html