在对拼装或者数据库下载的序列文件进行下一步分析时，我们通常会对序列进行去冗余操作，其中经常需要提取同一个‘gene’的最长转录本，所以动手用python写一个脚本。
一、基本思路
-
1.以Trinity拼装结果为例，分析其文件结构：
序列名：TRINITY_DN1000_c115_g5_i1
其中TRINITY_DN1000_c115_g5 是‘gene’名称， i1 表示该‘gene’的第一个转录本
>TRINITY_DN1000_c115_g5_i1 len=247 path=[31015:0-148 23018:149-246]
AATCTTTTTTGGTATTGGCAGTACTGTGCTCTGGGTAGTGATTAGGGCAAAAGAAGACAC
ACAATAAAGAACCAGGTGTTAGACGTCAGCAAGTCAAGGCCTTGGTTCTCAGCAGACAGA
AGACAGCCCTTCTCAATCCTCATCCCTTCCCTGAACAGACATGTCTTCTGCAAGCTTCTC
CAAGTCAGTTGTTCACAGGAACATCATCAGAATAAATTTGAAATTATGATTAGTATCTGA
TAAAGCA
-
2.设计提取思路：
技术路线图
其中需要注意的是字典的设置，key为 ‘gene’的名称，而Value为一个列表，一个基因的多个转录本分别为列表中的元素。
二、代码实现
废话不多说，下面开始代码。
-
1.import需要的模块和写注释：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#好像没啥需要特别导入的模块
-
2.打开文件并存入字典
re = {}
with open ('/Users/byeamx/Trinity.fasta') as f:
for line in f:
seq = []
if line.startswith('>'):
id = line.split(' ')[0].split('_') #切片分割序列名称
id = '_'.join(id[:4]) #合并切片的前5部分
else:
seq.append(line)
if id not in re:
re[id] = seq
else:
re[id] += seq
-
3.提取最长转录本
maxseq = {}
for k,v in re.items():
seq = max(v, key=len)
maxseq[k] = seq
-
4.输出为文件
with open('file_cluster.fa','w') as f:
for k,v in maxseq.items():
f.write( k +'\n')
tem = v.__str__()
f.write(tem)
三、一些思考
- 这个脚本只是最基础的功能，只是针对Trinity的转录本，在实际情况中可根据自己的需求更改，甚至加上批量运行语句。
- 通常提取转录本后还要用软件CD-HIT进一步去冗余，当然这是自己情况而定。
- BioPython等模块能方便的对序列进行各种操作，在掌握基础原理后完全可以用Biopython等模块简化我们的工作量。
- 其实在Linux下用awk等等就可以直接实现很多文本处理，我总觉得能在命令行能下直接解决的事情就没必要调用python等等，学好Linux编程也很重。
