#!/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"を出力してから上のエラーを吐き出します。