5

File 1 (contain only selected IDs):

AAX50250
AAX50251
AAX50252
AAX50253
AAX50254
AAX50257
AAX50258

File 2 (Contains sequences with the IDs). This is a FASTA file where lines beginning with as > are the sequence ID and the subsequent lines are the sequence itself.

> AAX50250.1 cds:annotated chromosome:ASM1212v1:Chromosome:1:1770:1 gene:CTA_0001.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGAGCATCAGGGGAGTAGGAGGCAACGGGAATAGTCGAATCCCTTCTCATAATGGGGAT
GGATCGAATCGCAGAAGTCAAAATACGAAGAATAAAGTTGAAGATCGAGTTCGTTCTCTA
TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT
> AAX50251.1 cds:annotated chromosome:ASM1212v1:Chromosome:1915:2187:-1 gene:CTA_0002.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical membrane associated protein
ATGCTTTGTAAAGTTTGTAGAGGATTATCTTCTCTTATTGTTGTTCTCGGAGCTATAAAC
ACTGGAATTTTAGGAGTAACAGGGTATAAGGTAAACCTACTTACTCACCTGCTTGGTGAA
GGAACCATGTGGACACAAGCAGCTTATGTAGTAACGGGAATCGCTGGGGTTATGGTCTGC
> AAX50252.1 cds:annotated chromosome:ASM1212v1:Chromosome:2388:2690:1 gene:CTA_0003.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:gatC description:glutamyl-tRNA(Gln) amidotransferase subunit C
ATGACAGAGTCATATGTAAACAAAGAAGAAATCATCTCTTTAGCAAAGAATGCTGCATTG
GAGTTGGAAGATGCCCACGTGGAAGAGTTCGTAACATCTATGAATGACGTCATTGCTTTA
ATGCAGGAAGTAATCGCGATAGATATTTCGGATATCATTCTTGAAGCTACAGTGCATCAT
TTCGTTGGTCCAGAGGATCTTAGAGAAGACATGGTGACTTCGGATTTTACTCAAGAAGAA
TAG
> AAX50262.1 cds:annotated chromosome:ASM1212v1:Chromosome:13318:13932:-1 gene:CTA_0013.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:ybbP description:hypothetical membrane spanning protein
ATGTTCGTAGGTATAACGTATTACACCACACCTCTGTTGGAGATAGCTTTAATTTGGGTG
GTCCTTAATTATTTGCTAAAGTTTTTCTGGGGAACAGGCGCCATGGACCTCGTCTTTGGC
TTGTTGTCTTTTCTTTGCCTATTTGTTCTAGCAGAAAAACTTCATCTCCCCGTTATTCGC
AATTTGATGCTTCATGTAGTGAATATTGCGGCTATCGTGGTATTTATTATCTTCCAACCA
GAAATTCGCCTTGCTCTCTCTAGGATACGCTTGTGTAGAGAGAAATTTGTCATCAATATG

So now I have to compare File 1 which contains the selected ids with File2 which has n number of sequences. We need to output the sequences associated with the IDs in file 1.

For example, the ID AAX50250 which is in file 1 should return:

>AAX50250.1 cds:annotated chromosome:ASM1212v1:Chromosome:1:1770:1 gene:CTA_0001.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGAGCATCAGGGGAGTAGGAGGCAACGGGAATAGTCGAATCCCTTCTCATAATGGGGAT
GGATCGAATCGCAGAAGTCAAAATACGAAGAATAAAGTTGAAGATCGAGTTCGTTCTCTA
TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT
ATGATCGAGAGCACCCCAACGAGTGGAGAGACGACAAGAGCTTCGCGTGGAGTATTCAGT
CGTTTCCAAAGAGGTTTAGGACGAGTAGCTGACAAAGTAAGACGAGCTGTTCAGCGTGCG
TGGAGTTCAGTCTCTATAAGAAGATCGTCTGCAACAAGAGCCGCAGAATCCAGATCAAGT

How can I do this?

  • Does each id starts with new line? – snoop Oct 3 '16 at 8:51
  • 1
    @snoop yes, lines starting with > are ID lines. See updated question. Nitha, in the future, please use the formatting tools to format your posts. Also, don't assume people know what FASTA is, you need to explain what the input file is. This isn't a site for bioinformaticians, only for Linux geeks, so don't expect us to know the specifics of this sort of thing :). – terdon Oct 3 '16 at 8:54
  • AAX50250.1 cds:a ..........TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT is this complete 1 sentence? – snoop Oct 3 '16 at 8:54
  • @snoop that's a complete entry, yes. The line starting with a > is the ID line, and the rest (until the next >, or the end of the file) is the sequence. – terdon Oct 3 '16 at 9:02
  • DNA processing with unix text tools, amazing – cat Oct 3 '16 at 12:43
6

You could write a script to do this directly, but I suggest you use the two little scripts below instead. I have spent many years working with this sort of data and they have been very, very useful. If you're planning to do a lot of FASTA file manipulation, I suggest you use them.

  1. FastaToTbl:

    #!/usr/bin/awk -f
    {
            if (substr($1,1,1)==">")
            if (NR>1)
                        printf "\n%s\t", substr($0,2,length($0)-1)
            else 
                printf "%s\t", substr($0,2,length($0)-1)
            else 
                    printf "%s", $0
    }END{printf "\n"}
    
  2. TblToFasta:

    #! /usr/bin/awk -f
    {
      sequence=$NF
    
      ls = length(sequence)
      is = 1
      fld  = 1
      while (fld < NF)
      {
         if (fld == 1){printf ">"}
         printf "%s " , $fld
    
         if (fld == NF-1)
          {
            printf "\n"
          }
          fld = fld+1
      }
      while (is <= ls)
      {
        printf "%s\n", substr(sequence,is,60)
        is=is+60
      }
    }
    

Save the scripts above in your $HOME/bin directory (if it doesn't exist, create it then log out and log back in to add it to your $PATH). Then make the scripts executable with chmod a+x ~/bin/TblToFasta ~/bin/FastaToTbl.

The FastaToTbl will convert FASTA sequence files to tbl format (the sequence ID, followed by a tab character and the sequence, all in one line). Once they're in tbl format and you have the ID and the sequence on the same line, you can use grep to search for your IDs. Then, you pass the output through TblToFasta toi convert back to FASTA:

$ FastaToTbl file2 | grep -wFf file1 | TblToFasta
>AAX50250.1 cds:annotated chromosome:ASM1212v1:Chromosome:1:1770:1 gene:CTA_0001.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein 
ATGAGCATCAGGGGAGTAGGAGGCAACGGGAATAGTCGAATCCCTTCTCATAATGGGGAT
GGATCGAATCGCAGAAGTCAAAATACGAAGAATAAAGTTGAAGATCGAGTTCGTTCTCTA
TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT
>AAX50251.1 cds:annotated chromosome:ASM1212v1:Chromosome:1915:2187:-1 gene:CTA_0002.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical membrane associated protein 
ATGCTTTGTAAAGTTTGTAGAGGATTATCTTCTCTTATTGTTGTTCTCGGAGCTATAAAC
ACTGGAATTTTAGGAGTAACAGGGTATAAGGTAAACCTACTTACTCACCTGCTTGGTGAA
GGAACCATGTGGACACAAGCAGCTTATGTAGTAACGGGAATCGCTGGGGTTATGGTCTGC
>AAX50252.1 cds:annotated chromosome:ASM1212v1:Chromosome:2388:2690:1 gene:CTA_0003.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:gatC description:glutamyl-tRNA(Gln) amidotransferase subunit C 
ATGACAGAGTCATATGTAAACAAAGAAGAAATCATCTCTTTAGCAAAGAATGCTGCATTG
GAGTTGGAAGATGCCCACGTGGAAGAGTTCGTAACATCTATGAATGACGTCATTGCTTTA
ATGCAGGAAGTAATCGCGATAGATATTTCGGATATCATTCTTGAAGCTACAGTGCATCAT
TTCGTTGGTCCAGAGGATCTTAGAGAAGACATGGTGACTTCGGATTTTACTCAAGAAGAA
TAG

The grep options used are:

   -F, --fixed-strings
          Interpret PATTERN as a list of fixed strings  (instead  of  regular
          expressions), separated by newlines, any of which is to be matched.

   -f FILE, --file=FILE
          Obtain  patterns  from  FILE, one per line.  If this option is used
          multiple times or is combined with the -e (--regexp) option, search
          for all patterns given.  The empty file contains zero patterns, and
          therefore matches nothing.
   -w, --word-regexp
          Select only those lines containing matches that form  whole  words.
          The  test  is  that  the  matching  substring must either be at the
          beginning of the  line,  or  preceded  by  a  non-word  constituent
          character.   Similarly, it must be either at the end of the line or
          followed by a  non-word  constituent  character.   Word-constituent
          characters are letters, digits, and the underscore.

So, the -f lets us give it a file to read the search [patterns from, the -F tells it to treat the patters as strings and not regular expressions (otherwise, if your IDs contain . which they often do, that will be taken as "match any character") and the -w ensures we only match the entire ID, so that foo will not match foobar.


Alternatively, here's a quick and dirty perl one liner to do what you want:

perl -e 'open(F,"File1.txt");while(<F>){/(\S+)/; $k{$1}++}; while(<>){if(/>\s*(\S+?)(\.| )/){if($k{$1}){$k=1}else{$k=0}; } print if $k==1;}' File2.fa

Or, slightly less dirty:

perl -e 'open(F,"File1.txt");
         while(<F>){
            /(\S+)/; 
            $k{$1}++
         } 
         while(<>){
            if(/>\s*(\S+?)(\.| )/){
                if($k{$1}){$k=1}
                else{$k=0}
            } 
            print if $k==1
        }' File2.fa
3

How about something like this?

$ awk -F'[ .]' 'NR==FNR{a[$0]; next}/^>/{p=$2 in a}p' file1 file2
> AAX50250.1 cds:annotated chromosome:ASM1212v1:Chromosome:1:1770:1 gene:CTA_0001.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein
ATGAGCATCAGGGGAGTAGGAGGCAACGGGAATAGTCGAATCCCTTCTCATAATGGGGAT
GGATCGAATCGCAGAAGTCAAAATACGAAGAATAAAGTTGAAGATCGAGTTCGTTCTCTA
TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT
> AAX50251.1 cds:annotated chromosome:ASM1212v1:Chromosome:1915:2187:-1 gene:CTA_0002.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical membrane associated protein
ATGCTTTGTAAAGTTTGTAGAGGATTATCTTCTCTTATTGTTGTTCTCGGAGCTATAAAC
ACTGGAATTTTAGGAGTAACAGGGTATAAGGTAAACCTACTTACTCACCTGCTTGGTGAA
GGAACCATGTGGACACAAGCAGCTTATGTAGTAACGGGAATCGCTGGGGTTATGGTCTGC
> AAX50252.1 cds:annotated chromosome:ASM1212v1:Chromosome:2388:2690:1 gene:CTA_0003.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:gatC description:glutamyl-tRNA(Gln) amidotransferase subunit C
ATGACAGAGTCATATGTAAACAAAGAAGAAATCATCTCTTTAGCAAAGAATGCTGCATTG
GAGTTGGAAGATGCCCACGTGGAAGAGTTCGTAACATCTATGAATGACGTCATTGCTTTA
ATGCAGGAAGTAATCGCGATAGATATTTCGGATATCATTCTTGAAGCTACAGTGCATCAT
TTCGTTGGTCCAGAGGATCTTAGAGAAGACATGGTGACTTCGGATTTTACTCAAGAAGAA
TAG

The script above can be expanded to:

awk -F'[ .]' '       # set the field separator to space or dot
        NR == FNR {  # for the lines of the first file
            a[$0]    # Store the line in an array
            next     # And skip to the next record
        }
        /^>/{        # For lines of the second file starting with `>`
            p = $2 in a # Set flag to true if index is present in array
        }
        p            # print line if flag is set
      ' file1 file2
  • @terdon and @ user000001 so whether this one line command solves my problem awk -F'[ .]' 'NR==FNR{a[$0]; next}/^>/{p=$2 in a}p' file1 file2 I tried it but I couldn't get any clear out put – Nitha Oct 3 '16 at 9:28
  • @Nitha yes, this should work on your example files but it might fail if your input isn't exactly as you show. For example, it will fail if you include the accession version (the .N) in your list of IDs. – terdon Oct 3 '16 at 10:26
  • @terdon, I have tried this command but unfortunately I am not getting the output. " awk -F'[ .]' 'NR==FNR{a[$0]; next}/^>/{p=$2 in a}p' CAI_Filtered_ID-A.txt (file1) CHT_A.fa (file2) >test.txt " but in output file test.txt I am not getting any out put. could you plz tell me where I am doing mistake. I have given files as i shown in example. File 1- list of IDs and File 2- consit of fasta sequences – Nitha Oct 3 '16 at 10:48
  • @Nitha we can't discuss this in the comments. Either edit your question and add this information, or ping me in chat. – terdon Oct 3 '16 at 10:51
  • @Nitha 1) please stop asking questions in comments! Ask a new question about this. 2) Of course there's no grep command for that. grep is a tool tat searches for patterns, it can't do any of the things you ask for. 3) what you need is a script. Start with my FastaToTbl command: FastaToTbl file.fa | awk _F"\t" '{print length($2)}' and try to write one yourself. Then, ask a question either here or on Unix & Linux about any specific problem you have. – terdon Oct 28 '16 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.