バイオインフォマティクス武者修行 #4 ~TSSfinder~ だれかたすけて

こんにちは。
TSSfinderという、高等生物のゲノムに存在するTSS(Transcription Start Site)を特定するソフトウェアがあったので、使おうとしてみました。
途中まではよかったのですが、Pythonの記法で１日以上かかってしまいました。いったん棚に上げてはいますが、あきらめきれず今に至ります。助け舟を求める意味も含めていったんここに投げます。誰か助けて下さい。

WEBもあるのにスタンドアロン版を使う意味としては、より大きなサイズのゲノム配列を扱えるので効率化を図れることが５０％、僕の負けず嫌いが５０％です。

TSSfinderを用いたTSS領域(Transcription Start Site)、すなわちゲノムのコアプロモーターの探索

#あきらめた

参考：https://github.com/tssfinder/tssfinder.github.io

真核生物ゲノムのコアプロモーターの予測を行う TSSfinder - macでインフォマティクス

GFF to BED conversion — AGAT documentation

#ゲノムのgff3ファイルから開始コドンの情報を含むbedファイルの作成
#AGAT(AGAT: Another Gff Analysis Toolkit to handle annotations in any GTF/GFF format)ツールの1つ、agat_convert_sp_gff2bed.pl を使う

#ちなみに 、five-prime UTR情報はなくなるが、開始コドンの情報を追加したgff3ファイルは以下のようにして作る

singularity exec /usr/local/biotools/g/genometools-genometools:1.6.2--py39h0c4336a_2 gt gff3 -tidy -sort -retainids
 Mesembryanthemum_crystallinum.BRAKERv2.1.5_GeMoMa_v1.9.gff3 | 
singularity exec /usr/local/biotools/a/aegean:0.16.0--he819bf6_2 canon-gff3 - > my_with_stops.gff3

#参考 ：https://www.biostars.org/p/236316/
#参考 ：https://manpages.debian.org/testing/aegean/canon-gff3.1.en.html

cd /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/GeMoMA/Final_GFF
singularity exec /usr/local/biotools/a/agat:1.0.0--pl5321hdfd78af_0 agat_convert_sp_gff2bed.pl --gff Mesembryanthemum_crystallinum.BRAKERv2.1.5_GeMoMa_v1.9.gff3 -o Mc_genome.bed

#cd ../../ #mkdir TSSfinder

cd /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder
cp /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/GeMoMA/Final_GFF/Mc_genome.bed ./

#conda create -n tssfinder python=3.6
conda activate tssfinder

#tssfinderのインストール #webではhttps://tssfinder.github.io/download.htmlにアクセス

cd /home/iceplant4561/Important_Software/
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/tssfinder.v1.0.0.zip
unzip tssfinder.v1.0.0.zip
cd tssfinder

#modelのダウンロード
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/athaliana.zip
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/oryza.zip
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/dmelanogaster.zip
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/hsapiens.zip
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/scerevisiae.zip
wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/ggallus.zip
unzip athaliana.zip
unzip oryza.zip

#私の研究では 、植物を扱うので、この２種を用いる。
#それ以外に 、植物種間の解析(cross validation)用のデータ、新たに分かった上記６種のTSSのデータ、サプリメントデータがある。必要に応じてダウンロードすること。

#tssfinderに必要なpythonスクリプトをインストールする

#pip install -r requirements.txt
#tssfinderをPATHに登録する
#echo "export PATH=$PATH:/home/iceplant4561/Important_Software/tssfinder" >> ~/.bashrc
#soource ~/.bashrc
#conda activate tssfinder
cd Important_Software/tssfinder/

#tssfinderというソフトはshellスクリプトである 。毎度のように間違えているので、書き換える。

vim /home/iceplant4561/Important_Software/tssfinder/tssfinder
#[vim⇒Ctrlキー+C長押し⇒:%dで全部消去⇒"i"を押してINSERTモード⇒下の(正)のスクリプトをコピー＆ペースト(Ctrl+C⇒右クリック)⇒Ctrlキー+C長押し⇒:wqで保存]

(誤)

#!/bin/bash
SCRIPT=`readlink -f -- $0`
SCRIPTPATH=`dirname $SCRIPT`
LANG=C.UTF-8 LC_ALL=C.UTF-8 LD_LIBRARY_PATH=$SCRIPTPATH/bin/lib /projects/aliemelo/resources/software/tssfinder/tssfinder.py $@

(正)
#!/bin/bash
SCRIPT=`readlink -f -- $0`
SCRIPTPATH=`dirname $SCRIPT`
LANG=en_US.utf8 LC_ALL=en_US.utf8 LD_LIBRARY_PATH=$SCRIPTPATH/bin/lib /home/iceplant4561/Important_Software/tssfinder/tssfinder.py $@

また、tssfinderのpythonスクリプトも一部変える

vim /home/iceplant4561/Important_Software/tssfinder/tssfinder.py

(誤)
#!/usr/bin/env python
import pandas as pd
from Bio import SeqIO
from subprocess import Popen, PIPE, STDOUT
import os
import click
CURR_DIR = os.path.dirname(os.path.realpath(file))
MYOP_PROM_BIN = os.path.join(CURR_DIR, "bin/cli/tssfinder")

(正)
#!/usr/bin/env python
import pandas as pd
from Bio import SeqIO
from subprocess import Popen, PIPE, STDOUT
import os
import click
CURR_DIR = os.path.dirname(os.path.realpath(file))
MYOP_PROM_BIN = os.path.join(CURR_DIR, "/home/iceplant4561/Important_Software/tssfinder/tssfinder")

cd /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder #戻る

#まさかの話だが、スクリプトが超絶中途半端でした。
書き直します。click関数の指定が甘すぎ。
参考：https://blog.imind.jp/entry/2019/05/25/001742

PythonのclickのType一覧 - iMind Developers Blog

Expected str, bytes or os.PathLike object, not TextIOWrapper | bobbyhadz

下のように書き直す

vim /home/iceplant4561/Important_Software/tssfinder/tssfinder.py

3行目
import sys

65行目
@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 = sys.argv[4]
outdir = output

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)

87行目
p = Popen("{} {}".format(MYOP_PROM_BIN, model).split(), stdout=PIPE, stdin=PIPE)

#予測

#練習
#参考 ：https://tssfinder.github.io/tutorials/01-quick-start.html

#練習用ファイルのダウンロード
#wget https://github.com/tssfinder/tssfinder.github.io/releases/download/v1.0.0/training_sets.zip

#unzip training_sets.zip
cd training_sets/
mkdir athaliana/output_athaliana.model_0 #output先をあらかじめ作っておく

tssfinder --model /home/iceplant4561/Important_Software/tssfinder/athaliana/athaliana.1/ 
--start athaliana/athaliana.dataset0.start.bed 
--genome athaliana/athaliana.tair10.fa 
--output athaliana/output_athaliana.model_0

./tssfinder.sh --model /home/iceplant4561/Important_Software/tssfinder/athaliana 
--start /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder/Mc_genome.bed 
--genome /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/data/Iceplant-genome_fasta_full_softmask.fasta 
--output /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder

tssfinder-train --model /home/iceplant4561/Agarie_group/ice_plant_genome_from_GSA/TSSfinder 
--start Mc_genome.bed 
--tata <TATA-BOX-BED> 
--tss <TSS-BED> 
--genome <GENOME-FASTA> 

※221125-221126の約１日以上かけて試しましたが、エラーだらけでまったく歯が立ちません。
最終的なスクリプトを載せます。
今のとここういうエラーが起こる

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)
TypeError: list indices must be integers or slices, not NoneType

#!/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()

誰か助けてください(´;ω;｀)