Запуск l oop для скрипта python, включающего linux и другие скрипты - PullRequest
Новая
       0

Запуск l oop для скрипта python, включающего linux и другие скрипты

1 голос
/ 23 марта 2020

Я использую простой l oop в скрипте python, чтобы повторять программу столько раз, сколько файлов в папке. Я разрабатываю этот скрипт, поэтому на данный момент у меня есть 3 файла на входе. Так что я ожидаю иметь 3 файла в качестве вывода. По причине, которую я не могу объяснить, я получаю только 1.

Я показываю вам весь свой код, но я указываю, где должна быть проблема

Это некоторые настройки, не беспокойтесь об этом 

import os

#GTDB-Tk requires an environmental variable named GTDBTK_DATA_PATH to be set to the directory containing the data

import_data = "export GTDBTK_DATA_PATH=/Users/monkiky/Desktop/GTDB/gtdbtk/release89"
print(import_data)

# Some third parties are not included in the PATH variable. Let's fix this problem.

third_parties = 'export PATH="/Users/monkiky/opt/anaconda3/lib/python3.7/site-packages/gtdbtk/release89/src:$PATH"'
print(third_parties)

Теперь некоторый код для запуска программы не беспокоит ни 

# First step: Identify The identify step calls genes using Prodigal,
# and uses HMM models and the HMMER package to identify the 120 bacterial
# and marker genes used for phylogenetic inference (Parks et al., 2018)

com1 = "gtdbtk identify --genome_dir /Users/monkiky/Desktop/GTDB/input --out_dir /Users/monkiky/Desktop/control/output"

# Second step: The align step concatenates the aligned marker genes and filters
# the concatenated MSA to approximately 5,000 amino acids.
# We will use a option "--skip_trimming" that provided the full concatenate we need

com2= "gtdbtk align --identify_dir /Users/monkiky/Desktop/control/output --skip_trimming --out_dir /Users/monkiky/Desktop/control/output.align"

path = "/Users/monkiky/Desktop/GTDB/input"
print("cd",path) # To tell the next loop for where it must work.

allfiles = os.listdir(path) # Where the genomes are.

Здесь, где мой код запускает программу, l oop для 

for myfile in allfiles:
    if myfile.endswith(".fna"):
        print(com1)
        print(com2)
        #com2 generates many files I script select what I want and some manipulation throught 
        #the next script.
        print("python /Users/monkiky/Desktop/control/getting_protein.py")
        # Remove all files we dont need
        print("rm -r /Users/monkiky/Desktop/control/output")
        print("mkdir /Users/monkiky/Desktop/control/output")
        print("rm -r /Users/monkiky/Desktop/control/output.align")
        print("mkdir /Users/monkiky/Desktop/control/output.align")

Getting_protein.py - это скрипт, который я покажу вам ниже 

import pandas as pd
import os

#Read the fasta document

def fasta_parser(myfile):
    with open(myfile) as f:

        header = ""
        seq = ""
        for line in f:
            if line[0] == ">":
                if seq != "":
                    yield (header[1:], seq)
                    header = line.strip()
                    seq = ""
                else:
                    header = line.strip()
            else:
                seq += line.strip()
        yield (header[1:], seq)





#Transform the document into a text string called "Sequence"
Sequence =[]
for header, seq in fasta_parser('/Users/monkiky/Desktop/control/output.align/gtdbtk.bac120.user_msa.fasta'):
    print(header,seq[:100])
    Sequence = seq




#The following code saves the length of each protein
f = open("/Users/monkiky/Desktop/GTDB/proteins_len.log", "r")
LEN= []
import re
for line in f:
    secuence = re.search('LEN:(\d+)', line)
    if secuence:
        LEN.append(secuence.group(1))

# To save the name of the protein

f = open("/Users/monkiky/Desktop/GTDB/proteins_len.log", "r")
Names =[]
for i, line in enumerate(f):
    if i > 19:
        Names.append(line.split(" ")[3])


Names = Names[:-9]# Last 9 lines in the docuement are not genes.
#Names
#['Ribosomal_S9:',
# 'Ribosomal_S8:',
# 'Ribosomal_L10:',
# 'GrpE:',
# 'DUF150:',
# 'PNPase:',
# 'TIGR00006:',
# ...


# To create a df with both lists

bac120 = {'protein_name':Names,'LEN':LEN}
df = pd.DataFrame(bac120)
# Calculate the protein sequences in the concatenated
df['LEN'] = df['LEN'].astype(int)
s = df['LEN'].cumsum()
df = df.assign(Start = s.shift(fill_value=0), End = s)
#print (df)

#protein_name  LEN  Start    End
#0     Ribosomal_S9:  121      0    121
#1     Ribosomal_S8:  129    121    250
#2    Ribosomal_L10:  100    250    350
#3             GrpE:  166    350    516
#4           DUF150:  141    516    657
#..              ...  ...    ...    ...
#115      TIGR03632:  117  40149  40266
#116      TIGR03654:  175  40266  40441
#117      TIGR03723:  314  40441  40755
#118      TIGR03725:  212  40755  40967
#119      TIGR03953:  188  40967  41155


#We added the name of the bacterium with the name of the protein
df['protein_name'] = header + ' ' + df['protein_name']



Lets create our fasta file using a dict



mydict = {}

for index,row in df.iterrows():
    mydict[row['protein_name']] =  Sequence[row['Start']:row['End']]


secuencias = [ v for v in mydict.values() ]
nombres = [k for k in mydict]

ofile = open("/Users/monkiky/Desktop/control/ultimate_output/concatenates/my_fasta.fasta", "w")

for i in range(len(secuencias)):

    ofile.write(">" + nombres[i] + "\n" +secuencias[i] + "\n")

ofile.close()


# Remove the "-" and change the name of the final file
import os
with open(r'/Users/monkiky/Desktop/control/ultimate_output/concatenates/my_fasta.fasta', 'r') as infile, open(r'/Users/monkiky/Desktop/control/ultimate_output/concatenates/my_fastaa.fasta', 'w') as outfile:
        data = infile.read()
        data = data.replace("-", "")
        outfile.write(data)
myfile = "/Users/monkiky/Desktop/control/ultimate_output/concatenates/my_fastaa.fasta"
path = '/Users/monkiky/Desktop/control/ultimate_output/concatenates/'
os.rename("/Users/monkiky/Desktop/control/ultimate_output/concatenates/my_fastaa.fasta", "/Users/monkiky/Desktop/control/ultimate_output/concatenates/" + str(header) + ".fasta")

Я действительно не знаю, в чем проблема. Я не знаю, что генерируется только один файл.

Когда я читаю процесс в терминале, я вижу, как в l oop для com1 и com2 работают только 2 раза Почему ??? когда должно быть три, и, наконец, генерируется один файл, когда должно быть три.

Если это поможет, я покажу вам, какой терминал показывает 

(base) monkikys-Mini:control monkiky$ ./commands 
[2020-03-23 13:17:09] INFO: GTDB-Tk v1.0.2
[2020-03-23 13:17:09] INFO: gtdbtk identify --genome_dir /Users/monkiky/Desktop/GTDB/input --out_dir /Users/monkiky/Desktop/control/output
[2020-03-23 13:17:09] INFO: Using GTDB-Tk reference data version r89: /Users/monkiky/Desktop/GTDB/gtdbtk/release89
[2020-03-23 13:17:09] INFO: Identifying markers in 1 genomes with 1 threads.
[2020-03-23 13:17:09] INFO: Running Prodigal V2.6.3 to identify genes.
==> Finished processing 1 of 1 (100.0%) genomes.
[2020-03-23 13:17:22] INFO: Identifying TIGRFAM protein families.
==> Finished processing 1 of 1 (100.0%) genomes.
[2020-03-23 13:17:28] INFO: Identifying Pfam protein families.
==> Finished processing 1 of 1 (100.0%) genomes.
[2020-03-23 13:17:28] INFO: Annotations done using HMMER 3.3 (Nov 2019)
[2020-03-23 13:17:29] INFO: Done.


##### Here com2 finishs

[2020-03-23 13:17:29] INFO: GTDB-Tk v1.0.2
[2020-03-23 13:17:29] INFO: gtdbtk align --identify_dir /Users/monkiky/Desktop/control/output --skip_trimming --out_dir /Users/monkiky/Desktop/control/output.align
[2020-03-23 13:17:29] INFO: Using GTDB-Tk reference data version r89: /Users/monkiky/Desktop/GTDB/gtdbtk/release89
[2020-03-23 13:17:29] INFO: Aligning markers in 1 genomes with 1 threads.
[2020-03-23 13:17:29] INFO: Processing 1 genomes identified as bacterial.
[2020-03-23 13:17:32] INFO: Read concatenated alignment for 23458 GTDB genomes.
==> Finished aligning 1 of 1 (100.0%) genomes.
[2020-03-23 13:17:36] INFO: Skipping custom filtering and selection of columns.
[2020-03-23 13:17:36] INFO: Creating concatenated alignment for 23459 GTDB and user genomes.
[2020-03-23 13:17:51] INFO: Creating concatenated alignment for 1 user genomes.
[2020-03-23 13:17:51] INFO: Done.


##### Here com2 finish egain

GCA_000010565.1_genomic GRRKNAIARVFAMPGEGRIIINNRPLSEYFGRKTLETIVRQPLDLTGTASRFDIMAKVQGGGISGQAGAIKLGIARALIQADPNLRPVLKKAGFLTRDPR

Возможно, проблема в том, что я использую скрипт, который в другом скрипте. Я не знаю. Я новичок в информатике, поэтому я считаю, что это правильный путь.

Любые советы по моей проблеме или второстепенной проблеме в моем коде приветствуются.

Кстати, я также новичок в этом сообществе StackOverflow, если вы обнаружите какую-либо ошибку, например, Как нужно спросить, пожалуйста, дайте мне знать. Рад улучшить и сделать это правильно.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Нет похожих вопросов

...