Linux | Python 命令行

之前学过一段时间python，发现python在处理生物学数据方面非常强大。最近在处理数据的过程中，最近摸索出来的一下不成熟的脚本，记录下来，以便下次使用。

一、根据基因id 提取序列

#!/usr/bin/env python
import sys

USAGE = "\nusage: python %s sequence.fasta  id_list  out.fasta\n" % sys.argv[0]

if len(sys.argv) != 4:
    print (USAGE)
    sys.exit()

def parseFasta(filename):
    fas = {}
    id = None
    with open(filename, "r") as fh:
        for line in fh:
            if line[0] == '>':
                header = line[1:].rstrip()
                id = header.split()[0]
                fas[id] = []
            else:
                fas[id].append(line.strip())
        for id,seq in fas.items():
            fas[id] = ''.join(seq)
    return fas

fas = parseFasta(sys.argv[1])
IDs = open(sys.argv[2], 'r')
OUT = open(sys.argv[3], 'w')

for line in IDs:
    id = line.strip('\n')
    OUT.write(">" + id + "\n" + fas[id] + "\n")

IDs.close()
OUT.close()