☰

sed: busca patrón y reemplaza por otro patrón en la misma línea

tag-icon • tag-icon sed bioinformatics
sed: busca patrón y reemplaza por otro patrón en la misma línea

Tengo un archivo con gene_id y nombres de genes en una línea. Quiero reemplazar la palabra despuésid_gencon la palabra despuésgeneo despuésproductoo despuésbrotar(si algo se perdió).

A continuación se muestra un ejemplo de una línea:

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";

Intenté hacerlo con sed:

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

Pero los resultados fueron incorrectos:

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] .*?;

Pero quiero guardar todas las líneas, pero con otra palabra después.id_gen, como esto:

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";

O así (si se perdió otro):

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";

Cualquier ayuda será muy apreciada.

Respuesta1

El siguiente script en Perl intenta hacer coincidir gene, producty sproten cada línea de entrada, en ese orden (es decir, prioriza el gen sobre el producto y el producto sobre sprot). Si uno de ellos coincide, extrae la palabra después de la coincidencia. Se supone que la palabra está entre comillas dobles.

Si se encuentra una coincidencia, reemplaza la palabra siguiente gene_idcon la palabra extraída.

La línea se imprime independientemente de que haya sido modificada o no.

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

Alternativamente, esto podría escribirse para usar un bucle para iterar sobre las palabras clave coincidentes:

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

En mi opinión, esta versión es mejor porque es más fácil de leer y comprender (en particular, el foreachbucle enfatiza que se trata de iterar sobre una lista de palabras). Más importante aún, evita repetir la $word = $1declaración: es menos probable que cometa un error si necesita cambiarla o agregar código adicional, si solo tiene que hacerlo una vez en lugar de tres. "No te repitas" no es tan importante en un programa pequeño y trivial como este, pero puede ser muy importante en programas más grandes. De todos modos, evitar/minimizar la repetición es un buen hábito de programación.

Si el orden de la coincidencia no fuera significativo (es decir, si no le importaba cuál se encontró, siempre y cuando lo fuera), entonces podría simplificar el script:

#!/usr/bin/perl

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

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

  print;
}

Independientemente de la versión del script que utilice, guárdelo como, por ejemplo replace.pl, y hágalo ejecutable con chmod +x replace.pl. O pruébalos todos como replace1.pl, replace2.pl, replace3.pl. Luego ejecútelo así:

$ ./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";

Respuesta2

Explotamos la propiedad de un hash de que si se aplican varios valores a una clave determinada, el último se convierte en el valor final.

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

Dado que los comandos son todos iguales y solo cambian los nombres, podemos fácilmente hacer que toda la tabla de operaciones sea controlada, en la que solo proporcionamos los nombres de los campos mientras que el bucle for se encargará del resto.

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

Respuesta3

Para completar la perlsolución, así es como lo harías con sed. No estoy seguro de cómo espera que funcione la sintaxis dada, pero en realidad necesita una expresión regular que coincida con la cadena.

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

gene_idy genese combinan por sí mismos. Una cadena entre comillas dobles es una concatanación de una comilla doble, cualquier número de caracteres que no sean comillas dobles ( [^"]*) y otra comilla doble. Finalmente tienes las cosas en el medio..*

Ahora necesitas colocar \(\)las piezas que necesitas reciclar en el reemplazo:

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

El par exterior cubre todo lo que debe dejarse intacto. Este se reutiliza como \1en el reemplazo. El par interno es la cadena que desea reutilizar gene_id.

Ahora, si desea tener producto sprotcomo reemplazos alternativos, puede usar cadenas alternativas de expresiones regulares extendidas:

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

pero este no preferirá genesobre productover sprot, sino que preferirá el último de ellos que esté presente. Si desea tener ese orden de preferencia, necesita pasos separados y comenzar con el último, para que pueda ser reemplazado por uno mejor:

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/'

O, si se sabe que el orden de geney productsprot` es fijo, primero puede extraer la identificación preferida mientras estaciona la línea real en el espacio de espera:

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

El marcador #puede ser cualquier cadena que se sepa que no forma parte del ID; para GNU sedpodrías usarlo \npara estar seguro. Entonces reemplaza la primera de dichas cadenas por el marcador y elimina el resto de la línea, luego elimina todo hasta el marcador, por lo que ahora solo queda la ID en el espacio del patrón. Luego, se agregará Gla línea original (que conservamos en el búfer de retención con ) y luego la ID (la parte antes de la nueva línea) reemplaza la posterior . De alguna manera es más fácil escribir que explicar.h"string"gene_id

información relacionada