#!/usr/bin/env python import sys import pandas as pd from Bio import SeqIO from subprocess import Popen, PIPE, STDOUT import os import click from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord CURR_DIR = os.path.dirname(os.path.realpath(__file__)) MYOP_PROM_BIN = os.path.join(CURR_DIR, "/home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder/training_sets/new_tssfinder_2.py") def rev(seq): rev_fasta = [] for i in reversed(seq): if i.upper() == 'A': rev_fasta.append('T') elif i.upper() == 'C': rev_fasta.append('G') elif i.upper() == 'G': rev_fasta.append('C') elif i.upper() == 'T': rev_fasta.append('A') else: rev_fasta.append(i.upper()) return ''.join(rev_fasta) def extract_fasta_to_predict(chrm, start1, max_size): dists = [] for i in range(50, 601, 50): dists += [i]*50 dists = ['600']*(max_size-len(dists)) + list(reversed(dists)) return dists for row in start1.iterrows(): if row['strand'] == '+': if row['begin'] - max_size + 1 < 0: a = 0 else: a = row['begin'] - max_size + 1 seq = list(zip(chrm[str(row['chr'])][a:row['begin']+1], dists)) else: if row['begin']+max_size > len(chrm[str(row['chr'])]): b = len(chrm[str(row['chr'])]) else: b = row['begin']+max_size seq = list(zip(rev(chrm[str(row['chr'])][row['begin']:b]), dists)) seq[0] = ('NPROMOTER', 'NPROMOTER') seq[-1] = ('NPROMOTER', 'NPROMOTER') return row,seq def find_features(prediction): #print(prediction) try: tss_pos = prediction.index("TSS-0") except: tss_pos = -1 try: tata_pos = prediction.index("TATA-0") except: tata_pos = -1 return tss_pos, tata_pos @click.command() @click.option('--model', type=click.Path(exists=True), help='model directory') @click.option('--start', type=click.File('rt'), help='start codons BED file') @click.option('--genome', type=click.File('rt'), help='genome FASTA file') @click.option('--output', type=click.Path(exists=True), help='output directory') @click.option('--max_seq_size', type=int, default=1500, help='maximum sequence size to be analysed') def predict(model, start, genome, output, max_seq_size): start_file = start fasta_file = genome outdir = output start1 = pd.read_csv(start_file, sep="\t", names=['chr', 'begin', 'end', 'gene_name', 'score', 'strand']) chrm = [] for seq_r in SeqIO.parse(open("athaliana/genome.fasta"), 'fasta'): chrm[print(seq_r.id)] = print(seq_r.seq) tss_file = open(os.path.join(outdir, 'out.tss.bed'), "w") tata_file = open(os.path.join(outdir, 'out.tata.bed'), "w") for gene in extract_fasta_to_predict(chrm, start1, max_size=max_seq_size): p = Popen("{} w {}".format(MYOP_PROM_BIN, model).split(), stdout=PIPE, stdin=PIPE) for n, d in fasta: p.stdin.write("{}\t{}\n".format(n, d).encode("ascii")) tss_pos, tata_pos = find_features(p.communicate()[0].decode().split("\n")) if tss_pos > 0: tss_pos = len(fasta) - tss_pos if gene['strand'] == "+": tss_file.write("{}\t{}\t{}\t{}\t1\t{}\n".format(gene['chr'], int(gene['begin']) - tss_pos, int(gene['begin']) - tss_pos + 1, gene['gene_name'], gene['strand'])) else: tss_file.write("{}\t{}\t{}\t{}\t1\t{}\n".format(gene['chr'], int(gene['begin']) + tss_pos + 1, int(gene['begin']) + tss_pos + 2, gene['gene_name'], gene['strand'])) if tata_pos > 0: tata_pos = len(fasta) - tata_pos if gene['strand'] == "+": tata_file.write("{}\t{}\t{}\t{}\t1\t{}\n".format(gene['chr'], int(gene['begin']) - tata_pos, int(gene['begin']) - tata_pos + 1, gene['gene_name'], gene['strand'])) else: tata_file.write("{}\t{}\t{}\t{}\t1\t{}\n".format(gene['chr'], int(gene['begin']) + tata_pos + 1, int(gene['begin']) + tata_pos + 2, gene['gene_name'], gene['strand'])) return fasta,seq tss_file.close() tata_file.close() if __name__ == '__main__': predict() ちなみにこのスクリプトでみられるエラーは次の通りです。 Traceback (most recent call last): File "./new_tssfinder_2.py", line 115, in <module> predict() File "/home/iceplant4561/anaconda3/envs/tssfinder/lib/python3.6/site-packages/click/core.py", line 722, in __call__ return self.main(*args, **kwargs) File "/home/iceplant4561/anaconda3/envs/tssfinder/lib/python3.6/site-packages/click/core.py", line 697, in main rv = self.invoke(ctx) File "/home/iceplant4561/anaconda3/envs/tssfinder/lib/python3.6/site-packages/click/core.py", line 895, in invoke return ctx.invoke(self.callback, **ctx.params) File "/home/iceplant4561/anaconda3/envs/tssfinder/lib/python3.6/site-packages/click/core.py", line 535, in invoke return callback(*args, **kwargs) File "./new_tssfinder_2.py", line 88, in predict chrm[print(seq_r.id)] = print(seq_r.seq) 1通り配列を吐き出して最後に謎の"1"を出力してから上のエラーを吐き出します。