☰

2文字のみで構成された行と1文字で始まり、特定の長さを満たす行の場合はgrep / awk / sed

tag-icon • tag-icon linux text-processing awk sed grep
2文字のみで構成された行と1文字で始まり、特定の長さを満たす行の場合はgrep / awk / sed

DNA塩基配列分析データについて知っている人がどれくらいになるのかわかりませんが、これが私のファイルの一部である場合（「>」で始まる行はIDで、文字で始まる行はDNA塩基配列です）：

>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

この行を削除したい：GGGGGGGAAAAAAA

シーケンスIDと一緒に（これを行うために使用できることを知っていますgrep -B1）。しかし、2文字だけで構成された行を削除する方法を知っている人はいますか？

また、5文字より短いシーケンスの場合は、そのIDとそのIDを削除したいと思います。すべてのIDが長いため、特定の長さよりも長い行を単にgrepすることはできませんgrep -v。文字（したがって「>」で始まらない行）であり、特定の長さよりも長くなります。

したがって、私の出力例は次のようになります。

>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

erl互換のgexpモジュールをgrep試してみてください：PCRE

  • 2文字の組み合わせを削除する：

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

    出力：

    >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

  • 5文字以下の組み合わせを削除してください。

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

    出力：

    >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

答え2

#!/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;
  }
}

答え3

ファイルを裏返してDNA配列をテストすることを検討することができます。テストに合格すると、この行を無視してNextワイヤー:

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

関連情報