GenBank Perl help...

akreibich07 · #1 (**permalink**) 05-24-2009

Hey guys,

I'm doing some Perl scripting for genomic data out of GenBank files...I have to extract the name of the plant, the file name, the number of bases, and all of the genes including their starting and ending positions...for example, with this GenBank file,

LOCUS NC_010093 153819 bp DNA circular PLN 07-MAY-2009
DEFINITION Acorus americanus chloroplast, complete genome.
ACCESSION NC_010093
VERSION NC_010093.1 GI:161622288
DBLINK Project:27981
KEYWORDS .
SOURCE chloroplast Acorus americanus
ORGANISM Acorus americanus
Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
Spermatophyta; Magnoliophyta; Liliopsida; Acoraceae; Acorus.
REFERENCE 1 (bases 1 to 153819)
AUTHORS Peery,R.M., Chumley,T.W., Kuehl,J.V., Boore,J.L. and Raubeson,L.A.
TITLE The complete chloroplast genome of Acorus americanus
JOURNAL Unpublished
REFERENCE 2 (bases 1 to 153819)
CONSRTM NCBI Genome Project
TITLE Direct Submission
JOURNAL Submitted (03-DEC-2007) National Center for Biotechnology
Information, NIH, Bethesda, MD 20894, USA
REFERENCE 3 (bases 1 to 153819)
AUTHORS Peery,R.M., Chumley,T.W., Kuehl,J.V., Boore,J.L. and Raubeson,L.A.
TITLE Direct Submission
JOURNAL Submitted (09-NOV-2007) Department of Biological Sciences, Central
Washington University, 400 E University Way, Ellensburg, WA
98926-7537, USA
COMMENT PROVISIONAL REFSEQ: This record has not yet been subject to final
NCBI review. The reference sequence was derived from EU273602.
COMPLETENESS: full length.
FEATURES Location/Qualifiers
source 1..153819
/organism="Acorus americanus"
/organelle="plastid:chloroplast"
/mol_type="genomic DNA"
/db_xref="taxon:263995"
gene complement(join(96591..97384,69098..69211))
/gene="rps12"
/locus_tag="AcamCp045"
/trans_splicing
/db_xref="GeneID:5777700"
CDS complement(join(96591..96616,97153..97384,69098..69211))
/gene="rps12"
/locus_tag="AcamCp045"
/trans_splicing
/note="trans-splices to 5' rps12 exon within the LSC"
/codon_start=1
/transl_table=11
/product="ribosomal protein S12"
/protein_id="YP_001586161.1"
/db_xref="GI:161622289"
/db_xref="GeneID:5777700"
/translation="MPTIKQLIRNTRQPIRNVTKSPALRGCPQRRGTCTRVYTITPKK
PNSALRKVARVRLTSGFEITAYIPGIGHNLQEHSVVLVRGGRVKDLPGVRYHIVRGTL
DAVGVKDRQQGRSKYGVKKPK"
misc_feature 1..83496
/note="large single copy region (LSC)"
gene complement(136..1197)
/gene="psbA"
/locus_tag="AcamCp001"
/db_xref="GeneID:5777757"

From this I would need to extract the NC number next to Locus for filename, the name next to organism, and the number in front of bp for bases. Then for every time there's a new gene, return what is next to gene=, and then the 2 numbers inside the complement parentheses for the start and end of that gene sequence.

Also, is there a way to make so if I put in the NC number I want, it spits back out the rest of the info...the website this is from is Nucleotide - Acorus americanus chloroplast, complete genome, so if I put in the number I wanted next to the NC, then it would give me the information I want as specified above for that specific NC number....?

I've just been having trouble getting the info extracted out of a GenBank file...any help in the coding for this would be great...thanks!

durden_tyler · #2 (**permalink**) 05-25-2009

Here's one way of doing it:

Code:

$ 
$ perl -ne 'if (/^LOCUS/) {
>   split; printf("Filename = %s\nBase = %s\n",$_[1],$_[2])
> } elsif (/organism/) {                                   
>   s/^\s*\/o/O/; s/=/ = /; print                 
> } elsif (/gene/) {                                       
>   print                                                  
> }' genbank.txt                                           
Filename = NC_010093                                       
Base = 153819                                              
Organism = "Acorus americanus"                             
     gene            complement(join(96591..97384,69098..69211))
                     /gene="rps12"                              
                     /gene="rps12"                              
     gene            complement(136..1197)                      
                     /gene="psbA"                               
                     /gene="psbA"                               
     gene            complement(1388..3959)                     
                     /gene="trnK"                               
                     /gene="trnK"                               
                     /gene="trnK"                               
                     /gene="trnK"                               
...
... (trimmed a lot of lines here for succinctness)
...
     gene            complement(153492..153566)
                     /gene="trnH"
                     /gene="trnH"
     gene            153733..153819
                     /gene="rps19"
$

akreibich07 · #3 (**permalink**) 05-25-2009

Hi tyler_durden,

Well, what's going on is I have over 100 GenBank files, each is it's own NC_ number...(this NC number is just a way to catalogue them). I want to be able to give that NC_ number as a command-line possibly and have it go to the internet, replace that website I posted with whatever NC_ number was given, go into that GenBank file, and return something like:

Name: ___________

FileName: ___________ (NC_######)

Bases: ________

Genes: rps (####, #####) - for example, this would be for the rps gene, and then a new line, and all the rest of the genes are given....if the word complement is present before it, then it should print a negative before the numbers in parentheses: rps -(####, #####)

If it's impossible to use the NC_ number as a command-line, it's okay for me to do it manually, just so I can have a program that can do the output of the name, filename, bases, and genes for me. It would save me loads and loads of time in my research, rather than going through and doing each one by hand. Thanks!

-akreibich07

durden_tyler · #4 (**permalink**) 05-25-2009

Quote:

Originally Posted by akreibich07

...
Name: ___________

(1) Name of what ?
(2) What's the value of this field for your example above ?

Quote:

Genes: rps (####, #####) - for example, this would be for the rps gene, and then a new line, and all the rest of the genes are given....if the word complement is present before it, then it should print a negative before the numbers in parentheses: rps -(####, #####)

This is not clear. You haven't addressed the occurrence of 2 different number ranges. Also, you haven't told what to do if the word "join" is present.

What should be the output for each of the following cases ? (Type in the required output for each case.)

Code:

(3)  gene 151670..151951
(4)  gene complement(join(96591..97384,69098..69211))
(5)  gene join(complement(69098..69211),139932..140725)
(6)  gene complement(153492..153566)

akreibich07 · #5 (**permalink**) 05-25-2009

(1) Name: _________ -- something like "Acorus Americanus"
(2) for situation (3) gene (151670..151951)
(4) gene -(96591..97384, 69098..69211)
(5) gene -(69098..69211), (139932..140725)
(6) gene -(153492..153566)

I'll check on the URL idea you presented...the most important thing to me right now is just getting a script that can be given a whole GenBank file and can give the output I specified. Even if it can only do it for the specific NC_ number I used as an example, if that one works, I can try to find ways to generalize it. Thanks again.

-akreibich07

durden_tyler · #6 (**permalink**) 05-26-2009

Quote:

Originally Posted by akreibich07

...the most important thing to me right now is just getting a script that can be given a whole GenBank file and can give the output I specified....

Ok, so what I did was that I went to the hyperlink: "http://www.ncbi.nlm.nih.gov/nuccore/NC_010093" and clicked "Download" -> "GenBank" and saved the resultant "sequences.gb" text file in my computer.

Assuming that you have figured out a way to download one "sequences.gb" for a given NC_ number, (either manually or via some other tool), here's what you can do:

Code:

$ 
$ head -5 sequences.gb
LOCUS       NC_010093             153819 bp    DNA     circular PLN 07-MAY-2009
DEFINITION  Acorus americanus chloroplast, complete genome.                    
ACCESSION   NC_010093                                                          
VERSION     NC_010093.1  GI:161622288                                          
DBLINK      Project:27981                                                      
$                                                                              
$ tail -5 sequences.gb                                                         
   153661 gattttggaa gaaaattcca attttcaaat acaaaaactt taaactttag gtaaaaaaag    
   153721 gaggaatcag ctgtgacacg ttcattaaaa aaaaatcctt ttgtggccaa tcatttattg    
   153781 gaaaaaattg agaagctcaa catgaaagag gagaaagaa                           
//                                                                             

$ 
$ cat process_genbank.pl
#!/usr/bin/perl -w      
my $element;            
my $in = 0;             
$gbfile = "sequences.gb";
open(GB, $gbfile) or die "Can't open $gbfile: $!";
while (<GB>) {                                    
  chomp($line = $_);                              
  if (/^locus/i) {                                
    @x = split /\s+/,$line;                       
    $fname = $x[1];                               
    $base = $x[2];                                
  } elsif (/ organism /i) {                       
    $line =~ s/\s*ORGANISM\s*//g;                 
    $oname = $line;
  } elsif (/ gene /) {
    @x = split/\s+/,$line;
    if ($x[2] =~ /^\d+..\d+$/) {
      $x[2] = "($x[2])";
    } elsif ($x[2] =~ /^complement\(\d+..\d+\)$/) {
      $x[2] =~ s/complement/-/;
    } elsif ($x[2] =~ /^complement\(join\((\d+..\d+),(\d+..\d+)\)\)$/) {
      $x[2] = "-($1,$2)";
    } elsif ($x[2] =~ /^join\(complement\((\d+..\d+)\),(\d+..\d+)\)$/) {
      $x[2] = "-($1),($2)";
    }
    $element = $x[2];
    $in = 1;
  } elsif (/^\s+\/gene="(.*)"/ and $in == 1) {
    $element = $1.":".$element;
    push @genearray, $element;
    $element = "";
    $in = 0;
  }
}
close(GB) or die "Can't close $gbfile: $!";

# Now start printing the information we gathered
print "Name: $oname\n";
print "FileName: $fname\n";
print "Bases: $base\n";
print "Genes:";
foreach $item (@genearray) {
  @x = split/:/, $item;
  printf("\t%-10s  %s\n",$x[0],$x[1]);
}

$
$

More UNIX and Linux Forum Topics You Might Find Helpful
Thread	Thread Starter	Forum	Replies	Last Post
perl -write values in a file to @array in perl	meghana	Shell Programming and Scripting	27	06-07-2009 06:05 PM
How to Turn perl one-liners into full perl script?	EDALBNUG	UNIX for Dummies Questions & Answers	1	02-04-2009 10:49 AM
[Perl] Accessing array elements within a sed command in Perl script	userix	Shell Programming and Scripting	2	10-03-2008 01:05 PM
[PERL] Running unix commands within Perl Scripts	userix	Shell Programming and Scripting	1	05-28-2008 07:06 PM
Perl: Run perl script in the current process	vino	Shell Programming and Scripting	10	12-09-2005 10:45 AM

Bits Awarded / Charged to durden_tyler for this Post
Date	User	Comment	Amount
05-25-2009	Neo	N/A	50,000

Bits Awarded / Charged to durden_tyler for this Post
Date	User	Comment	Amount
05-25-2009	Neo	N/A	30,000

Thread Tools	Search this Thread
Show Printable Version Email this Page	Search this Thread: Advanced Search
Display Modes	Rate This Thread
Switch to Linear Mode Hybrid Mode Switch to Threaded Mode	Rate This Thread:


	Powered by