Início / coding / Perguntas / 77700117
Accepted
pedro
Asked: 2023-12-22 03:24:39 +0800 CST

String reversa em campos específicos com condição

  • 772

Eu tenho esse arquivo:

m64071_220512_054244/12584899/ccs rev pet047-10055 ACGTGCGACCTTGTGA TTGAGGGTTCAAACGTGCGACCTTGTGA
m64071_220512_054244/128321000/ccs rev pet047-10055 ACGTGCGACCTTGTGA TTGAGGGTTCAAACGTGCGACCTTGTGA
m64071_220512_054244/132186699/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/134874748/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA

Eu preciso tre reverseos campos $4 e $5 somente se$2==rev

Esperar :

m64071_220512_054244/12584899/ccs rev pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/128321000/ccs rev pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/132186699/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/134874748/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA

Tentei :

perl -lpe 'if(/rev/) {$rev=/rev/;next}; if ($rev) {$F[4,5]=~tr/ATGC/TACG/; $F[4,5]=reverse $F[4,5]; print "@F"}' file

Eu também tentei usar o awk ( Execute o comando bash dentro do awk e imprima a saída do comando )

awk '{
            if($2==rev)
        {
            cmd1="echo \047" $4 "\047 | rev | tr \047ATGC\047 \047TACG\047" 
            cmd2="echo \047" $5 "\047 | rev | tr \047ATGC\047 \047TACG\047"
            newVar1=((cmd1 | getline line) > 0 ? line : "failed") 
            newVar2=((cmd2 | getline line) > 0 ? line : "failed")
            close(cmd)
            print $1, $2, $3, newVar1, newVar2
        }
        else {print}
}' file
bash

4 respostas

  1. Best Answer

    Para seguir a tentativa da pergunta:

    perl -w -lanE'
    if ($F[1] eq "rev") { 
        for (@F[3,4]) { tr/ATGC/TACG/; $_ = reverse $_ } 
    } 
    say "@F"
' file

    Posso colocar isso em uma linha (ou copiar e colar como está na maioria dos shells), espalhei-o para facilitar a leitura. Ou coloque em um programa, claro, principalmente se houver mais o que fazer.

    Editar - comentários sobre o código da pergunta

    • Para dividir a string de entrada no @Farray (" autosplit "), é necessário -asinalizar

    • Como você está imprimindo explicitamente o que é necessário, use -no sinalizador, não-p

    • Os campos 4 e 5 na linha são os elementos 3 e 4 da matriz

    • Presumo que $F[4,5]você queira dizer os dois elementos da matriz (que devem ser 3,4). Isso, então, deveria ser @F[3,4]- e com -wbandeira, como aviso, ouviríamos sobre isso

    • Mais importante ainda, não podemos vincular uma expressão regular ou trpadrão a uma lista, mas apenas a um único escalar. Para aplicar isso tra vários itens, é necessário iterar sobre eles, como acima.

    • O código no bloco desse loop sobre os @F[3,4]elementos é executado em cada elemento da matriz e o modifica no local (já que $_é apenas um alias para os elementos da matriz), primeiro com tr(revertendo e) atribuindo. Tudo o que também pode ser escrito como

      $_ = reverse tr/ATGC/TACG/r  for @F[3,4];

      O rmodificador trretorna a string alterada, que é então reverse-ed e atribuída de volta ao elemento da matriz atualmente processado, por meio do $_alias para ela

    • 3

  2. Aqui está uma maneira de fazer isso em Python:

    #!/usr/bin/env python

import io
import sys

RECORDS_STR = '''m64071_220512_054244/12584899/ccs rev pet047-10055 ACGTGCGACCTTGTGA TTGAGGGTTCAAACGTGCGACCTTGTGA
m64071_220512_054244/128321000/ccs rev pet047-10055 ACGTGCGACCTTGTGA TTGAGGGTTCAAACGTGCGACCTTGTGA
m64071_220512_054244/132186699/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/134874748/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA'''

'''
fast pure-Python reverse complement, courtesy of Devon Ryan
ref. https://bioinformatics.stackexchange.com/a/3585/776
'''
DNA_TABLE = str.maketrans("ACTGactg", "TGACtgac")
def reverse_complement(seq):
    return seq.translate(DNA_TABLE)[::-1]

def main():
    records = io.StringIO(RECORDS_STR) # replace with sys.stdin etc.
    for line in records:
        elems = line.rstrip().split()
        if elems[1] == 'rev':
            elems[3] = reverse_complement(elems[3])
            elems[4] = reverse_complement(elems[4])
        sys.stdout.write('{}\n'.format('\t'.join(elems)))

if __name__ == "__main__":
    main()

    Saída:

    m64071_220512_054244/12584899/ccs   rev pet047-10055    TCACAAGGTCGCACGT    TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/128321000/ccs  rev pet047-10055    TCACAAGGTCGCACGT    TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/132186699/ccs  fwd pet047-10055    TCACAAGGTCGCACGT    TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/134874748/ccs  fwd pet047-10055    TCACAAGGTCGCACGT    TCACAAGGTCGCACGTTTGAACCCTCAA
    • 2

  3. No TXR Lisp, usando a awkmacro:

    (awk ((equal [f 1] "rev")
      (each ((i 3..5))
        (upd [f i] (mapcar (relate "ATGC" "TACG")) reverse)))
     (t))

    txr rev.tl data
m64071_220512_054244/12584899/ccs rev pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/128321000/ccs rev pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/132186699/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA
m64071_220512_054244/134874748/ccs fwd pet047-10055 TCACAAGGTCGCACGT TCACAAGGTCGCACGTTTGAACCCTCAA

$ diff <(txr rev.tl data) expected
$
    • 1

  4. Se você quiser experimentar o Bipython (eu recomendo fortemente quando você trabalha com sequências de ácidos nucléicos), aqui está uma maneira de fazer isso:

    from Bio.Seq import Seq
import re

with open("input_file.txt") as f:
    for line in f:
        line = line.rstrip() # remove new line
        fields = re.split('\s', line)
        if fields[1] != 'rev':
            print(line)
        else:
            dna = Seq(fields[3])
            print(' '.join(fields[0:3]), dna.reverse_complement(), fields[4])
    • 0

relate perguntas