3

Not sure how many people are familiar with DNA sequencing data, but if this is part of my file (lines starting with ">" are IDs, and lines starting with a letter are DNA sequences):

>NB501013:9:HJJ75BGXX:4:13609:24076:18015/2
GGGGGGGAAAAAAA
>NB501013:9:HJJ75BGXX:4:21602:19346:16945/2
CTCGTCGCATCACAAAGGGAT
>NB501013:9:HJJ75BGXX:3:11407:17650:13229/2
CCGCGGGCCGGTGCGGGGGTTTTTTTGTTTTTTTGGTTACAACGGGTGGG
>NB501013:9:HJJ75BGXX:3:13509:1817:13239/2
CAGCCC
>NB501013:9:HJJ75BGXX:4:22611:20567:13384/2
GAATA

I would want to remove the line: GGGGGGGAAAAAAA

Along with its sequencing ID (I know you can do that using grep -B1). But does anyone know how to remove the lines that are only composed up of two letters?

Also, for sequences that are shorter than 5 letters, I would want to remove those along with their IDs, I can't simply grep for lines longer than a certain length because all the IDs are pretty long, so I need to somehow use grep -v on lines that start with a letter (so doesn't start with ">") and longer than a certain length.

Therefore, my sample output would be:

>NB501013:9:HJJ75BGXX:4:21602:19346:16945/2
CTCGTCGCATCACAAAGGGAT
>NB501013:9:HJJ75BGXX:3:11407:17650:13229/2
CCGCGGGCCGGTGCGGGGGTTTTTTTGTTTTTTTGGTTACAACGGGTGGG
>NB501013:9:HJJ75BGXX:3:13509:1817:13239/2
CAGCCC
  • It would have been thoughtful to provide data where the DNA sequence is actually, you know, only 2 letters. Like what you ask for in the title of your question. – glenn jackman Jun 2 '16 at 17:42
  • Hello, sorry this is my first time posting on here, I didn't realize the protocol. But I'm a little confused, I did provide a DNA sequence that is made up of only two letters (the first DNA sequence). I specifically didn't post a sequence that is only GA because I wanted to remove a sequence composed of two letters but those two letters can appear more than once. Also, I did include an input sample and described my output, so I'm not sure what I should do differently next time. Please let me know if I did anything wrong! – ehze Jun 2 '16 at 17:56
3

Give grep with Perl Compatible REgexp module a try:

  • to remove two-letters combinations:

    pcregrep -Mv '>.*\n([ACGT])\1*([ACGT])\2*(\1|\2)*$' file
    

    output:

    >NB501013:9:HJJ75BGXX:4:21602:19346:16945/2
    CTCGTCGCATCACAAAGGGAT
    >NB501013:9:HJJ75BGXX:3:11407:17650:13229/2
    CCGCGGGCCGGTGCGGGGGTTTTTTTGTTTTTTTGGTTACAACGGGTGGG
    >NB501013:9:HJJ75BGXX:3:13509:1817:13239/2
    CAGCCC
    >NB501013:9:HJJ75BGXX:4:22611:20567:13384/2
    GAATA
    
  • to remove combination of 5-letters or less:

     pcregrep -Mv '>.*\n[ACGT]{1,5}$' file
    

    output:

    >NB501013:9:HJJ75BGXX:4:13609:24076:18015/2
    GGGGGGGAAAAAAA
    >NB501013:9:HJJ75BGXX:4:21602:19346:16945/2
    CTCGTCGCATCACAAAGGGAT
    >NB501013:9:HJJ75BGXX:3:11407:17650:13229/2
    CCGCGGGCCGGTGCGGGGGTTTTTTTGTTTTTTTGGTTACAACGGGTGGG
    >NB501013:9:HJJ75BGXX:3:13509:1817:13239/2
    CAGCCC
    
1
#!/usr/bin/env perl
#
# Usage: thisscriptname < someinputfile

use strict;
use warnings;

while (1) {
  exit if eof;
  # rash assumption there are always pairs of ID and sequence lines
  # NOTE these contain a newline, so many need chomp() depending
  # on what you do with them...
  my $id = readline;
  my $seq = readline;

  # calculate unique sequence letters via hash (is there also a U
  # or something? been a few decades since AP bio...)
  my %chars;
  $chars{$_}++ for $seq =~ m/[ATGC]/g;

  # business logic time!
  if (keys %chars > 2 and length $seq >= 5) {
    print $id;
    print $seq;
  }
}
1

You might consider reversing the file, testing the DNA sequence, and if the test passes, ignore this line and the next line:

tac file |
  awk '!/^>/ && (length($1) < 5 || $1 == "GGGGGGGAAAAAAA") {getline; next} 1' |
  tac

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.