Как редактировать скрипт Python? - PullRequest
Новая
       27

Как редактировать скрипт Python?

0 голосов
/ 06 декабря 2018

Я испытываю трудности при попытке отредактировать рабочий скрипт на Python.

У меня есть 2 файла:

  1. .txt-файл, содержащий идентификаторы
  2. .fasta-файл, в котором есть последовательности Fasta с их идентификаторами.

Цель этого сценария - сравнить 2 файла, и как только идентификатор из первого файла совпадает с последовательностью, а его идентификатор - из второго файла, на выходе должен быть идентификатор, полная последовательность иего идентификатор.

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

введите описание изображения здесь

Это скрипт: 

with open('uniprot_reviewed_taxonomy_9606.fasta', 'r') as f:
    for line in f.readlines():
        line = line.replace("\n", "")
        if line.startswith('>'):
            full_name = line.split('|')
            accession_x = full_name[1]
            print(accession_x)
        else:
           print (line)

            with open('homosapiens_output1.txt', 'r') as f1:
                for line1 in f1.readlines()[1:]:  # ignores the first line
                    line1 = line1.replace("\n", "")

                    full_name1 = line1.split(' ')
                    accession_y = full_name1[0].replace(" ", "")
                    accession_z = full_name1[1].replace(" ", "")
                    main_accession = accession_x + " " + accession_z + " " + accession_y + " " + line

                    if accession_x == accession_z:
                        print(main_accession)

Так что вы можете помочь мне, отредактировав скрипт такможет на выходе будет идентификатор, последовательность Фаста и его идентификатор?

1 Ответ

0 голосов
/ 06 декабря 2018

Вот небольшой тестовый пример, как это сделать.

Код: 

# create test fasta:
test_content = """>tr|Q53XC5|Q53XC5_HUMAN Bone morphogenetic protein 4 OS=Homo sapiens OX=9606 GN=BMP4 PE=2 SV=1
MIPGNRMLMVVLLCQVLLGGASHASLIPETGKKKVAEIQGHAGGRRSGQSHELLRDFEAT
LLQMFGLRRRPQPSKSAVIPDYMRDLYRLQSGEEEEEQIHSTGLEYPERPASRANTVRSF
HHEEHLENIPGTSENSAFRFLFNLSSIPENEVISSAELRLFREQVDQGPDWERGFHRINI
YEVMKPPAEVVPGHLITRLLDTRLVHHNVTRWETFDVSPAVLRWTREKQPNYGLAIEVTH
LHQTRTHQGQHVRISRSLPQGSGNWAQLRPLLVTFGHDGRGHALTRRRRAKRSPKHHSQR
ARKKNKNCRRHSLYVDFSDVGWNDWIVAPPGYQAFYCHGDCPFPLADHLNSTNHAIVQTL
VNSVNSSIPKACCVPTELSAISMLYLDEYDKVVLKNYQEMVVEGCGCR
>tr|A8K571|A8K571_HUMAN Bone morphogenetic protein 7 (Osteogenic protein 1), isoform CRA_b OS=Homo sapiens OX=9606 GN=BMP7 PE=2 SV=1
MHVRSLRAAAPHSFVALWAPLFLLRSALADFSLDNEVHSSFIHRRLRSQERREMQREILS
ILGLPHRPRPHLQGKHNSAPMFMLDLYNAMAVEEGGGPGGQGFSYPYKAVFSTQGPPLAS
LQDSHFLTDADMVMSFVNLVEHDKEFFHPRYHHREFRFDLSKIPEGEAVTAAEFRIYKDY
IRERFDNETFRISVYQVLQEHLGRESDLFLLDSRTLWASEEGWLVFDITATSNHWVVNPR
HNLGLQLSVETLDGQSINPKLAGLIGRHGPQNKQPFMVAFFKATEVHFRSIRSTGSKQRS
QNRSKTPKNQEALRMANVAENSSSDQRQACKKHELYVSFRDLGWQDWIIAPEGYAAYYCE
GECAFPLNSYMNATNHAIVQTLVHFINPETVPKPCCAPTQLNAISVLYFDDSSNVILKKY
RNMVVRACGCH
>tr|A8K660|A8K660_HUMAN Adiponectin OS=Homo sapiens OX=9606 GN=ADIPOQ PE=2 SV=1
MLLLGAVLLLLALPGHDQETTTQGPGVLLPLPKGACTGWMAGIPGHPGHNGAPGRDGRDG
TPGEKGEKGDPGLIGPKGDIGETGVPGAEGPRGFPGIQGRKGEPGEGAYVYRSAFSVGLE
TYVTIPNMPIRFTKIFYNQQNHYDGSTGKFHCNIPGLYYFAYHITVYMKDVKVSLFKKDK
AMLFTYDQYQENNVDQASGSVLLHLEVGDQVWLQVYGEGERNGLYADNDNDSTFTGFLLY
HDTN
"""

with open('test_sequences.fasta', 'w') as f:
    f.write(test_content)

# create test ids:
test_ids = 'A8K660\nQ53XC5\n'

with open('test_ids.txt', 'w') as f:
    f.write(test_ids)

# Load all sequences and store them in dict:

with open('test_sequences.fasta', 'r') as f:
    lines = f.read().split('>')

sequences = {}
for seq in lines:
    if seq:
        id_ = seq[: seq.find('\n')].split('|')[1]
        seq = seq[seq.find('\n')+1:]
        sequences[id_] = seq

# import ids:
with open('test_ids.txt', 'r') as f:
    ids = f.readlines()
    ids = [id_.strip() for id_ in ids]  # remove \n from id end

# checking ids
# and filter out those that are not in sequences dict

filtered_ids = [id_ for id_ in ids if id_ in sequences.keys()]

# writing new file with filtered sequences:

with open('filtered_ids.txt', 'w') as f:
    for id_ in filtered_ids:
        f.write('>|' + id_ + '\n')
        f.write(sequences[id_])

# the final function:

def ids_filter(ids_file, seq_file, out_file):
    with open(seq_file, 'r') as f:
        lines = f.read().split('>')

    sequences = {}
    for seq in lines:
        if seq:
            id_ = seq[: seq.find('\n')].split('|')[1]
            seq = seq[seq.find('\n')+1:]
            sequences[id_] = seq

    with open(ids_file, 'r') as f:
        ids = f.readlines()
        ids = [id_.strip() for id_ in ids]

    filtered_ids = [id_ for id_ in ids if id_ in sequences.keys()]

    with open(out_file, 'w') as f:
        for id_ in filtered_ids:
            f.write('>|' + id_ + '\n')
            f.write(sequences[id_])
...