Keep current and previous line only if current line fulfills a given condition

Viewed 279

I have a file which looks like this:

>4RYF_1
MAENTKNENITNILTQKLIDTRTVLIYGEINQELAEDVSKQLLLLESISNDPITIFINSQGGHVEAGDTIHDMIKFIKPTVKVVGTGWVASAGITIYLAAEKENRFSLPNTRYMIHQPAGGVQGQSTEIEIEAKEIIRMRERINRLIAEATGQSYEQISKDTDRNFWLSVNEAKDYGIVNEIIENRDGLKMASWSHPQFEK
>4RYF_2
MNLIPTVIEQTSRGERAYDIYSRLLKDRIIMLGSAIDDNVANSIVSQLLFLDAQDPEKDIFLYINSPGGSISAGMAIYDTMNFVKADVQTIGMGMAASMGSFLLTAGANGKRFALPNAEIMIHQPLGGAQGQATEIEIAARHILKIKERMNTIMAEKTGQPYEVIARDTDRDNFMTAQEAKDYGLIDDIIINKSGLKGHHHHHH

I want to keep the sequence and previous line only if the sequence has a given length. For selecting only lines with that condition I use:

awk 'length($0) > 50 && length($0) <=800)' sample.txt

But how can I keep lines starting with > as well if this condition is met?

7 Answers

Yet another awk solution:

awk '/^>/ { header = $0; next } length > 50 && length <= 800 { print header ORS $0 }'

Here is one-liner:

LANG=C grep -B1 '^.\{51,800\}$' < sample.txt

The command was really slow with LANG=en_US.UTF-8 that I set by default, so using LANG=C instead.

man grep tells you that '-B NUM' means ' Print NUM lines of leading context before matching lines.'.

'^' means start of line
'.' means any character
'{51,800}' means we want between 51 and 800 of the previous thing
'$' means end of line.

Or in other words, we want to match lines that are between 51 and 800 characters, and print it and the previous line.

Would you please try the following:

awk -v RS='>' -F'\n' '
    length($2) > 50 && length($2) <= 800 {printf ">%s", $0}
' sample.txt
  • Assigning RS to '>' tells awk to split the file on > into records, treating the header line and the sequence line in the same record.

  • Assigning FS to '\n' splits the record to the header and sequence, each assigning $1 to the header and $2 to the sequence.

  • As the leading > is chopped off as a delimiter, we need to prepend it when printing the matched records.

A potential solution with AWK is:

awk '!/^>/ {next}; {getline s}; length(s) > 50 && length(s) <= 800 { print $0 "\n" s }' example.fasta

e.g. if example.fasta contains

>4RYF_1
WLSVNEAKDYGIVNEIIENRDGLKMASWSHPQFEK
>4RYF_2
MNLIPTVIEQTSRGERAYDIYSRLLKDRIIMLGSAIDDNVANSIVSQLLFLDAQDPEKDIFLYINSPGGSISAGMAIYDTMNFVKADVQTIGMGMAASMGSFLLTAGANGKRFALPNAEIMIHQPLGGAQGQATEIEIAARHILKIKERMNTIMAEKTGQPYEVIARDTDRDNFMTAQEAKDYGLIDDIIINKSGLKGHHHHHH
>1000_chars
YiJOgeCApTkcJWxIuvooOxuqVnPdSLtOQmUfnzpBvcpYKyCvelFwKgMchYFnlvuZwVxNcnSvGcACsMywDQVvYBAiaIesQkLkYNsExRbqKPZIPnCRMAFHLmIzxIBqLwoNEPSKMZCTpwbbQCNrHSrbDMtCksTjvQsMeAkoudRGUJnPpQTEzwwnKoZBHtpMSIQBfYSPDYHwKktvCiFpewrsdDTQpqBajOWZkKURaKszEqDmdYMkzSAkMtlkXPfHroiTbyxZwzvrrMSXMRSavrBdgVYZanudjacRHWfpErJMkomXpzagXIzwbaeFgAgFnMxLuQHsdvZysqAsngkCZILvVLaFpkWnOpuYensROwkhwqUdngvlTsXBoCBwJUENUFgVdnSnxVOvfksyiabglFPqmSwhGabjNZiWGyvktzSDOQNGlEvoxhJCAOhxVAtZfyimzsziakpzfIszSWYVgKZTHatWSfttHYTkvgafcsVmitfEfQDuyyDAAAoTKpuhLrnHVFKgmEsSgygqcNLQYkpnhOosKiZJKpDolXcxAKHABtALqVXoVcSHpskrpWPrkkZLTpUXkENhnesmoQjonLWxkpcuJrOosXKNTDNuZaWIEtrDILXsIFTjAnrnwJBoirgNHcDURwDIzAXJSLPLmWkurOhWSLPrIOyqNvADBdIFaCGoZeewKleBHUGmKFWFcGgZIGUdOHwwINZqcOClPAjYaLNdLgDsUNCPwKMrOXJEyPvMRLaTJGgxzeoLCggJYTVjlJpyMsoCRZBDrBDckNMhJSQWBAxYBlqSpXnpmLeEJYirwjfCqZGBZdgkHzWGoAMxgNKHOAvGXsIbbuBjeeORhZaIrruBwDfzgTICuwWCAhCPqMqkHrxkQMZbXUIavknNhuIycoDssXlOtbSWsxVXQhWMyDQZWDlEtewXWKBPUcHDYWWgyOerbnoAxrnpsCulOxqxdywFJFoeWNpVGIPMUJSWwvlVDWNkjIBMlXPi

It will only print

>4RYF_2
MNLIPTVIEQTSRGERAYDIYSRLLKDRIIMLGSAIDDNVANSIVSQLLFLDAQDPEKDIFLYINSPGGSISAGMAIYDTMNFVKADVQTIGMGMAASMGSFLLTAGANGKRFALPNAEIMIHQPLGGAQGQATEIEIAARHILKIKERMNTIMAEKTGQPYEVIARDTDRDNFMTAQEAKDYGLIDDIIINKSGLKGHHHHHH

Edit

The method that I would recommend to better handle edge-cases is to use purpose-built bioinformatics software, e.g. seqkit

seqkit seq -m 50 -M 800 example.fasta
>4RYF_2
MNLIPTVIEQTSRGERAYDIYSRLLKDRIIMLGSAIDDNVANSIVSQLLFLDAQDPEKDI
FLYINSPGGSISAGMAIYDTMNFVKADVQTIGMGMAASMGSFLLTAGANGKRFALPNAEI
MIHQPLGGAQGQATEIEIAARHILKIKERMNTIMAEKTGQPYEVIARDTDRDNFMTAQEA
KDYGLIDDIIINKSGLKGHHHHHH

Is perl an option?

perl -nle '$prev && print if length() >50 and length() < 800 && print $prev; $prev = $_' input_file

$prev - Create a variable which will hold every line. When the length condition is met, and there has been a previous line $prev, then it prints the condition matched in $prev and prints the last line.

$prev = $_ Assigns the current line to the prev line variable


If the upper limit 800 is not essential, could sed be an option?

$ sed -En '/>/ {N;/[a-zA-Z0-9]{50,}/p}' input_file

/>/ - Match > and read into the pattern space

N; Run the condition on the next line after the match and append that to the pattern space also:

{50,} - If the length is 50 or more

\1/p - Return it and print


Output

>4RYF_2
MNLIPTVIEQTSRGERAYDIYSRLLKDRIIMLGSAIDDNVANSIVSQLLFLDAQDPEKDIFLYINSPGGSISAGMAIYDTMNFVKADVQTIGMGMAASMGSFLLTAGANGKRFALPNAEIMIHQPLGGAQGQATEIEIAARHILKIKERMNTIMAEKTGQPYEVIARDTDRDNFMTAQEAKDYGLIDDIIINKSGLKGHHHHHH

With your shown samples, please try following awk code. Written and tested with GNU awk.

awk -v RS= '
{
  val=""
  delete arr
  while(match($0,/>[^\n]*\n*[^\n]*/)){
    val=substr($0,RSTART,RLENGTH)
    split(val,arr,"\n")
    if(length(arr[2])>50 && length(arr[2])<=800){
      print val
    }
    $0=substr($0,RSTART+RLENGTH)
  }
}
'  Input_file

If only the next line should meet the length restrictions, you can match and store the line that starts with > in a variable, for example previous

Then for the next line, check for the length and if the previous line is not empty.

If is is not, print the previous and the current line.

At the end, set the previous variable to an empty string.

awk '{
  if (/^>/) {
    previous = $0
    next
  }
  if (length(previous) != 0 && length($0) > 50 && length($0) <= 800) {
    print previous ORS $0
  }
  previous=""
}' sample.txt

See an AWK demo

Related