2

I would like to use grep in order to copy from one file to another all lines that are between lines /protein_id= untill the end of the protein sequence shown. For example, from this input:

 CDS             448..1269
                 /gene="nptII"
                 /note="neomycin phosphotransferase II"
                 /codon_start=1
                 /product="kanamycin resistance protein"
                 /protein_id="AAQ05967.1"
                 /db_xref="GI:33320494"
                 /translation="MAITLSATSLPISARIRAGSPAAWVERLFGYDWAQQTIGCSDAA
                 VFRLSAQGRPVLFVKTDLSGALNELQDEAARLSWLATTGVPCAAVLDVVTEAGRDWLL
                 LGEVPGQDLLSSHLAPAEKVSIMADAMRRLHTLDPATCPFDHQAKHRIERARTRMEAG
                 LVDQDDLDEEHQGLAPAELFARLKARMPDGEDLVVTHGDACLPNIMVENGRFSGFIDC
                 GRLGVADRYQDIALATRDIAEELGGEWADRFLVLYGIAAPDSQRIAFYRLLDEFF"
  regulatory      1443..2148

I would want this output:

                 /protein_id="AAQ05967.1"
                 /db_xref="GI:33320494"
                 /translation="MAITLSATSLPISARIRAGSPAAWVERLFGYDWAQQTIGCSDAA
                 VFRLSAQGRPVLFVKTDLSGALNELQDEAARLSWLATTGVPCAAVLDVVTEAGRDWLL
                 LGEVPGQDLLSSHLAPAEKVSIMADAMRRLHTLDPATCPFDHQAKHRIERARTRMEAG
                 LVDQDDLDEEHQGLAPAELFARLKARMPDGEDLVVTHGDACLPNIMVENGRFSGFIDC
                 GRLGVADRYQDIALATRDIAEELGGEWADRFLVLYGIAAPDSQRIAFYRLLDEFF"

Note that the input could vary in that the line starting with regulatory could be replaced with something else. What is unchanging is that the sequence is given in capital letters and ends with a ". Is it possible with grep?

  • 2
    Karel's answer is perfect if you insist on grep..although should mention an alternate easier sed way: sed '/db_xref/,/regulatory_class/!d' file.txt – heemayl Mar 25 '15 at 12:09
  • No it isn't. The question has been changed, so I deleted my answer. – karel Mar 25 '15 at 15:25
3

No, grep can't match across multiple lines. You could do it with pcregrep as shown by @karel, but not pure grep. Instead, since you know that the protein sequences will always be in UPPER CASE and will end with ", you could match that:

  1. sed

    sed -n '/\/protein_id=/,/^\s*[[:upper:]]\+"\s*$/{p}' two_seq.txt
    

    The sed pattern /foo/,/bar/{p} means "print all lines between foo and bar. The -n suppresses normal output so only the lines requested are printed. Note that the / of /protein_id= needs to be escaped (\/) because the / is part of the match operator. The second pattern is a bit more complex, it looks for 0 or more spaces at the beginning of the line (^\s*), then one or more upper case letters followed by a double quote ([[:upper:]]") and then 0 or more whitespace characters until the end of the line (\s*$).

  2. Perl

    perl -ne 'print if m#/protein_id=# ... m#[A-Z]+"\s*$#' file.flat 
    

    The same idea here, the ... operator specifies a range and lines between the two patterns are printed.

  3. awk

    awk '/\/protein_id=/{a=1}; a==1{print} /^\s*[[:upper:]]+"\s*$/{a=0}' file.flat 
    

    Here, we're setting the variable a to 1 if the line matches the first pattern and to 0 if it matches the last. Then, we tell awk to print if a is 1. Since the print is called before a is set to 0 for the second pattern, this will include the line containing the second pattern as well.

  • Thanks, I tried all but its not working 100%, in the out file I still have unnecessary lines: /regulatory_class="terminator" /db_xref="GI:33320495" /translation="MREAVIAEVSTQLSEVVGVIERHLEPTLLAVHLYGSAVDGGLKP HSDIDLLVTVTVRLDETTRRALINDLLETSASPGESEILRAVEVTIVVHDDIIPWRYP AKRELQFGEWQRNDILAGIFEPATIDIDLAILLTKAREHSVALVGPAAEELFDPVPEQ DLFEALNETLTLWNSPPDWAGDERNVVLTLSRIWYSAVTGKIAPKDVAADWAMERLPA QYQPVILEARQAYLGQEDRLASRADQLEEFVHYVKGEITKVVGK" 10706..11324 /note="right border" Why is that, is something missing in the sed perl or awk command lines? – Urszula Uciechowska Mar 25 '15 at 13:51
  • @UrszulaUciechowska that's strange. What is your input file? It looks like some kind of flat file format, could you give me a link to the original file so I can test it? You probably have extra characters between the matching patterns. – terdon Mar 25 '15 at 13:53
  • yes, I have some extra lines. I was thinking that it should work. I put the txt file here sendspace.com/file/2o72m3 its two_seq.txt. – Urszula Uciechowska Mar 25 '15 at 14:17
  • @UrszulaUciechowska ah, yes. That's because the input file has multiple occurrences of /db_xref= and /regulatory_class= and I am guessing you only want the ones describing the CDS. It might be easier if you could join me in this chat room so you can explain your requirements ore precisely. You can then come and edit your question to clarify. – terdon Mar 25 '15 at 14:23

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.