Python - dictionary maken van een '.fasta'-file

[Noordpoollicht]

Deze vraag is eigenlijk een soort vervolg op de vorige vraag die ik gesteld had.

Ik moet een dictionary maken van een '.fasta'-file van een paar proteïnen in een organisme.

Ik kom steeds op het probleem dat de verschillende lines in de '.fasta'-file niet worden samengevoegd in één sequentie zoals het hoort. Dit heb ik al, maar wat er nu staat geeft de values als '[]' weer (en dat weet ik ook niet hoe ik heb gedaan haha):

Python:

#Een dictionary aanmaken:
dictionary = {}

#empty
ID = ""
sequentie = ""

with open("6EP.fasta","r") as f:
    for line in f:

        if line.startswith(">"):
            ID_plit = line.strip().split("|")
            ID = ID_gesplitst[1].strip()
            if ID not in dictionary:
                dictionaryID] = []
            continue
        ??

print(dictionary)

De '.fasta'-file is opgebouwd uit FASTA-formaten, dit zijn een op tekst gebaseerd formaat voor het weergeven van nucleotidesequenties of aminozuursequenties (eiwitsequenties). De ID van het proteïne wordt standaard weergeven met het '>' teken en daarachter de hele specifieke ID (bestaande uit accesion number, proteïne naam, etc.). Daaronder wordt dan de sequentie van dat proteïne weergeven.

Het lastige is ook dat niet elke proteïne een even grote sequentie heeft en al de aparte lines van de sequenties worden gescheiden door een 'enter'. Eén sequentie heeft dus meerdere lines in het bestand, en deze is bij elke proteïne van verschillende grootte.

Ik moet dus door het hele bestand een loop doen zodanig dat die elke keer de ID, die met '>' begint. En dan de sequentie eronder samenvoegt tot één element (niet aparte lines) en die toevoegt als de value. Dit moet elke keer herhaalt worden per ID en sequentie paar.

Ik ben een absolute beginner.

Alvast bedankt voor jullie hulp!

[Voor 77% gewijzigd door Noordpoollicht op 18-12-2022 23:27]

[eheijnen]

Het bestand begint met een > en vervolgens doe je een split op |.

Kun je aub. aan de topic start uitgebreid en duidelijk toevoegen hoe het bestand is opgebouwd.
Wat waar staat en wat daarvan doorzocht moet worden naar de waarde(s) waarin je geïnteresseert bent?

[Noordpoollicht]

.

[Voor 99% gewijzigd door Noordpoollicht op 18-12-2022 16:41]

[eheijnen]

Als ik naar deze 2 regels kijk dan zouden dat zaken zijn die je zoekt. In ieder geval in dit voorbeeld.

>sp|A1BCD2|EF123_APTFO protein OS=Aptenodytes forsteri OX=9233
>tr|A0A087QGI9|A0A087QGI9_APTFO Neuroblast differentiation-associated protein AHNAK (Fragment)

In principe kun je per regel door het bestand lopen (itereren) en per regel kijken of deze bv. begint met > of dat bevat.

Dan is de vraag tot hoever moet er gelezen worden. Maw. waar begint de tekst en waar eindigd die. Daar zijn functies voor.

code:

Hier kun je de positie van > bepalen (index 0)
En zou je kunnen zoeken naar het volgende scheidingsteken | (index 3) 
dan weet je dat je van 1 t/m 2 sp of tr kunt uitlezen.

>sp|A1BCD2|EF123_APTFO protein OS=Aptenodytes forsteri OX=9233

Een regeleinde is een "\n" die kun je ook detecteren. En ook vervangen met bv. niets / een lege string "". Als je dat doet krijg je 1 regel tekst

Zoeken naar een karakter in een string om de positie daarvan te bepalen.
https://docs.python.org/2.7/library/string.html#string.find

[superduper]

Is dit een opdracht die je per se zo moet insteken of is dit je eign design? Voor de DNA/RNA/Prot zijn diverse tools beschikbaar om een genoom build + index te maken in een enorm efficientere manier dan even een dictinary van maken. Dit gaat je namelijk zo enorm veel geheugen kosten..

[eheijnen]

@superduper
Fasta formaat is (zonet opgezocht) een standaard
Wikipedia: FASTA format

[Noordpoollicht]

Op zich is de verandering van de ID naam niet de hoofdzaak want ik kan inderdaad achteraf ook gewoon gebruik maken van string indexing om de juiste weergave te krijgen.

De hoofdzaak is om een dictionary aan te maken die de ID weergeeft als key en de sequentie als value. Alleen kom ik dan op het probleem dat de sequentie telkens wordt gescheiden door de 'enters' en niet als één element in de values van de dictionary terecht komen.

Elke line die dus met '>' begint moet in de ID key en al de lines eronder (de sequentie dus) moeten samengevoegt worden en in de sequentie value. Dit wil ik doen zodanig de dictionary elke ID koppelt aan de sequentie die er bijhoort (dictionary = ['ID' : 'sequentie']).

[Voor 116% gewijzigd door Noordpoollicht op 18-12-2022 23:22]

[superduper]

eheijnen schreef op zondag 18 december 2022 @ 16:58:
@superduper
Fasta formaat is (zonet opgezocht) een standaard
Wikipedia: FASTA format

Daar weet ik alles van; werk er al zo'n 25 jaar mee Daarom vraag ik TS ook of hij daadwerkelijk het wiel opnieuw wil uitvinden, en er een vierkant van wil maken.

[eheijnen]

@Noordpoollicht
Daarom kun je die enters ("\n" ) vervangen met niets (een lege string ) of wat je graag hebt.
Daarna ga je aan de gang met de string te doorzoeken.

Dus eerst opschonen en dan op zoek...
https://docs.python.org/2...tring.html#string.replace

[Voor 13% gewijzigd door eheijnen op 18-12-2022 17:38]

[Noordpoollicht]

@eheijnen
Bedankt voor de hulp! Ik ga nog wat bij moeten leren over python.

[Voor 80% gewijzigd door Noordpoollicht op 18-12-2022 23:22]

[eheijnen]

Als je hier niet eerder mee gewerkt hebt. Of tenminste een bestand inlezen, doorlopen en daar gericht informatie uithalen dan zul je hier best moeite mee hebben.

Is het zo dat je hier in de toekomst baat bij hebt dan zou ik wat meer tijd in python steken. Misschien is dat ook de bedoeling van de professor.

[Noordpoollicht]

@superduper
Het moet inderdaad verplicht in FASTA-formaat gemaakt worden.

[Voor 117% gewijzigd door Noordpoollicht op 18-12-2022 23:23]

[eheijnen]

Bv.
https://www.w3schools.com/python/default.asp
Of op Geeks for Geeks

Zijn er medestudenten die daar al wat handiger in zijn....praat daar eens mee..

[Noordpoollicht]

.

[Voor 101% gewijzigd door Noordpoollicht op 18-12-2022 23:23]

[Noordpoollicht]

@eheijnen
Dankje voor uw hulp!

[Voor 101% gewijzigd door Noordpoollicht op 18-12-2022 23:24]

[Renzmeister]

Ik heb het even geprobeerd met een defaultdict i.p.v. een standaard dictionary. Is dit ongeveer wat je nodig hebt?

Python:

from collections import defaultdict

dictionary = defaultdict(str)

with open("test.fasta", "r") as file:
    for line in file:
        if line.startswith(">"):
            id = line.strip().split("|")[1]
        else:
            dictionary[id] += line.strip()

Edit: Door defaultdict(str) te gebruiken, voorkom je een KeyError in regel 10. Je kan het ook met een gewone dictionary oplossen, maar dit is wellicht wat 'cleaner'.

[Voor 16% gewijzigd door Renzmeister op 18-12-2022 20:02]

[Noordpoollicht]

@Renzmeister
Dat is bijna het antwoord dat ik zocht! Super bedankt!
Alleen is de eerste ID voor één of andere reden verdwenen is als ik de dictionary.keys() afprint.

edit : blijkbaar print het de eerste sequentie van dictionary.values() ook niet af.

[Voor 28% gewijzigd door Noordpoollicht op 18-12-2022 23:24]

[Renzmeister]

Ja, mijn voorbeeld werkt alleen op de voorbeelden die je in je topicstart hebt geplaatst. Ik zie nu dat het FASTA formaat iets uitgebreider is. Wellicht is de eerste header in een ander formaat dan ">tr|code|etc..."?

[Noordpoollicht]

@Renzmeister
Normaal gezien begint elke identifier met >tr|code|etc...

[Voor 93% gewijzigd door Noordpoollicht op 18-12-2022 23:25]

[T.Kreeftmeijer]

Inhoudelijk is dit ook niet mijn sterkste kant.
Ik ben scheikundestudent en heb best wat programmeerervaring, maar dit soort dingen doe ik niet dagelijks. Het is wel een uitdagende opdracht in ieder geval.

Maar als ik zo even snel op het internet zoek, dan vind ik in ieder geval het volgende:
https://pypi.org/project/fasta/
https://onestopdataanalysis.com/read-fasta-file-python/
Is er een reden dat standaard packages niet gebruikt mogen worden? Ik zag het zo gauw niet staan.

Ik ben wel benieuwd waar dit is, niet dat het goed of slecht is. Het is tenslotte bonus voor de mensen die net wat meer weten of kunnen. Een beetje uitdaging is niet erg.

In ieder geval veel succes ermee.

[Noordpoollicht]

@T.Kreeftmeijer
Dankje voor uw antwoord.
De opdracht moet inderdaad gewoon in standaard python gebeuren.

[Voor 30% gewijzigd door Noordpoollicht op 18-12-2022 23:25]

[Noordpoollicht]

@Renzmeister
Ik heb het proberen schrijven in een standaard dictionary:

Python:

dictionary = {}
accesion_number = ''
sequentie = ''
with open("6EP.fasta", "r") as f:
    for line in f:
        if line.startswith(">"):
            dictionary[accesion_number] = sequence
            sequence = ""
        if line.startswith(">"):
            accesion_number = line.split("|")[1]
        else:
            sequentie = sequentie + line.strip()
    dictionary[accesion_number] = sequentie

Helaas gaf het niet de gewenste output.

[Voor 83% gewijzigd door Noordpoollicht op 18-12-2022 23:28]

[Renzmeister]

Ehhh ja...Het fasta bestand dat je probeert te openen is geen plain text file, maar een RTF bestand... Ik denk dat dat niet klopt! Heb je het zo aangeleverd gekregen? Wat je in ieder geval even kan doen, is het bestand openen in Word, en dan alle inhoud kopieren en plakken in Kladblok, en dit opslaan als txt file. Maar als je het zo aangeleverd hebt gekregen, zou ik eens nagaan of dat wel klopt, want een RTF bestand parsen lijkt me een beetje ver buiten scope van je opdracht.

[Noordpoollicht]

@Renzmeister
Ah ja, ik zie net dat ik het verkeerde bestand ermee heb geopend.

[Voor 41% gewijzigd door Noordpoollicht op 18-12-2022 23:28]

[Renzmeister]

Yup. Kijk nog eens goed naar je if-statements (waarom heb je er 2? Ze zijn identiek...) Je initieert je variabelen met lege strings. De eerste keer dat je script een ">" tegenkomt in het bestand, doet hij letterlijk:

Python:

dictionary[''] = ''

Ik denk dat je daar nog eens naar de volgorde van je statements moet kijken.

[Ben(V)]

OK ik heb een een simpel voorbeeld voor je gemaak

code:

dictionary = {}

with open("6EP.fasta", "r") as f:
    for Line in f:
        if Line.startswith('>'):
            sequentie = Line[1:-1]
            dictionary[sequentie] = ''
        else:
            dictionary[sequentie] += Line[:-1]

SortedDict = dict(sorted(dictionary.items()))
for Ids, Seq in SortedDict.items():
    print( f"{Ids} = {Seq}")

Uitleg:
Line 4:
Lees de file per regel.

Line 5, 6 en 7:
Als regel begint met ">" dan is het een sequence en maak je een nieuw leeg item in de dict aan.
Het eerste en het laatste karakter laat je weg dat zijn de ">" en de "\n" (end of line ) die wil je niet ( Line[1:-1] )

Line 9:
Alle regels die niet met een "<" beginnen zijn onderdeel van de informatie en die plak je eraan vast met ( "+=") en uiteraard haal je de end of line eraf ( Line[:-1] ).

Line 11:
Hier sorteer je de dict op values
Eerst splits je de dict op in key, value ( dictonary.items() )
Dan sorteer je ze ( sorted() )
En maakt er weer een dict van ( dict() )

Line 12 en 13
Loop door de gesorteerde dict heen en print ze netjes uit.

[Voor 6% gewijzigd door Ben(V) op 19-12-2022 11:08]