从字符串中间删除字符

Question

我有一个带有 RX: 字段的 SAM 文件，其中包含 12 个碱基，中间用 - 分隔，即 RX:Z:CTGTGC-TCGTAA

我想从此字段中删除连字符，但我不能简单地从整个文件中删除所有连字符，因为读取的名称包含它们，例如 1713704_EP0004-T

大部分时间都在尝试 tr, 但这只是从文件中删除所有连字符。:

tr -d '"-' < sample.fq.unaln.umi.sam > sample.fq.unaln.umi.re.sam

输入是一个超过 10,000,000 行的大型 SAM 文件，如下所示：

1902336-103-016_C1D1_1E-T:34    99  chr1    131341  36  146M    =   131376  182 GGACAGGGAGTGTTGACCCTGGGCGGCCCCCTGGAGCCACCTGCCCTGAAAGCCCAGGGCCCGCAACCCCACACACTTTGGGGCTGGTGGAACCTGGTAAAAGCTCACCTCCCACCATGGAGGAGGAGCCCTGGGCCCCTCAGGGG  NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN  MC:Z:147M   MD:Z:83T62cD:i:4    cE:f:0  PG:Z:bwa    RG:Z:A  MI:Z:34 NM:i:1  cM:i:3  MQ:i:36 UQ:i:45 AS:i:141    XS:i:136    RX:Z:CTGTGC-TCGTAA

期望的输出（即最后一个字段）

1902336-103-016_C1D1_1E-T:34    99  chr1    131341  36  146M    =   131376  182 GGACAGGGAGTGTTGACCCTGGGCGGCCCCCTGGAGCCACCTGCCCTGAAAGCCCAGGGCCCGCAACCCCACACACTTTGGGGCTGGTGGAACCTGGTAAAAGCTCACCTCCCACCATGGAGGAGGAGCCCTGGGCCCCTCAGGGG  NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN  MC:Z:147M   MD:Z:83T62cD:i:4    cE:f:0  PG:Z:bwa    RG:Z:A  MI:Z:34 NM:i:1  cM:i:3  MQ:i:36 UQ:i:45 AS:i:141    XS:i:136    RX:Z:CTGTGCTCGTAA

如何解决这个问题？

Answer 1

awk

awk '{sub(/-/,"",$NF)}1' file

正是您所需要的。

说明

从可以看出，您只关心最后一个字段。
NF 是记录包含的字段总数，因此 $NF 是最后一个字段。
sub(/-/,"",$NF) 将最后一个字段中的 - 替换为空字符串，使更改持久化。

GNU sed

出于同样的原因，

sed -Ei 's/^(.*)-//' file

会起作用。它还有一个额外的优势，那就是它可以执行就地编辑。

说明

-E 选项启用扩展正则表达式引擎。
(.*) 是一种贪婪搜索，它将匹配任何字符 (.) 任意次数 (*)。因为它是贪婪的，它会匹配任何东西直到最后一个连字符。
() 使 sed 记住匹配的内容。
在替换部分，我们只放入匹配的部分</code>（<code>1因为我们只有一对括号，注意你可以有多少个括号）没有连字符，从而有效地将它从它应该出现的最后一个字段中删除。

_{注意： GNU awk 支持 -i inplace，但我不确定从哪个版本开始。}

Answer 2

此模式在您要编辑的许多记录上，并且总是在行尾？如果是 -

sed -E 's/^(.*)(\s..:.:......)-(......\s*)$//' < sample.fq.unaln.umi.sam > sample.fq.unaln.umi.re.sam

Answer 3

最好的解决方案是使用 BAM 而不是 SAM 文件，并使用适当的 BAM parser/writer 库，例如 htslib。

如果没有，您可以通过在可选标签（第 12 列及以上）中搜索正则表达式 ^RX:Z: 来拼凑一些内容。

虽然可以使用列，但使用 sed 很难。相反，这是在 awk 中执行此操作的方法：

awk -F '[[:space:]]*' '{
    for (i = 12; i <= NF; i++) {
        if ($i ~ /^RX:Z:/) gsub("-", "", $i)
    }
}
1' file.sam

这里有一个与 Perl “one-liner”大致等效的解决方案：

perl -ape '
    for (@F[11..(scalar @F)]) {
        s/-//g if /^RX:Z:/;
    }
    $_ = join("\t", @F);
' file.sam

要在原始文件中执行替换，您可以将选项 -i.bak 传递给 perl（这将创建一个备份 file.sam.bak；如果您不想备份, 省略扩展名）。

Answer 4

我已经使用 pysam 解决了这个问题，它更快、更安全并且需要更少的磁盘 space 因为不需要 sam 文件。还不完美，还在学习中python，pysam用了半天

import pysam
import sys
from re import sub

# Provide a bam file
if len(sys.argv) == 2:
    assert sys.argv[1].endswith('.bam')

# Makes output filehandle
inbamfn = sys.argv[1]
outbamfn = sub('.bam$', '.fixRX.bam', inbamfn)

inbam = pysam.Samfile(inbamfn, 'rb')
outbam = pysam.Samfile(outbamfn, 'wb', template=inbam)

# Counters for reads processed and written
n = 0
w = 0

# .get_tag() retrieves RX tag from each read
for read in inbam.fetch(until_eof=True):
    n += 1
    umi = read.get_tag('RX')
    assert umi is not None
    umifix = umi[:6] + umi[7:]
    read.set_tag('RX', umifix, value_type='Z')
    if '-' in umifix:
        print('Hyphen found in UMI:', umifix, read)
        break
    else:
        w += 1
        outbam.write(read)

inbam.close()
outbam.close()

print ('Processed', n, 'reads:\n',
       w, 'UMIs written.\n',
       str(int((w / n) * 100)) + '% of UMIs fixed')

从字符串中间删除字符

Remove character from the middle of a string

regex

sed

bioinformatics

samtools