☰

sed: найти шаблон и заменить его на другой шаблон в той же строке

tag-icon • tag-icon sed bioinformatics
sed: найти шаблон и заменить его на другой шаблон в той же строке

У меня есть файл с gene_id и именами генов в одной строке. Я хочу заменить слово послеген_idсо словом послегенили послепродуктили послеспрот(если что-то из этого пропущено).

Вот пример строки:

chrM    Gnomon  CDS 8345    8513    .   +   1   gene_id "cds-XP_008824843.3"; transcript_id "cds-XP_008824843.3"; Parent "rna-XM_008826621.3"; Dbxref "GeneID:103728653_Genbank:XP_008824843.3"; Name "XP_008824843.3"; end_range "8513,."; gbkey "CDS"; gene "semaphorin-3F"; partial "true"; product "semaphorin-3F"; protein_id "XP_008824843.3"; sprot "sp|Q13275|SEM3F_HUMAN";
chrM    StringTie   exon    2754    3700    .   +   .   gene_id "cds-YP_007626758.1"; transcript_id "cds-YP_007626758.1"; Parent "gene-ND1"; Dbxref "Genbank:YP_007626758.1,Gene "ID:15088436"; Name "YP_007626758.1"; Note "TAAstopcodoniscompletedbytheadditionof3'AresiduestothemRNA"; gbkey "CDS"; gene "ND1"; product "NADHdehydrogenasesubunit1"; protein_id "YP_007626758.1"; transl_except "(pos:3700..3700%2Caa:TERM)"; transl_table "2";

Я попробовал сделать это с помощью sed:

sed -E 's/[^gene_id] .*?;/[^gene] .*?;|[^sprot] .*?;|[^product] .*?;/g'

Но результаты оказались неверными:

chrM    Gnomon  CDS 8345    8513    .   +   1   gene_id "cds-XP_008824843.3"[^gene] .*?;|[^sprot] .*?;|[^product] .*?;
chrM     StringTie       exon    2754    3700    .       +       .       gene_id "cds-YP_007626758.1"[^gene] .*?;|[^sprot] .*?;|[^product] .*?;

Но я хочу сохранить всю строку, но с другим словом послеген_id, так:

chrM    Gnomon  CDS 8345    8513    .   +   1   gene_id "semaphorin-3F"; transcript_id "cds-XP_008824843.3"; Parent "rna-XM_008826621.3"; Dbxref "GeneID:103728653_Genbank:XP_008824843.3"; Name "XP_008824843.3"; end_range "8513,."; gbkey "CDS"; gene "semaphorin-3F"; partial "true"; product "semaphorin-3F"; protein_id "XP_008824843.3"; sprot "sp|Q13275|SEM3F_HUMAN";
chrM     StringTie       exon    2754    3700    .       +       .       gene_id "ND1"; transcript_id "cds-YP_007626758.1"; Parent "gene-ND1"; Dbxref "Genbank:YP_007626758.1,Gene "ID:15088436"; Name "YP_007626758.1"; Note "TAAstopcodoniscompletedbytheadditionof3'AresiduestothemRNA"; gbkey "CDS"; gene "ND1"; product "NADHdehydrogenasesubunit1"; protein_id "YP_007626758.1"; transl_except "(pos:3700..3700%2Caa:TERM)"; transl_table "2";

Или вот так (если кто-то пропустил):

chrM    Gnomon  CDS 8345    8513    .   +   1   gene_id "sp|Q13275|SEM3F_HUMAN"; transcript_id "cds-XP_008824843.3"; Parent "rna-XM_008826621.3"; Dbxref "GeneID:103728653_Genbank:XP_008824843.3"; Name "XP_008824843.3"; end_range "8513,."; gbkey "CDS"; gene "semaphorin-3F"; partial "true"; product "semaphorin-3F"; protein_id "XP_008824843.3"; sprot "sp|Q13275|SEM3F_HUMAN";
chrM     StringTie       exon    2754    3700    .       +       .       gene_id "ND1"; transcript_id "cds-YP_007626758.1"; Parent "gene-ND1"; Dbxref "Genbank:YP_007626758.1,Gene "ID:15088436"; Name "YP_007626758.1"; Note "TAAstopcodoniscompletedbytheadditionof3'AresiduestothemRNA"; gbkey "CDS"; gene "ND1"; product "NADHdehydrogenasesubunit1"; protein_id "YP_007626758.1"; transl_except "(pos:3700..3700%2Caa:TERM)"; transl_table "2";

Любая помощь будет очень признательна.

решение1

Следующий скрипт perl пытается сопоставить gene, product, и sprotв каждой входной строке в указанном порядке (т. е. он отдает приоритет gene над product и product над sprot). Если один из них совпадает, он извлекает слово после совпадения. Предполагается, что слово заключено в двойные кавычки.

Если совпадение найдено, то слово после него заменяется gene_idизвлеченным словом.

Строка печатается независимо от того, была ли она изменена или нет.

#!/usr/bin/perl

while (<>) {
  my $word = '';

  if (m/\b(?:gene)\s+("[^"]*")/) {
    $word = $1;
  } elsif (m/\b(?:product)\s+("[^"]*")/) {
    $word = $1;
  } elsif (m/\b(?:sprot)\s+("[^"]*")/) {
    $word = $1;
  };

  if ($word) {
    s/\bgene_id\s+(?:"[^"]*")/gene_id $word/
  };

  print;
}

В качестве альтернативы это можно было бы записать так, чтобы использовать цикл для перебора соответствующих ключевых слов:

#!/usr/bin/perl

while (<>) {
  my $word = '';

  foreach my $match (qw(gene product sprot)) {
    if (m/\b(?:$match)\s+("[^"]*")/) {
      $word = $1;
      last; # first match wins, exit this loop
    }
  };

  if ($word) {
    s/\bgene_id\s+(?:"[^"]*")/gene_id $word/
  };

  print;
}

На мой взгляд, эта версия лучше, потому что ее легче читать и понимать (в частности, цикл foreachподчеркивает, что речь идет об итерации по списку слов). Что еще важнее, она избегает повторения оператора $word = $1— вы с меньшей вероятностью сделаете ошибку, если вам нужно будет изменить его или добавить дополнительный код, если вам нужно сделать это только один раз, а не трижды. «Не повторяйтесь» не так уж важно в такой тривиальной маленькой программе, как эта, но может быть очень важно в более крупных программах. В любом случае, избегать/минимизировать повторения — это хорошая привычка программирования.

Если порядок сопоставления не имеет значения (т.е. если вам все равно, какой из них будет найден, главное, чтобы он был найден), то вы можете упростить скрипт:

#!/usr/bin/perl

while (<>) {
  my ($word) = m/\b(?:gene|product|sprot)\s+("[^"]*")/;

  if ($word) {
    s/\bgene_id\s+(?:"[^"]*")/gene_id $word/
  };

  print;
}

Независимо от того, какую версию скрипта вы используете, сохраните его как, например replace.pl, и сделайте его исполняемым с помощью chmod +x replace.pl. Или попробуйте их все как replace1.pl, replace2.pl, replace3.pl. Затем запустите его так:

$ ./replace.pl input.txt 
chrM    Gnomon  CDS 8345    8513    .   +   1   gene_id "semaphorin-3F"; transcript_id "cds-XP_008824843.3"; Parent "rna-XM_008826621.3"; Dbxref "GeneID:103728653_Genbank:XP_008824843.3"; Name "XP_008824843.3"; end_range "8513,."; gbkey "CDS"; gene "semaphorin-3F"; partial "true"; product "semaphorin-3F"; protein_id "XP_008824843.3"; sprot "sp|Q13275|SEM3F_HUMAN";
chrM    StringTie   exon    2754    3700    .   +   .   gene_id "ND1"; transcript_id "cds-YP_007626758.1"; Parent "gene-ND1"; Dbxref "Genbank:YP_007626758.1,Gene "ID:15088436"; Name "YP_007626758.1"; Note "TAAstopcodoniscompletedbytheadditionof3'AresiduestothemRNA"; gbkey "CDS"; gene "ND1"; product "NADHdehydrogenasesubunit1"; protein_id "YP_007626758.1"; transl_except "(pos:3700..3700%2Caa:TERM)"; transl_table "2";

решение2

Мы используем свойство хеша, заключающееся в том, что если к данному ключу применено несколько значений, то последнее становится окончательным значением.

perl -lpe 'my($l,%h)=($_);
  $h{gene_id}=$_ for map {
     $l =~ /\b$_\s+(".*?");/
  } reverse qw(gene product sprot);
  s/\bgene_id\s+\K".*?";/$h{gene_id};/;
' your_file_genes

Поскольку все команды одинаковы и меняются только имена, мы можем легко сделать всю таблицу операций управляемой, просто указав имена полей, а цикл for позаботится обо всем остальном.

for i in gene product sprot;do
  cat - <<\_FMT_ |\
  sed -e "s/%s/$i/"
s/(\<gene_id\s+)"[^"]*"(.*\s%s\s+("[^"]*"))/\1\3\2/;t
_FMT_
done | sed -Ef - your_file_genes

решение3

Чтобы завершить perlрешение, вот как это сделать с sed. Я не уверен, как вы ожидаете, что ваш синтаксис будет работать, но на самом деле вам нужно регулярное выражение для сопоставления строки

... gene_id "remove me" ... some other stuff gene "replacement" ... more stuff
    =======                                  ====
    gene_id   "[^"]*"        .*              gene    "[^"]*"

gene_idи geneсопоставляются сами по себе. Строка в двойных кавычках — это конкатанация двойных кавычек, любого количества символов, которые не являются двойными кавычками ( [^"]*) и еще одной двойной кавычки. Наконец, у вас есть то, что находится между.*

Теперь вам нужно разложить \(\)детали, которые необходимо переработать при замене:

sed 's/gene_id "[^"]*"\(.* gene \("[^"]*"\)\)/gene_id \2\1/'

Внешняя пара охватывает все, что должно остаться нетронутым. Это используется повторно, как \1при замене. Внутренняя пара — это строка, которую вы хотите использовать повторно как gene_id.

Теперь, если вы хотите иметь productили sprotв качестве альтернативных замен, вы можете использовать альтернативные строки расширенных регулярных выражений:

sed -E 's/gene_id "[^"]*"(.*(gene|product|sprot) ("[^"]*"))/gene_id \3\1/'

но это не будет предпочитать geneover productover sprot, а предпочтет последний из них, который присутствует. Если вы хотите иметь этот порядок предпочтений, вам нужны отдельные шаги и начать с последнего, поэтому его можно заменить лучшим:

sed 's/gene_id "[^"]*"\(.* sprot \("[^"]*"\)\)/gene_id \2\1/
     s/gene_id "[^"]*"\(.* product \("[^"]*"\)\)/gene_id \2\1/
     s/gene_id "[^"]*"\(.* gene \("[^"]*"\)\)/gene_id \2\1/'

Или, если известно, что порядок gene, productи sprot` фиксирован, вы можете сначала извлечь предпочтительный идентификатор, одновременно помещая фактическую линию в пространство ожидания:

sed -E 'h;s/(sprot|product|gene) ("[^"]*").*/#\2/;s/.*#//;G;s/(.*)\n(.*gene_id )"[^"]*"/\2\1/'

Маркером #может быть любая строка, которая, как известно, не является частью идентификатора; для GNU sedвы можете использовать , \nчтобы быть уверенным. Поэтому вы заменяете первую из указанных строк маркером и удаляете остальную часть строки, затем удаляете все до маркера, так что теперь в пространстве шаблона остается только идентификатор. Затем, с Gисходной строкой (которую мы сохранили в буфере удержания с помощью h) будет добавлена, а затем идентификатор (часть до новой строки) заменяет "string"после gene_id. Как-то проще написать, чем объяснить.

Связанный контент