1

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()
4
  • 発生しているエラーの種類は何でしょうか?(Value Error とか Type Error とかのこと)
    – merino
    Commented 2022年11月25日 23:54
  • ちなみに動かないと言われている書き直していない初期状態でのエラーはどんなものだったか、それは何が原因か、公開されている動作環境との違いはどんなものか、といった情報も出来る限り追記してみてください。
    – kunif
    Commented 2022年11月26日 0:31
  • 1
    こちら バイオインフォマティクス武者修行 #4 ~TSSfinder~ だれかたすけて が内容的に質問者さんご本人かもしれない記事のようですが、質問と記事の内容を同期させるかリンクを貼って重点部分を解説するとか工夫してみてください。あとnoteの記事のPythonプログラムのインデントが一部分失われて実態と合っていないので修正した方が良いでしょう。それからnoteの記事もそうですが、使っているPC・OS・言語・IDE・仮想環境とかの種別や版数の情報も記述した方が良いでしょう。
    – kunif
    Commented 2022年11月26日 10:04
  • rev() 関数(逆鎖を求める関数?)内の return ''.join(rev_fasta) ですが、インデントの深さが for ループの内部と同じになっています。そのため、リードの最初の塩基だけが返されることになります。
    – metropolis
    Commented 2022年11月27日 4:59

1 件の回答 1

0

質問の全体的な内容がバイオインフォマティクス武者修行 #4 ~TSSfinder~ だれかたすけてと同じであると仮定すると、現在発生している問題は質問への転記不足?である最後の行でしょう。

  File "./new_tssfinder_2.py", line 88, in predict
    chrm[print(seq_r.id)] = print(seq_r.seq)
TypeError: list indices must be integers or slices, not NoneType

note記事でも幾つか見受けられますが、print(...)の戻り値はNoneであり、その結果を有効な値として使用することは出来ません。
だからひと通り配列を吐き出しての部分は上記のprint(seq_r.id)print(seq_r.seq)が行われた結果かもしれません。ただしループの最初で直ぐにエラーになるはずなのでひと通りというのは謎ですが。
(あるいは質問やnote記事には書かれなかったか、使っている関数/メソッドが自動的にか、データを表示しているのかも? そして最後に謎の"1"を出力というのがループの最初のprint(seq_r.seq)の部分で、その結果(None)をchrm[print(seq_r.id)] =で代入しようとしてエラーになっているのかも?)

note記事の下のように書き直すに書かれた以下の部分(も少し間違っている訳ですが)と対比させると:

start = pd.read_csv(start_file, sep="\t", names=['chr', 'begin', 'end', 'gene_name', 'score', 'strand'])

chrm = {}
for f in SeqIO.parse(open(fasta_file), 'fasta'):
    name = print(f.id)
    chrm[name] = str(f.seq)

ソースコードの以下の部分を

    chrm = []
    for seq_r in SeqIO.parse(open("athaliana/genome.fasta"), 'fasta'):
        chrm[print(seq_r.id)] = print(seq_r.seq)

こちらのように変更して試してみてはどうでしょう?

    chrm = {}  #### リストではなく辞書として初期化
    for seq_r in SeqIO.parse(open("athaliana/genome.fasta"), 'fasta'):
        chrm[str(seq_r.id)] = str(seq_r.seq)  #### print() ではなく str() を使う

他にも色々と問題が潜んでいそうですが、取り敢えずは何かしらの変化があると思われます。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。