如何根据2个输入文件输出文件避免MemoryErrors python

Question

我有两个以制表符分隔的大文件，如下所示：

第一个文件：inputFile1.txt

#CHR    #POS    #ID    #REF    #ALT    #IND1    #IND2    #IND3    #IND4
chr1    1000    .      A       C       0/1:3:2  0/0:2:8  0/1:8:6  1/1:4:1
chr1    1050    .      G       A       0/0:2:8  0/1:1:2  1/1:4:1  0/1:8:6
[...]

第二个文件：inputFile2.txt

#CHR    #POS    #REF    #ALT    #IND5    #IND6    #IND7    #IND8    #IND9
chr1    1000    A       T       0/1      0/0      0/1      1/1      0/1
chr1    2000    T       A       0/0      0/1      1/1      0/1      1/1
[...]

两个文件包含相同的染色体 (#CHR)，但位置更多 (#POS) 和不同的个体 (#IND)。

我的目的是根据上述文件创建一个新文件。此输出文件如下所示：

所需的输出文件：outputFile.txt

#CHR    #POS    #REF    #ALT    #IND1    #IND2    #IND3    #IND4    #IND5    #IND6    #IND7    #IND8    #IND9
chr1    1000    A       C,T     0/1      0/0      0/1      1/1      0/2      0/0      0/2      2/2      0/2 
chr1    1050    G       A       0/0      0/1      1/1      0/1      0/0      0/0      0/0      0/0      0/0     
chr1    2000    T       A       0/0      0/0      0/0      0/0      0/0      0/1      1/1      0/1      1/1
[...]

基本上，此输出文件与 inputFile2.txt 具有相同的格式，包括两个文件的所有个体 (#IND)。它还包含两个文件的所有位置，但每个 #IND 的信息可能会根据几个条件而有所不同：

1) 如果#POS 仅存在于一个inputFiles 中，则在缺少的#IND 的信息中添加0/0

2）如果两个inputFiles中都存在#POS，则根据几个条件修改#IND的内容，为简单起见，这里不再赘述。

我已经使用字典编写了一个脚本，它适用于小文件。但是当我处理大文件时，会出现 MemoryErrors

我的下一次尝试是创建以下代码结构，但失败了。

import sys;

outputFile = open("outputFile.txt", 'w')

with open("inputFile1.txt", 'r') as f:
    for line in csv.reader(f, delimiter="\t"):
        if line[0].startswith('#'): ## header
            with open("inputFile2.txt", 'r') as f2:
                for line2 in csv.reader(f2, delimiter="\t"):    
                    if line2[0].startswith('#'):  ## header
                        outputFile.write("\t".join(line2))
                        for index in range(len(line)):
                            if index > 4:
                                outputFile.write("\t"+line[index])
                        outputFile.write("\n")
                    else:
                        continue
        else:  ## positions from inputFile1

在这个“其他”中，我想检查从 inputFile1.txt 读取的特定 #POS 是否在 inputFile2.txt 的 #POS 列中（但没有创建任何列表或字典，因为它会出现 MemoryErrors） .我已经尝试使用循环读取 else 中的 inputFile2.txt 并检查位置是否相同，但我错过了 inputFile2.txt 中的位置在 inputFile1.txt 中不存在的情况，所以我无法添加那些到 outputFile.txt。

我应该用另一种策略来解决这个任务吗？或者可能是 awk？

Answer 1

我的脚本

#!/usr/bin/env python
import csv
import collections

if __name__ == '__main__':
    db = collections.defaultdict(dict)

    with open('inputFile1.txt') as f:
        reader1 = csv.DictReader(f, delimiter='\t')
        for line in reader1:
            # The for loop will "normalize" the #INDx fields
            for i in range(1, 5):
                k = '#IND{}'.format(i)
                line[k] = line[k].split(':')[0]
            key = (line['#CHR'], line['#POS'])
            del line['#ID'] # Remove this column
            db[key] = line

    with open('inputFile2.txt') as infile:
        reader2 = csv.DictReader(infile, delimiter='\t')
        for line2 in reader2:
            key = (line2['#CHR'], line2['#POS'])
            line1 = db.get(key)
            if line1 is not None:
                line2['#ALT'] = line1['#ALT'] + ',' + line2['#ALT']
            db[key].update(line2)

    with open('outputFile.txt', 'wb') as outfile:
        fieldnames = ['#CHR', '#POS', '#REF', '#ALT',
            '#IND1', '#IND2', '#IND3', '#IND4',
            '#IND5', '#IND6', '#IND7', '#IND8', '#IND9']
        writer = csv.DictWriter(outfile, fieldnames, delimiter='\t', restval='0/0')
        writer.writeheader()
        for k in sorted(db):
            v = db[k]
            writer.writerow(v)

讨论

在我的解决方案中，我使用了csv.DictReader 和csv.DictWriter，因此我可以按名称处理每一列。以下是关于我的变量的一些说明：

• line、line1、line2：每个都是字典，其中的键是各个列的名称，例如 {'#CHR': 'chr1', '#POS':'1000', ... }。您可以按名称引用每一列，例如：line['#ALT']
• db 是一个大字典，其中每个键都是一个元组（#CHR、#POS），每个值是一个单独的行（如上所述，这是另一个字典）。

算法如下：

1. 打开文件 1. 在这一步中，我们去掉#IND1 ... #IND4 列中以冒号开头的文本。我们还删除了 #ID 列，因为输出不需要它。
2. 打开文件 2 并更新列。使用update() 电话更新很容易。但是，必须特别注意#ALT 列。
3. 打开输出文件进行写入。在这一步中，我们通过使用fieldnames 指定列的顺序。 restval 字段告诉 csv 模块使用 0/0 填写任何空白列。然后我们循环遍历db，首先按染色体排序，然后按位置（这就是我将键创建为元组的原因：使排序更容易）。如果您不关心排序，只需删除 db 周围的 sort(...) 函数即可。对于db 中的每个元素，我们将其写入输出文件。

我对我的解决方案唯一关心的是内存使用情况。如果您有非常大的输入文件，您可能会遇到内存错误。告诉我进展如何。