根据另一个文件的列中的值过滤文件

Question

我有一个包含职位列表（第 1 + 2 列）和与这些职位相关的值的文件：

文件1.txt：

1 20   A   G
4 400  T   C
1 12   A   T
2 500  G   C

还有一些位置相同的另一个文件。可能有多行与 File1.txt 中的位置相同

文件2.txt

#CHR    POS       Count_A Count_C Count_G Count_T
1       20        0       18      2       0
4       400       0       0       0       1
1       12        0       7       0       40
4       400       0       1       0       1
5       50        16      0       0       0
2       500       9       0       4       0

我需要输出一个 File1.txt 版本，不包括任何同时满足这两个条件的行：

1：如果位置（第 1+2 列）在 File1.txt 和 File2.txt 中匹配 2：如果 File2.txt 中与该位置的 File1.txt 第 4 列中的字母(A,G,C,T) 匹配的列中的计数 > 0。

因此对于上面的示例，不会输出 File1.txt 的第一行，因为在 file2.txt 中匹配的行（基于前 2 列：1 20），第 4 列有字母 G，为此File2.txt 中的行 Count_G 列 > 0。

此示例将输出的唯一行是：

2 500 G C

对我来说，特别棘手的部分是 file2.txt 中可以有多个匹配的行，如果 File2.txt 中的相应列在 File2 中的一行中 >0，我想排除 File1.txt 中的行。文本。这意味着在上面的示例中，不会包含 File1.txt 的第 2 行，因为 Count_C 在 File2.txt 中第二次出现时 > 0 (Count_C = 1)。

我不确定这种过滤是否可以一步完成。在 File1.txt 中输出行列表是否更容易，其中 File2.txt 中 File1.txt 中第 4 列中的字母的计数> 0。然后使用此列表与 File1.txt 进行比较并删除两个文件中出现的任何行？

我之前使用下面的代码根据另一个文件中的值过滤了一个文件，但这是针对 file2.txt 中只有一列值要过滤的情况。我不确定如何进行条件过滤，以便根据 file1.txt 的第 4 列中的字母检查右列

我当前的代码是在 python 中，但欢迎任何解决方案：

f2 = open('file2.txt', 'r')
d2 = {}
for line in f2.split('\n'):
    line = line.rstrip()
    fields = line.split("\t")
    key = (fields[0], fields[1])
    d2[key] = int(fields[2])

f1 = open('file1.txt', 'r')
for line in file1.split('\n'):
    line = line.rstrip()
    fields = line.split("\t")
    key = (fields[0], fields[1])
    if d2[key] > 1000:
        print line

我认为我之前的解决方案已经非常冗长，并且觉得可能有一个简单的工具可以解决我不知道的这类问题。

Answer 1

我使用 Perl 来解决这个问题。首先，它将 File2 加载到由 char、pos 和核苷酸作为键的哈希表中，值是核苷酸的编号。然后，处理第二个文件。如果哈希表中的字符、位置和核苷酸的值不为零，则不会打印。

#!/usr/bin/perl
use warnings;
use strict;

my %gt0;
open my $in2, '<', 'File2.txt' or die $!;
<$in2>;  # Skip the header.
while (<$in2>) {
    my %count;
    (my ($chr, $pos), @count{qw{ A C G T }}) = split;
    $gt0{$chr}{$pos}{$_} = $count{$_} for qw( A C G T );
}

open my $in1, '<', 'File1.txt' or die $!;
while (<$in1>) {
    my ($chr, $pos, undef, $c4) = split;
    print unless $gt0{$chr}{$pos}{$c4};
}

Answer 2

您的代码对我来说似乎很不错。您也许可以编辑

d2[key] = int(fields[2])

和

if d2[key] > 1000:
        print line

他们让我有点困惑。

我会这样做：

f2 = open('file2.txt', 'r')
d2 = {}
for line in f2.split('\n'):
    fields = line.rstrip().split("\t")
    key = (fields[0], fields[1])
    d2[key] = {'A':int(fields[2]),'C':int(fields[3]),'G':int(fields[4]),
      'T':int(fields[5])}

f1 = open('file1.txt', 'r')
for line in f1.split('\n'):
    line = line.rstrip()
    fields = line.split("\t")
    key = (fields[0], fields[1])
    if (key not in d2) or (d2[key][str(fields[2])] == 0 and d2[key][str(fields[3])] == 0):
        print(line)

编辑： 如果您有任意数量的字母（和文件 2 中的列），只需概括我硬编码的 d2 中的字典。简单的。让我们添加 2 个字母：

col_names = ['A','C','G','T','K','L']

for a,i in zip(fields[2:],range(len(fields[2:]))):
    d2[key][col_names.index(i)] = a