我有这样的序列(超过9000):
>TsM_000224500
MTTKWPQTTVTVATLSWGMLRLSMPKVQTTYKVTQSRGPLLAPGICDSWSRCLVLRVYVDRRRPGGDGSLGRVAVTVVETGCFGSAASFSMWVFGLAFVVTIEEQLL
>TsM_000534500
MHSHIVTVFVALLLTTAVVYAHIGMHGEGCTTLQCQRHAFMMKEREKLNEMQLELMEMLMDIQTMNEQEAYYAGLHGAGMQQPLPMPIQ
>TsM_000355900
MESGEENEYPMSCNIEEEEDIKFEPENGKVAEHESGEKKESIFVKHDDAKWVGIGFAIGTAVAPAVLSGISSAAVQGIRQPIQAGRNNGETTEDLENLINSVEDDL
包含">"的行是ID,带字母的行是氨基酸(aa)序列。我需要删除(或移动到另一个文件)的序列低于40 aa和超过4000 aa。然后,生成的文件应该只包含这个范围内的序列(>= 40aa和<= 4K aa)。
我试着写了下面的脚本:
def read_seq(file_name):
with open(file_name) as file:
return file.read().split('n')[0:]
ts = read_seq("/home/tiago/t_solium/ts_phtm0less.txt")
tsf = open("/home/tiago/t_solium/ts_secp-404k", 'w')
for x in range(len(ts)):
if ([x][0:1] != '>'):
if (len([x]) > 40 or len([x]) < 4000):
tsf.write('%sn'%(x))
tsf.close()
print "OK!"
我已经做了一些修改,但我得到的都是空文件或所有的+9000序列。
在for循环中,由于使用range()
(即0,1,2,3,4...
), x
是一个迭代整数。试试这个:
for x in ts:
这将给你ts
中的每个元素作为x
同样,你不需要在x
周围加上括号;Python可以自己遍历字符串中的字符。当你在字符串周围加上括号时,你将它放入一个列表中,因此,如果你尝试,例如,获取x
: [x][1]
中的第二个字符,Python将尝试获取x
所在列表中的第二个元素,并且会遇到问题。
编辑:要包含id,请尝试:
注意:我还将if (len(x) > 40 or len(x) < 4000)
更改为if (len(x) > 40 and len(x) < 4000)
-使用and
而不是or
会给你你正在寻找的结果。
for i, x in enumerate(ts): #NEW: enumerate ts to get the index of every iteration (stored as i)
if (x[0] != '>'):
if (len(x) > 40 and len(x) < 4000):
tsf.write('%sn'%(ts[i-1])) #NEW: write the ID number found on preceding line
tsf.write('%sn'%(x))
试试这个,简单易懂。它不会将整个文件加载到内存中,而是逐行迭代文件。
tsf=open('output.txt','w') # open the output file
with open("yourfile",'r') as ts: # open the input file
for line in ts: # iterate over each line of input file
line=line.strip() # removes all whitespace at the start and end, including spaces, tabs, newlines and carriage returns.
if line[0]=='>': # if line is an ID
continue # move to the next line
else: # otherwise
if (len(line)>40) or (len(line)<4000): # if line is in required length
tsf.write('%sn'%line) # write to output file
tsf.close() # done
print "OK!"
仅供参考,如果在unix环境中工作,您也可以使用awk作为一行解决方案:
cat yourinputfile.txt | grep -v '>' | awk 'length($0)>=40' | awk 'length($0)<=4000' > youroutputfile.txt