1

我有一个带有 RX: 字段的 SAM 文件，其中包含 12 个碱基，中间用-ie分隔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

我该如何解决这个问题？

4

4 回答 4

4

awk

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

是你需要的。

解释

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

GNU sed

出于同样的原因，

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

将工作。它还有一个额外的优势，即它可以执行就地编辑。

解释

  • 该-E选项启用扩展的正则表达式引擎。
  • 这(.*)是一个贪婪的搜索，它将匹配任何字符（.）任意次数（*）。对于贪婪的事实，它将匹配到最后一个连字符的任何内容。
  • 使记住匹配()的内容。sed
  • 在替换部分，我们只放置了匹配的部分\1（1因为我们只有一对括号，请注意，您可以拥有任意数量的连字符），从而有效地将其从应该出现的最后一个字段中删除。

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

于 2019-05-01T15:06:06.523 回答
2

我已经使用 pysam 解决了这个问题，它更快、更安全并且需要更少的磁盘空间，因为不需要 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')
于 2019-05-10T12:55:29.987 回答
1

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

^RX:Z:缺少它，您可以通过在可选标签（第 12 列及以上）中搜索正则表达式来拼凑一些东西。

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

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

这是一个与 Perl“单线”大致等效的解决方案：

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

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

于 2019-05-09T10:37:00.673 回答
0

此模式存在于您要编辑的许多记录上，并且始终位于行尾？如果是这样的话 -

sed -E 's/^(.*)(\s..:.:......)-(......\s*)$/\1\2\3/' < sample.fq.unaln.umi.sam > sample.fq.unaln.umi.re.sam
于 2019-05-01T15:07:50.450 回答