this is not printing the sequence at the end.. thanks for the reply
Helo all, I wanted to parse aEMBL format like file to fasta. i cannot use bioperl because this is not complete EMBL format. so please suggest me how to get this done..
ID US74811111-0005
OO giensis
OS giensis
SN US74811111
PT I-003, a gene and methods for its use
PA NIX CORPORATION RESEARCH TRIANGLE PARK, NC
PI Carozzi; Nadine (Raleigh, NC); Hargiss; Tracy (Cary, NC); Koziel; Michael G. (Raleigh, NC); Duck; Nicholas B. (Apex, NC); Carr; Brian (Raleigh, NC);
PR 20030828 US20030498518P; 20040826 US20040926819; 20070620 US20070765494;
PE US200304985AN 20070765494
P1 Compositions and methods and seeds are provided.
MDNNPNINECIPYNCLSNPEVEVLGGERIETGYTPIDISLSLTQFLLSEFVPGAGFVLGLVDIIWGIFGPSQWDAFPVQIEQLINQRIEEFARNQAISRLEGLSNLYVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRENPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE
//
The output should be in fasta format which consists of lines starting with ID, PT, PA and Sequence. "//" the two slashes are dividing lines between two EMBL genes.
>US74811111-0005 ; I-003, a gene and methods for its use ; NIX CORPORATION RESEARCH TRIANGLE PARK, NC
MDNNPNINECIPYNCLSNPEVEVLGGERIETGYTPIDISLSLTQFLLSEFVPGAGFVLGLVDIIWGIFGPSQWDAFPVQIEQLINQRIEEFARNQAISRLEGLSNLYVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRENPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE
Like this i have 50,000 sequences in a single file which should be converted to fasta format
6 answers
the following awk script should do the job:
/^ID/ {printf(">%s;",$0); next;}
/^(PT|PA)/ {printf(" %s;",$0); next;}
/^\/\// {printf("\n"); next;}
/^ / {printf("\n%s",substr($0,5)); next;}
{
/* ignore default */
}
END {
printf("\n");
}
> awk -f file.awk file.txt
>ID US74811111-0005 ; PT I-003, a gene and methods for its use; PA NIX CORPORATION RESEARCH TRIANGLE PARK, NC;
MDNNPNINECIPYNCLSNPEVEVLGGERIETGYTPIDISLSLTQFLLSEFVPGAGFVLGLVDIIWGIFGPSQWDAFPVQIEQLINQRIEEFARNQAISRLEGLSNLYVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRENPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE
I've added a 'END' close in the script.
I've added a 'END' statement in the script.
works good but if i have sequence in multiple lines its printing only last line in sequence.. say for example..
"LLGTFDECYPTYLYQKIDESKLKAYTRYQLRGYIEDSQDLEIYLIRYNAKHETVNVPGTGSLWPLSAQSPIGKCGEPNRC APHLEWNPDLDCSCRDGEKCAHHSHHFSLDIDVGCTDLNEDLGVWVIFKIKTQDGHARLGNLEFLEEKPLVGEALARVKR AEKKWRDKREKLEWETNIVYKEAKESVDALFVNSQYDQLQADTNIAMIHAADKRVHSIREAYLPELSVIPGVNAAIFEEL EGRIFTAFSLYDARNVIKNGDFNNGLSCWNVKGHVDVEEQNNQRSVLVVPEWEAEVSQEVRVCPGRGYILRVTAYKEGYG EGCVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRE NPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE"
works good but if i have sequence in multiple lines its printing only last line in sequence.. say for example.. " LLGTFDECYPTYLYQKIDESKLKAYTRYQLRGYIEDSQDLEIYLIRYNAKHETVNVPGTGSLWPLSAQSPIGKCGEPNRC APHLEWNPDLDCSCRDGEKCAHHSHHFSLDIDVGCTDLNEDLGVWVIFKIKTQDGHARLGNLEFLEEKPLVGEALARVKR AEKKWRDKREKLEWETNIVYKEAKESVDALFVNSQYDQLQADTNIAMIHAADKRVHSIREAYLPELSVIPGVNAAIFEEL EGRIFTAFSLYDARNVIKNGDFNNGLSCWNVKGHVDVEEQNNQRSVLVVPEWEAEVSQEVRVCPGRGYILRVTAYKEGYG EGCVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRE NPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE"
Can you use this script to process multiple input files and output to multiple files such as awk -f file.awk *.txt >> *.fasta
This Python script should work? Takes input file as an argument e.g. embllike2fasta.py infile.txt:
import sys
dict = {}
infile = open(sys.argv[1], "r")
outfile = open("outfile.fas", "w")
while 1:
line = infile.readline()
if not line:
break
parts = line.split(None, 1)
if len(parts) == 1 and parts[0] == "//":
outfile.write(">" + dict["ID"] + " ; " + dict["PT"] + " ; " + dict["PA"] + "\n")
outfile.write(dict["seq"] + "\n")
elif len(parts) == 1:
dict["seq"] = parts[0].strip()
else:
dict[parts[0]] = parts[1].strip()
outfile.close()
infile.close()
I haven't tested it, and just hacked it up in the browser on my iPad, so it might not work? Plus it's late :-S lol!
Update: Tested it on my machine this morning, made a few quick edits! Very hacky, but it works :D
How about this? It's remarkably hacky, but should do the trick.
#!/usr/bin/env perl
use strict;
use warnings;
$/ = "//\n";
while (<>) {
chomp;
my %record;
foreach my $line (split /\n/) {
my ($key, $val) = unpack("A4A*", $line);
$record{$key} .= $val;
}
print ">$record{ID} $record{PT} $record{PA}\n$record{''}\n";
}
The important bits are lines 4:($/ = "//n"), which sets the system to read one record at a time instead of a line, and 9:(unpack("A4A*", $line)), which splits the string into the EMBL key and the value.
I am getting this error when i execute this
Can't modify single ref constructor in scalar assignment at embl2fasta.pl line 4, near ""//n";"
I am getting this error when i execute this "Can't modify single ref constructor in scalar assignment at eml2fasa.pl line 4, near ""//n";" "
Hmm. That error normally only occurs when you attempt to assign a value to a reference (e.g. $/ = "//n"; ) If you can show me lines 3-5, I might have a better idea.
I dont know for some reason this code is not working for this.. do you have any reason?? i apologize for this
ID 013789-0068
PS TBD
OO huringiensis
OS ringiensis
OX
SI 68
RA
RL 2010. OKAYAMA UNIVERSITY,JAPAN LAMB CO LTD
FT source 1..1176
MT
AC 67106
SV
CT
PN 013789
PT PROTEIN PRODUCTION METHOD, FUSION PROTEIN, AND ANTISERUM
PA AMA UNIVERSITY,JAPAN LAMB CO LTD.
PI HAYAKAWA TORU (JP) SAKAI, HIROSHI, HAYAKAWA, TORU
P8
P4 10013789
P5 0
PC International Classification: \nUS Classification: \nEuropean Classification: C12N15/62; C07K14/47A25
PR 80199166;
PE 199166
AN 09JP63603
KC 1
P1 ng the DNA into a host bacterium to transform the host bacterium; and (c) causing the expression of the fusion protein in the transformed host bacterium.; The method may further comprise a step of removing the peptide chain (B) from the fusion protein. \n \n
P7
P9 112
PO
PM 10013789;
PB 10013789
PQ 10013789;
EM esentative
W1 PRT
D1 0204
D2 0217
D3 0730
D4 0801
D5 0204
HL [L[P9_GQ;0;3,WO2010013789,45,67]] [L[PM_PN_GQNUC;0;12,WO2010013789]] [L[PQ_PN_GQNUC;0;12,WO2010013789]]
CC mer C1-1-f FH Key Location/Qualifiers Copyright (c)Inc. 2011
LS Application
L2 Publ. Of int. appl. w4
MDNNPNINECIPYNCLSNPEVEVLGGERIETGYTPIDISLSLTQFLLSEFVPGAGFVLGLVDIIWGIFGPSQWDAFPVQI
EQLINQRIEEFARNQAISRLEGLSNLYQIYAESFREWEADPTNPALREEMRIQFNDMNSALTTAIPLLAVQNYQVPLLSV
YVQAANLHLSVLRDVSVFGQRWGFDAATINSRYNDLTRLIGNYTDYAVRWYNTGLERVWGPDSRDWVRYNQFRRELTLTV
LDIVALFSNYDSRRYPIRTVSQLTREIYTNPVLENFDGSFRGMAQRIEQNIRQPHLMDILNSITIYTDVHRGFNYWSGHQ
ITASPVGFSGPEFAFPLFGNAGNAAPPVLVSLTGLGIFRTLSSPLYRRIILGSGPNNQELFVLDGTEFSFASLTTNLPST
IYRQRGTVDSLDVIPPQDNSVPPRAGFSHRLSHVTMLSQAAGAVYTLRAPTFSWQHRSAEFNNIIPSSQITQIPLTKSTN
LGSGTSVVKGPGFTGGDILRRTSPGQISTLRVNITAPLSQRYRVRIRYASTTNLQFHTSIDGRPINQGNFSATMSSGSNL
QSGSFRTVGFTTPFNFSNGSSVFTLSAHVFNSGNEVYIDRIEFVPAEVTFEAEYDLERAQKAVNELFTSSNQIGLKTDVT
DYHIDQVSNLVECLSDEFCLDEKQELSEKVKHAKRLSDERNLLQDPNFRGINRQLDRGWRGSTDITIQGGDDVFKENYVT
LLGTFDECYPTYLYQKIDESKLKAYTRYQLRGYIEDSQDLEIYLIRYNAKHETVNVPGTGSLWPLSAQSPIGKCGEPNRC
APHLEWNPDLDCSCRDGEKCAHHSHHFSLDIDVGCTDLNEDLGVWVIFKIKTQDGHARLGNLEFLEEKPLVGEALARVKR
This is not an answer. Please use comments for discussion or edit the original question.
Or just:
egrep -e '^(ID| )' | sed -e 's/^ID />/g' -e 's/^ //g'
Oh, I see you wanted the contents of PT there too. I leave it as an excercise for the reader! :-)
i dont know how to use sed :(
i ran it but that was not even close to what i need, its just adding ">" symbol in place of ID and rest of input file stays as is
Strange. The egrep in front should pick out only ID and sequence lines from the file, so I don't see how it could retain the rest of the file. The sed stuff replaces ID, but it should also remove initial spaces on the sequence data. But it's just a quick hack anyway, you're probably better off with some other solution.
I guess I should add that I only tested it on the example in the question, so if that isn't representative of the actual data, all bets are off, of course.
This one-liner should do it:
perl -ne 'if ($_ =~ /^ID (.+)/) { print ">$1" } elsif ($_ =~ /^P[TA] (.+)/) { print "; $1 " } elsif ($_ =~ /^\s{4}(.+)/) { print "\n$1\n" }'
Output:
>US74811111-0005 ; I-003, a gene and methods for its use ; NIX CORPORATION RESEARCH TRIANGLE PARK, NC
MDNNPNINECIPYNCLSNPEVEVLGGERIETGYTPIDISLSLTQFLLSEFVPGAGFVLGLVDIIWGIFGPSQWDAFPVQIEQLINQRIEEFARNQAISRLEGLSNLYVTIHEIENNTDELKFSNCVEEEIYPNNTVTCNDYTVNQEEYGGAYTSRNRGYNEAPSVPADYASVYEEKSYTDGRRENPCEFNRGYRDYTPLPVGYVTKELEYFPETDKVWIEIGETEGTFIVDSVELLLMEE
not working for my dataset.. :(
Log in to answer this question.