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
                Essential Discoveries and Tools in Epitranscriptomics
                by seqadmin


                The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist on Modified Bases...
                Yesterday, 07:01 AM
              • seqadmin
                Current Approaches to Protein Sequencing
                by seqadmin


                Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
                04-04-2024, 04:25 PM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, 04-11-2024, 12:08 PM
              0 responses
              39 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 04-10-2024, 10:19 PM
              0 responses
              41 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 04-10-2024, 09:21 AM
              0 responses
              35 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 04-04-2024, 09:00 AM
              0 responses
              55 views
              0 likes
              Last Post seqadmin  
              Working...
              X