在對(duì)拼裝或者數(shù)據(jù)庫(kù)下載的序列文件進(jìn)行下一步分析時(shí)，我們通常會(huì)對(duì)序列進(jìn)行去冗余操作最铁，其中經(jīng)常需要提取同一個(gè)‘gene’的最長(zhǎng)轉(zhuǎn)錄本，所以動(dòng)手用python寫(xiě)一個(gè)腳本垮兑。
一冷尉、基本思路
-
1.以Trinity拼裝結(jié)果為例，分析其文件結(jié)構(gòu)：
序列名：TRINITY_DN1000_c115_g5_i1
其中TRINITY_DN1000_c115_g5 是‘gene’名稱甥角， i1 表示該‘gene’的第一個(gè)轉(zhuǎn)錄本
>TRINITY_DN1000_c115_g5_i1 len=247 path=[31015:0-148 23018:149-246]
AATCTTTTTTGGTATTGGCAGTACTGTGCTCTGGGTAGTGATTAGGGCAAAAGAAGACAC
ACAATAAAGAACCAGGTGTTAGACGTCAGCAAGTCAAGGCCTTGGTTCTCAGCAGACAGA
AGACAGCCCTTCTCAATCCTCATCCCTTCCCTGAACAGACATGTCTTCTGCAAGCTTCTC
CAAGTCAGTTGTTCACAGGAACATCATCAGAATAAATTTGAAATTATGATTAGTATCTGA
TAAAGCA
-
2.設(shè)計(jì)提取思路：
技術(shù)路線圖
其中需要注意的是字典的設(shè)置网严，key為 ‘gene’的名稱识樱，而Value為一個(gè)列表嗤无，一個(gè)基因的多個(gè)轉(zhuǎn)錄本分別為列表中的元素。
二怜庸、代碼實(shí)現(xiàn)
廢話不多說(shuō)当犯，下面開(kāi)始代碼。
-
1.import需要的模塊和寫(xiě)注釋：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#好像沒(méi)啥需要特別導(dǎo)入的模塊
-
2.打開(kāi)文件并存入字典
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.提取最長(zhǎng)轉(zhuǎn)錄本
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)
三割疾、一些思考
- 這個(gè)腳本只是最基礎(chǔ)的功能嚎卫，只是針對(duì)Trinity的轉(zhuǎn)錄本，在實(shí)際情況中可根據(jù)自己的需求更改宏榕，甚至加上批量運(yùn)行語(yǔ)句拓诸。
- 通常提取轉(zhuǎn)錄本后還要用軟件CD-HIT進(jìn)一步去冗余，當(dāng)然這是自己情況而定麻昼。
- BioPython等模塊能方便的對(duì)序列進(jìn)行各種操作奠支，在掌握基礎(chǔ)原理后完全可以用Biopython等模塊簡(jiǎn)化我們的工作量。
- 其實(shí)在Linux下用awk等等就可以直接實(shí)現(xiàn)很多文本處理抚芦，我總覺(jué)得能在命令行能下直接解決的事情就沒(méi)必要調(diào)用python等等倍谜，學(xué)好Linux編程也很重迈螟。
