主页 / coding / 问题 / 77700117
Accepted
pedro
Asked: 2023-12-22 03:24:39 +0800 CST

根据条件反转特定字段中的字符串

  • 772

我有这个文件：

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

我需要tr和reverse字段 $4 et $5 仅当$2==rev

预计 ：

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

我试过 ：

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

我还尝试使用 awk （在 awk 中执行 bash 命令并打印命令输出）

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 个回答

  1. Best Answer

    要遵循问题的尝试：

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

    可以将其放在一行中（或将其复制粘贴到大多数 shell 中），为了便于阅读，我将其展开。或者，当然，将其放入程序中，特别是如果还有更多事情要做的话。

    编辑——对问题中代码的评论

    • 要将输入字符串分解为@F数组（“自动拆分”），需要-a标志

    • 由于您明确打印需要的内容，因此使用-n标志，而不是-p

    • 该行中的字段 4 和 5 是数组元素 3 和 4

    • 我假设$F[4,5]你指的是两个数组元素（应该是 3,4）。那么，应该是@F[3,4]——并且带有-w标志，对于警告，我们会听到它

    • 更重要的是，我们不能将正则表达式或tr模式绑定到列表，而只能绑定到单个标量。为了将其应用于tr多个项目，需要对其进行迭代，如上所示。

    • 元素循环块中的代码@F[3,4]在数组的每个元素上运行，并就地修改它（因为$_只是数组元素的别名），首先是tr然后通过（反转和）分配。所有这些也可以写成

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

      修饰符r使tr返回更改后的字符串，然后对其进行reverse-ed，然后通过$_别名将其分配回当前处理的数组元素

    • 3

  2. 这是在 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()

    输出：

    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. 在TXR Lisp 中，使用awk宏：

    (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. 如果你想尝试 Bipython（当你处理核酸序列时我强烈推荐），这里有一种方法：

    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

相关问题