<p>I wrote a Perl script to do this. Just like <code>zgrep</code>, it can search through as many files as you give it. You can use the <code>-c</code> switch to tell it how much context (how many characters) you want around the match. So <code>-c10</code> will search for your target string and 10 characters on either side of it. If you don't use the <code>-c</code>, it will default to 155, which should include the entire sequence for most cases. </p>

<p>To use it, save as <code>alignString.pl</code> somewhere in your <code>$PATH</code> (<code>~/bin</code> should work if you're using Ubuntu, for example), make it executable with <code>chmod a+x ~/bin/alignString</code> and then run it like this:</p>

<pre><code>alignString.pl GAAAATCAG file1.fastq.gz file2.fastq.gz file3.fastq.gz
</code></pre>

<p>For example:</p>

<p><a href="https://i.sstatic.net/hARo8.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/hARo8.png" alt="example output"></a></p>

<p>The code:</p>

<pre class="lang-perl prettyprint-override"><code>#!/usr/bin/env perl

use warnings;
use strict;
use Term::ANSIColor;
use Getopt::Std;

my %opts;
getopts('c:',\%opts);
my <span class="math-container">$context = $</span>opts{'c'} || 50;

## This subroutine does the alignment
sub align{
  ## Read the line and the target pattern
  my (<span class="math-container">$line, $</span>pat) = @_;
  ## This will find 0 to <span class="math-container">$context characters on either side
  ## of the search pattern. Note that it will find the first
  ## occurrence of the pattern.
  $</span>line =~ /(.{0,<span class="math-container">$context})$</span>pat(.{0,<span class="math-container">$context})/;
  ## Save the strings to the left and right of the pattern
  my $</span>left = <span class="math-container">$1;
  my $</span>right = <span class="math-container">$2;
  ## If the left hand string is shorter than the value of $</span>context,
  ## pad it with spaces so it can be aligned.
  while (length(<span class="math-container">$left) &lt; $</span>context) {
    <span class="math-container">$left = " " . $</span>left;
  }
  ## Return the padded left hand string, the pattenr colored red
  ## and the right hand string unchanged.
  return "<span class="math-container">$left" . colored($</span>pat,"red") . "$right";
}

## The search pattern is the 1st argument
my $pattern = shift;
## Any remaining arguments are taken as target files to search
my @fastqFiles = @ARGV;

## Run the zgrep command
open(my <span class="math-container">$grepStream, '-|', "zgrep $</span>pattern @fastqFiles");
## For each line returned by zgrep
while (my <span class="math-container">$line = &lt;$</span>grepStream&gt;) {
  ## Align this line
  <span class="math-container">$line = align($</span>line, <span class="math-container">$pattern);
  print "$</span>line\n";
}
</code></pre>

CC BY-SA 4.0
2018-10-27T20:10:33.390Z
wonderful! I've resisted learning perl for a long time but it seems useful sometimes...
conchoecia
CC BY-SA 4.0
2018-10-28T11:02:56.423Z
Beatiful! :) I didn't know about Term::ANSIColor, it is so simple to use! Your left padding while loop can probably be shortened to: `$left = " " x ($context - length($left)) if length($left) < $context`;
Peter Menzel
CC BY-SA 4.0
2018-10-28T11:07:27.640Z
@PeterMenzel thanks :) and yes, various bits of this can be golfed but I was trying to keep it as clear to non Perl people as possible.
terdon
CC BY-SA 4.0
2019-10-23T04:58:41.230Z
@terdon - over a year later I found an edge case where this doesn't work. If you search for a homopolymer like 'TTTTT', if there is a case where there are multiple matches (such as in a longer homopolymer 'TTTTTTTTTT'), this perl script will align to the rightmost match- not the leftmost
conchoecia
CC BY-SA 4.0
2018-10-27T19:24:55.807Z
2018-10-28T00:21:11.260Z
2018-10-28T00:21:11.260Z
terdon
<p>Perhaps, <code>grep</code> is not the best tool to use in this case, but it should be in principle possible by using <code>grep</code> &amp; <code>sed</code>. Here is an example showing three symbols around a match.</p>

<pre class="lang-bash prettyprint-override"><code>zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC
</code></pre>

<p>Here are some explanations:</p>

<ol>
<li><code>zcat</code> will just decompress a file.</li>
<li><code>grep</code> will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.</li>
<li>first <code>sed</code> command will add 3 symbols before all the matches.</li>
<li>second <code>sed</code> command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of <code>grep -Eo '.{3}GATCGATC.{0,3}'</code>.</li>
<li>Finally, now that all matches are aligned, <code>grep</code> will color your sequence.</li>
</ol>

<p>This code is possible to adapt to include more symbols of context.</p>

<p>I would recommend doing this in another language, though. E.g., in Python:</p>

<pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'


filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')
</code></pre>

<h1>UPD</h1>

<p>A more advanced Python version inspired by <a href="https://bioinformatics.stackexchange.com/a/5305/191">this answer</a> with support for multiple files, gzip, etc.</p>

<p><a href="https://i.sstatic.net/nHhRS.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/nHhRS.png" alt="screenshot of match script"></a></p>

<pre class="lang-python prettyprint-override"><code>#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored


def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context 
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')


# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')



if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)
</code></pre>

CC BY-SA 4.0
2018-10-27T18:23:10.967Z
2018-10-28T09:54:05.243Z
2018-10-28T09:54:05.243Z
Iakov Davydov
<p>Grepping out short sequences from a fasta or fastq file is a really useful way to look at sequencing data. Using the option <code>--color=always</code> makes this even more useful, as you can visualize where the sequences appear in sequencing reads. For example:</p>

<pre><code>zgrep 'GATCGATC' --color=always My_Hiseq_Data.fq.gz
</code></pre>

<p><a href="https://i.sstatic.net/4evAS.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4evAS.png" alt="grep color output"></a></p>

<p>To make this more useful, how can I align the output of <code>grep --color=always</code> such that the beginning of the highlighted parts are aligned?</p>

<p><a href="https://i.sstatic.net/wFRmj.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/wFRmj.png" alt="aligned output"></a></p>

<phylogenetics><fasta><fastq><illumina><bash>
How to align output of grep --color=always? (To QC fasta/fastq files)
# How to align output of grep --color=always? (To QC fasta/fastq files)

(Asked by: username_0 on 2018-10-26)

Grepping out short sequences from a fasta or fastq file is a really useful way to look at sequencing data. Using the option --color=always makes this even more useful, as you can visualize where the sequences appear in sequencing reads. For example:
zgrep 'GATCGATC' --color=always My_Hiseq_Data.fq.gz


To make this more useful, how can I align the output of grep --color=always such that the beginning of the highlighted parts are aligned?

**Tags:** phylogenetics, fasta, fastq, illumina, bash

---

## 2 answers

### Answer by username_3 on 2018-10-27. Score: 6 | Accepted answer

Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match.
zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC

Here are some explanations:

zcat will just decompress a file.
grep will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.
first sed command will add 3 symbols before all the matches.
second sed command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of grep -Eo '.{3}GATCGATC.{0,3}'.
Finally, now that all matches are aligned, grep will color your sequence.

This code is possible to adapt to include more symbols of context.
I would recommend doing this in another language, though. E.g., in Python:
#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'


filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')

UPD
A more advanced Python version inspired by this answer with support for multiple files, gzip, etc.

#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored


def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context 
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')


# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')



if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)


---

### Answer by username_1 on 2018-10-27. Score: 2

I wrote a Perl script to do this. Just like zgrep, it can search through as many files as you give it. You can use the -c switch to tell it how much context (how many characters) you want around the match. So -c10 will search for your target string and 10 characters on either side of it. If you don't use the -c, it will default to 155, which should include the entire sequence for most cases. 
To use it, save as alignString.pl somewhere in your $PATH (~/bin should work if you're using Ubuntu, for example), make it executable with chmod a+x ~/bin/alignString and then run it like this:
alignString.pl GAAAATCAG file1.fastq.gz file2.fastq.gz file3.fastq.gz

For example:

The code:
#!/usr/bin/env perl

use warnings;
use strict;
use Term::ANSIColor;
use Getopt::Std;

my %opts;
getopts('c:',\%opts);
my $context = $opts{'c'} || 50;

## This subroutine does the alignment
sub align{
  ## Read the line and the target pattern
  my ($line, $pat) = @_;
  ## This will find 0 to $context characters on either side
  ## of the search pattern. Note that it will find the first
  ## occurrence of the pattern.
  $line =~ /(.{0,$context})$pat(.{0,$context})/;
  ## Save the strings to the left and right of the pattern
  my $left = $1;
  my $right = $2;
  ## If the left hand string is shorter than the value of $context,
  ## pad it with spaces so it can be aligned.
  while (length($left) < $context) {
    $left = " " . $left;
  }
  ## Return the padded left hand string, the pattenr colored red
  ## and the right hand string unchanged.
  return "$left" . colored($pat,"red") . "$right";
}

## The search pattern is the 1st argument
my $pattern = shift;
## Any remaining arguments are taken as target files to search
my @fastqFiles = @ARGV;

## Run the zgrep command
open(my $grepStream, '-|', "zgrep $pattern @fastqFiles");
## For each line returned by zgrep
while (my $line = <$grepStream>) {
  ## Align this line
  $line = align($line, $pattern);
  print "$line\n";
}

> **Comments:**
> **username_0**: wonderful! I've resisted learning perl for a long time but it seems useful sometimes...  
> **username_2**: Beatiful! :) I didn't know about Term::ANSIColor, it is so simple to use! Your left padding while loop can probably be shortened to: `$left = " " x ($context - length($left)) if length($left) < $context`;  
> **username_1**: @PeterMenzel thanks :) and yes, various bits of this can be golfed but I was trying to keep it as clear to non Perl people as possible.  
> **username_0**: @terdon - over a year later I found an edge case where this doesn't work. If you search for a homopolymer like 'TTTTT', if there is a case where there are multiple matches (such as in a longer homopolymer 'TTTTTTTTTT'), this perl script will align to the rightmost match- not the leftmost  

---

