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
2
  • 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. Commented Jun 2, 2016 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
    Commented Jun 2, 2016 at 17:56

3 Answers 3

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

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .