Selv om løsningen fra https://bioinformatics.stackexchange.com/users/96/daniel-standage skal fungere (etter justering for mulig python3-inkompatibilitet), er følgende kortere og mindre minne krevende metode som bruker biopython:

  #! / usr / bin / env pythonfra Bio import SeqIOimport syslastGene = Nonelongest = (None, None) for rec in SeqIO.parse (sys.argv [1], "fasta"): gene = ".". Join (rec.id.split (".") [: - 1]) l = len (rec) hvis lastGene ikke er Ingen: hvis gen == lastGene: hvis lengst [0] < l: lengst = (l, rec) annet: lastGene = gen SeqIO.write (lengst [1], sys.stdout, "fasta") lengst = (l, rec) else: lastGene = gen longest = (l, rec) SeqIO.write (longest [1], sys.stdout, "fasta")  

Hvis du lagret dette som filter.py , så ville filter.py original.fa > subset.fa være kommandoen å bruke.

Her er en løsning i R. Kan bli veldig treg med store filer. Fungerer for eksemplet du la ut.

  bibliotek (Biostrings) ## les fasta i som Biostrings objectfasta.s <- readDNAStringSet ("sample.fa") ## få lesenavnene (i saken din har isoform info) names.fasta <- names (fasta.s) ## trekker bare ut det aktuelle genet og isoform id (delt navn med periodesymbolet) gene.iso <- sapply (names.fasta, function ( j) cbind (unlist (strsplit (j, '\\.')) [2: 3])) ## konverter til god data.frame = transponere resultat fra forrige trinn og legg til relevante kolonnen namesgene.iso.df <- data .frame (t (gene.iso)) colnames (gene.iso.df) <- c ('gen', 'isoform') ## og lengden på isoformsgene.iso.df $ bredde <- bredde (fasta.s) ## split data.frame in list with entry for each genegene.iso.df.split <- split (gene.iso.df, gene.iso.df $ gen) ## optional for å beholde all informasjonen, men egentlig trenger du bare indekser ## gene.iso.df.split.best <- lapply (gene.iso.df.split, funksjon (x) x [rekkefølge (x $ bredde) [1],]) ## trekk ut den lengste isoform-IDen for hvert gen (i tilfelle uavgjort er det bare å ta den første) best.id <- sapply (gene.iso.df.split, function (x) row.names (x) [order (x $ bredde) [1]]) ## delmengde originalen din leser med delsettfasta.s.best <- fasta.s [best.id] ## eksporter ny hurtigfil som inneholder lengste isoform per genewriteXStringSet (fasta.s, filepath = 'sample_best_isoform .fasta ')  

Sent til festen her, men jeg liker å prøve å unngå å skrive manus når noe magisk kommandolinje vil gjøre. Det er god praksis å indeksere FASTA , så bruk den.

  samtools faidx <myfasta.fa>  

Litt awk og sort kan bestemme den største isoformen for hvert gen (og de trenger ikke engang å bli sortert etter navn i kilden FASTA ).

  awk -F '[\ t.]' '{skriv ut $ 1, $ 2, $ 3, $ 4}' <myfasta.fa>.fai | sorter -k4nr, 4 | sorter -uk1,2 | kutt -f1-3 -d '' | tr '' '.' > selection.ls  

• -F [\ t.] spesifiserer både fai -avgrensere, og hele stopper i kontignavnene dine ( Doug_NoIndex_L005_R1_001_contig_2.g7.t1 ) som avgrensere for awk.
• For hver linje skriver vi ut fire kolonner. $ 1 til $ 3 er de tre komponentene i kontinnavnet ditt (etter å ha delt på . ), og kolonne $ 4 er kolonnen i indeksen for sekvenslengde.
• Vi bruker sort for å organisere linjene etter $ 4 (contiglengden). Omvendt sorter med -r slik at de lengste genene vises først.
• Vi sorterer deretter etter kolonner $ 1 og $ 2 (den kontig og gennavn ( g1 ...)). -u slipper deretter enhver gjentagelse av kolonne $ 1 og $ 2 par (det vil si alle linjene etter det første paret - den lengste lengden).
• Til slutt setter vi sammen $ 1 til $ 3 med cut og tr for å gi deg en liste over sekvensnavn.

En sløyfe med vår gamle venn samtools faidx kan sette dine valgte sekvenser ut til en ny FASTA .

  mens du leser kontig; do samtools faidx <myfasta.fa> $ contig >> selection.fa; done < selection.ls  

Ikke glem å indeksere det, de er nyttige!

  samtools faidx selection.fa  

En tidligere kollega (Josep Avril) har skrevet et par veldig nyttige små manus som konverterer fasta til tbl ( seqID<TAB>Sequence ) og tilbake igjen. Disse er ekstremt praktiske for denne typen ting (og er inkludert på slutten av dette svaret). Ved å bruke dem kan du konvertere fasta til en sekvens per linjeformat, beholde den lengste sekvensen med et enkelt awk-skript og konvertere tilbake til fasta.

Jeg antar at delen etter den siste . i ID-linjen bør fjernes (slik at Doug_NoIndex_L005_R1_001_contig_2.g7.t1 og Doug_NoIndex_L005_R1_001_contig_2.g7.t2 begge tilordnes til Doug_NoIndex_L005__ / kode>). Hvis det er en riktig antagelse, bør dette fungere for deg:

  $ FastaToTbl file.fa | sed 's / \. [^.] * \ t / \ t /' | awk -F "\ t" '{if (length ($ 2) >length (a [$ 1])) {a [$ 1] = $ 0}} END {for (i in a) {print a [i]}}' | TblToFasta>Doug_NoIndex_L005_R1_001_contig_2.g7 atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg  

Et mulig problem med denne tilnærmingen er at det er behov for å holde en sekvens per ID i minnet. Hvis du har en enorm fasta-fil og ikke mye minne, kan det være et problem. Hvis det er tilfelle, kan du først sortere filen for å sikre at lignende ID-er alltid vises sammen, og deretter skrive ut hver linje så snart du finner neste ID:

  FastaToTbl file.fa | sed 's / \. [^.] * \ t / \ t /' | LC_ALL = C sorter -t $ '\ t' -k1,1 | awk -F "\ t" '{if (NR>1 && prev! = $ 1) {skriv ut [$ 1]; forrige = $ 1}
hvis (lengde ($ 2) >length (a [$ 1])) {a [$ 1] = $ 0}} END {print a [$ 1]} '| TblToFasta  

• FastaToTbl

    #! / Usr / bin / awk -f {if (substr ($ 1,1,1) == ">") hvis (NR>1) printf "\ n% s \ t", substr ($ 0,2, lengde ($ 0) -1) annet printf "% s \ t ", substr ($ 0,2, lengde ($ 0) -1) annet printf"% s ", $ 0} SLUT {printf" \ n "}  

• TblToFasta

    #! / usr / bin / awk -f {sekvens = $ NF ls = lengde (sekvens) er = 1 fld = 1 mens (fld < NF) {if (fld == 1) {printf ">"} printf "% s" , $ fld hvis (fld == NF-1) {printf "\ n"} fld = fld + 1} mens (er < = ls) {printf "% s \ n", substr (sekvens, er, 60) er = er + 60}}  

Denne løsningen fungerer for ditt eksempel og sannsynligvis for alle Fastas fra Augustus, men kjørelengde vil variere utover det.

  #! / usr / bin / env python fra __future__ import print_functionimport sysdef parse_fasta (data): "" "Stjålet skamløst fra http://stackoverflow.com/a/7655072/459780." "" navn, seq = Ingen, [] for linje i data: line = line.rstrip () if line.startswith ('>'): if name: yield (name, '' .join (seq)) name, seq = line, [] else: seq.append (line) if name : yield (name, '' .join (seq)) isoforms = dict () for defline, sequence in parse_fasta (sys.stdin): geneid = '.'. join (defline [1:]. split ('.') [: -1]) hvis genid i isoformer: otherdefline, othersequence = isoformer [geneid] hvis len (sekvens) > len ​​(othersequence): isoformer [geneid] = (defline, sekvens) annet: isoformer [geneid] = (defline, sekvens) f eller defline, sekvens i isoforms.values ​​(): skriv ut (defline, sequence, sep = '\ n')  

Lagre som longest.py , og påkalle på kommandolinjen slik.

  python longest.py < intput.fasta > output.fasta  

Linjen geneid = '.'. join (defline [1:]. split ('.') [: - 1]) er nøkkelen. La meg bryte det ned.

• defline [1:] : ignorere det første tegnet i defline ( > -symbolet)
• defline [1:]. split ('.') : del streng på . symbol
• defline [ 1:]. Split ('.') [: - 1] : ignorere den siste verdien etter splittelsen (isoformnavnet)
• '.'. Join (defline [ 1:]. Split ('.') [: - 1]) : bli med i split verdiene igjen

Denne verdien vil være den samme for alle isoformer for det samme genet , så vi bruker det i skriptet for å holde rede på den lengste sekvensen assosiert med hvert gen-ID.