Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • perl script to exact sequences by name list

    Dear All,
    i am a perl beginner. Recently i wrote a script to extract multiple sequences from a fasta file. The following is an example.

    FASTA FILE:
    >gene1
    cgcggccccccccccccccccccc
    >gene2
    ggggggggggggggggggcgcgcg
    >gene3
    tagctacgatagcgtcagatgcaacg
    >gene4
    attttttttttttttttttttttttttttaaaaaaa

    LIST FILE:
    >gene4
    >gene3

    MY SCRIPT:
    #!/usr/bin/perl
    use strict;
    use warnings;

    open SEQ, "seq.txt";
    my @seq=<SEQ>;
    open LIST, "list.txt";
    my @list=<LIST>;

    my $i=0;
    while ($i<$#seq){
    my $j=0;
    while ($j<$#list){
    if ($seq[$i] eq $list[$j]){
    print "$seq[$i]\t$seq[$i+1]";
    }
    $j++;
    }
    $i +=2;
    }
    close SEQ;
    close LIST;

    THE OUTPUT RESULT:
    >gene4
    attttttttttttttttttttttttttttaaaaaaa

    Could anyone help to correct my script or make suggestions? Thank you very much for your help!!

  • #2
    Originally posted by pony2001mx View Post


    Code:
    	my $j=0;
    		while ($j<$#list){
    		if ($seq[$i] eq $list[$j]){
    		print "$seq[$i]\t$seq[$i+1]";
    		}
    		$j++;
    	}
    In your second while loop, you should have

    Code:
          while($j <= $#list){
    Alternatively, you could use a foreach loop to go through the @list array.

    Comment


    • #3
      This isn't quite the right way to do this. It's way too brittle.

      If you have @list, do

      Code:
      my %hash = map {$_, 1} @list.
      Or, I think slurping a whole file at one go can be dicey if the file is huge, so it's a bit safer to do

      Code:
      my %hash = ();
      while (<LIST>) {
           $hash{$_}++;
      }
      Those will both make a hash where each key is a gene name, and all the values are one.

      Then go through the sequences file, ask what $hash{$seq_name} is. If it's 1, you want that sequence, because its on your list. If it doesn't exist, that seq name was not on your list.

      Comment


      • #4
        Here is a script I put together some time ago that follows swbarnes logic. It reads the list file and stores the IDs in a hash. It then parses through the FASTA file one record at a time and checks if the ID of the current record is in the list. It operates in two modes, include (i) or exclude (e) (defaults to include). In include mode it will output those IDs which are in your list file; in exclude mode it will only output those IDs not in your list file.

        Usage:

        Code:
        % subSetFasta_Simple.pl <-f|--fasta sequenceFile.fasta> <-l|--list listFile.txt> [-m|--mode <i|e>]
        
        subSetFasta_simple.pl takes three arguments:
             -f|--fasta the name of the FASTA input file
             -l|--list the name of the list file. One ID per line, do not include ">"
             -m|--mode either 'i' or 'e' (default 'i')
        Attached Files

        Comment


        • #5
          or just

          Code:
          grep -A 1 -f listFile.txt file.fasta > list.fasta
          Doesn't work if there are linebreaks in the sequences. Also, ridiculously slow..
          savetherhino.org

          Comment


          • #6
            Sometime in the distant past, I installed meme, which has a fasta-fetch program with the following syntax:

            Code:
                    fasta-fetch <fasta> [-f <file> | [<id>]+] [-c] [-s <suf>] 
                            [-off <off>] [-len <len>]
            
                            <fasta>         name of FASTA sequence file
                            [-f <file>]     file containing list of sequence identifiers
                            [<id>]          sequence identifier
                            [-c]            check that first word in fasta id matches <id>
                            [-s <suf>]      put each sequence in a file named "after" the
                                            sequence identifier with ".<suf>" appended; 
                                            pipes in file names are changed to underscores
                            [-off <off>]    print starting at position <off>; default: 1
                            [-len <len>]    print up to <len> characters for each seq; 
                                            default: print entire sequence
            
                    Note: Assumes and index file has been made using fasta-make-index.
                    Sequence identifiers must be same as made by fasta-make-index.
            
                    Reads sequence identifiers from the command line, from
                    a file and from standard input, in that order.
            
                    Fetch sequences from a FASTA sequence file and 
                    write to standard output.
            In response to your perl code, it would be better to process your files line-by-line, and use hashes for lookup, because then your code is more generalisable in the future (e.g. for extracting raw reads out of a 100GB FASTA file from a NGS sequencer):

            Code:
            my %sequenceIDs = ();
            while(<handle>){
              if(/>(.*?) /){
                $sequenceIDs{$1} = 1;
              }
            }
            instead of the "try to read everything into memory at once" approach:
            Code:
            @sequenceIDs = <handle>;
            ... in other words, what swbarnes said.
            Last edited by gringer; 11-11-2013, 08:17 PM.

            Comment


            • #7
              THANKS a lot for your suggestions, which are really helpful!

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Strategies for Sequencing Challenging Samples
                by seqadmin


                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                03-22-2024, 06:39 AM
              • seqadmin
                Techniques and Challenges in Conservation Genomics
                by seqadmin



                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                Avian Conservation
                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                03-08-2024, 10:41 AM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 06:37 PM
              0 responses
              11 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, Yesterday, 06:07 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-22-2024, 10:03 AM
              0 responses
              51 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-21-2024, 07:32 AM
              0 responses
              67 views
              0 likes
              Last Post seqadmin  
              Working...
              X