#下列內(nèi)容保存到list.txt
gene46   NC_000964.3  42917  43660  +
NP_387934.1  NC_000964.3    59504  60070    +
yfmC  NC_000964.3  825787   826734  -
cds821  NC_000964.3  885844  886173 -

# -*- coding: utf-8 -*-
#初始傳遞命令
input_file = 'Bacillus_subtilis.str168.fasta'
list_file = 'list.txt'
output_file = 'gene.fasta'
name_detail = True
##讀取文件
#讀取基因組序列
seq_file = {}
with open(input_file, 'r') as input_fasta:
    for line in input_fasta:
        line = line.strip()
        if line[0] == '>':
            seq_id = line.split()[0]
            seq_file[seq_id] = ''
        else:
            seq_file[seq_id] += line
input_fasta.close()
#讀取列表文件
list_dict = {}
with open(list_file, 'r') as list_table:
    for line in list_table:
        if line.strip():
            line = line.strip().split('\t')
            list_dict[line[0]] = [line[1], int(line[2]) - 1, int(line[3]), line[4]]
list_table.close()
##截取序列并輸出
#定義函數(shù)，用于截取反向互補
def rev(seq):
    base_trans = {'A':'T', 'C':'G', 'T':'A', 'G':'C', 'N':'N', 'a':'t', 'c':'g', 't':'a', 'g':'c', 'n':'n'}
    rev_seq = list(reversed(seq))
    rev_seq_list = [base_trans[k] for k in rev_seq]
    rev_seq = ''.join(rev_seq_list)
    return(rev_seq)
#截取序列并輸出
output_fasta = open(output_file, 'w')
for key,value in list_dict.items():
    if name_detail:
        print('>' + key, '[' + value[0], value[1] + 1, value[2], value[3] + ']', file = output_fasta)
    else:
        print('>' + key, file = output_fasta)
    seq = seq_file['>' + value[0]][value[1]:value[2]]
    if value[3] == '+':
        print(seq, file = output_fasta)
    elif value[3] == '-':
        seq = rev(seq)
        print(seq, file = output_fasta)
output_fasta.close()</pre>

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#導(dǎo)入模塊，初始傳遞命令、變量等
import argparse
parser = argparse.ArgumentParser(description = '\n該腳本用于在基因組特定位置截取序列，需額外輸入記錄有截取序列信息的列表文件', add_help = False, usage = '\npython3 seq_select.py -i [input.fasta] -o [output.fasta] -l [list]\npython3 seq_select.py --input [input.fasta] --output [output.fasta] --list [list]')
required = parser.add_argument_group('必選項')
optional = parser.add_argument_group('可選項')
required.add_argument('-i', '--input', metavar = '[input.fasta]', help = '輸入文件，fasta 格式', required = True)
required.add_argument('-o', '--output', metavar = '[output.fasta]', help = '輸出文件，fasta 格式', required = True)
required.add_argument('-l', '--list', metavar = '[list]', help = '記錄“新序列名稱/序列所在原序列ID/序列起始位置/序列終止位置/正鏈（+）或負鏈（-）”的文件，以 tab 作為分隔', required = True)
optional.add_argument('--detail', action = 'store_true', help = '若該參數(shù)存在，則在輸出 fasta 的每條序列 id 中展示序列在原 fasta 中的位置信息', required = False)
optional.add_argument('-h', '--help', action = 'help', help = '幫助信息')
args = parser.parse_args()
##讀取文件
#讀取基因組序列
seq_file = {}
with open(args.input, 'r') as input_fasta:
    for line in input_fasta:
        line = line.strip()
        if line[0] == '>':
            seq_id = line.split()[0]
            seq_file[seq_id] = ''
        else:
            seq_file[seq_id] += line
input_fasta.close()
#讀取列表文件
list_dict = {}
with open(args.list, 'r') as list_file:
    for line in list_file:
        if line.strip():
            line = line.strip().split('\t')
            list_dict[line[0]] = [line[1], int(line[2]) - 1, int(line[3]), line[4]]
list_file.close()
##截取序列并輸出
#定義函數(shù)，用于截取反向互補
def rev(seq):
    base_trans = {'A':'T', 'C':'G', 'T':'A', 'G':'C', 'a':'t', 'c':'g', 't':'a', 'g':'c'}
    rev_seq = list(reversed(seq))
    rev_seq_list = [base_trans[k] for k in rev_seq]
    rev_seq = ''.join(rev_seq_list)
    return(rev_seq)
#截取序列并輸出
output_fasta = open(args.output, 'w')
for key,value in list_dict.items():
    if args.detail:
        print('>' + key, '[' + value[0], value[1] + 1, value[2], value[3] + ']', file = output_fasta)
    else:
        print('>' + key, file = output_fasta)
    seq = seq_file['>' + value[0]][value[1]:value[2]]
    if value[3] == '+':
        print(seq, file = output_fasta)
    elif value[3] == '-':
        seq = rev(seq)
        print(seq, file = output_fasta)
output_fasta.close()

#python3 seq_select.py -h
python3 seq_select.py -i Bacillus_subtilis.str168.fasta -l list.txt -o gene.fasta --detail