TSSfinderが動かないので自分で書き直していますがエラーが止まりません。やっとこの辺りまで来ました。以下のスクリプトのどこがおかしいかどなたかご教示ください。
このスクリプトでみられるエラーは次の通りです。
ひと通り配列を吐き出して最後に謎の"1"を出力してから下のエラーを吐き出します。
エラーメッセージ:
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)
現状のコード:
#!/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()
rev()
関数(逆鎖を求める関数?)内のreturn ''.join(rev_fasta)
ですが、インデントの深さが for ループの内部と同じになっています。そのため、リードの最初の塩基だけが返されることになります。