#!/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()
TSSfinderが動かないので自分で書き直していますがエラーが止まりません。やっとこの辺りまで来ました。以下のスクリプトのどこがおかしいかどなたかご教示ください。
ちなみにこのこのスクリプトでみられるエラーは次の通りです。
ひと通り配列を吐き出して最後に謎の"1"を出力してから下のエラーを吐き出します。
エラーメッセージ:
1通り配列を吐き出して最後に謎の"1"を出力してから上のエラーを吐き出します。現状のコード:
#!/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()