-2

我有一个看起来像这样的文件：

Chr-坐标覆盖

chr1    236968289   2
chr1    236968318   2
chr1    236968320   2
chr1    236968374   2
chr1    237005709   2
chr14   22086843    2
chr14   22086846    2
chr14   22086849    2
chr14   22086851    4
chr2    5078129 2
chr2    5341758 2
chr2    5342443 2

我想操纵它来获得：

chr-开始-结束-平均覆盖距离

chr1  236968289  236968374 2    85    
chr14 22086843   22086851  2.5  8
chr2  5078129    5078129   2    0
chr2  5341758    5342443   2    685

我想要：如果 chr 与前一个 chr不同，或者坐标之间的差异大于 1000：它会打印输出，如图所示。与 chr、起始坐标、结束坐标、平均覆盖范围以及开始和结束之间的距离。

为此，我编写了以下代码：

cov=open("coverage.txt")
oldchr="chr55"      #dummy starting data
oldcoordinate=1
sumcoverage=0
startcoordinate=0

try:
   while True:
    line=next(cov).split("\t",2)
    newchr=line[0]
    newcoordinate=int(line[1])       #read informations from file
    newcoverage=int(line[2].strip())

    if oldchr != newchr or newcoordinate - oldcoordinate > 1000:
        distance=oldcoordinate-startcoordinate
        averagecoverage=sumcoverage/distance
        merge=oldchr+'\t'+str(startcoordinate)+'\t'+str(oldcoordinate)+'\t'+str(averagecoverage)+'\t'+str(distance) 
        print merge
        startcoordinate=newcoordinate
        sumcoverage=0

    oldchr=newchr
    oldcoordinate=newcoordinate     #replace old with new chr and coordinates
    sumcoverage=sumcoverage+newcoverage


except(StopIteration):
    print ""

我无法理解为什么它不能正常工作。我得到的错误是获得“平均覆盖率”的除法试图除以 0，所以在许多情况下，“距离”（ 距离 = oldcoordinate-startcoordinate）等于 0。这不应该发生在输入中文件绝不是两行具有相同坐标的情况。我无法看到错误在哪里。我希望有人可以帮助我，提前谢谢你。

4

1 回答 1

0

您可以使用 Python 的groupby函数根据chr列对条目进行分组。该csv库还可以更轻松地处理文件：

from itertools import groupby
import csv

def display_block(block):
    average_coverage = sum(x[2] for x in block) / float(len(block))
    print block[0][0], "\t", block[0][1], "\t", block[-1][1], "\t", average_coverage, "\t", block[-1][1] - block[0][1]


with open('coverage.txt', 'rb') as f_coverage:
    for chr, entries in groupby(csv.reader(f_coverage, delimiter='\t'), lambda x: x[0]):
        entries = [(e[0], int(e[1]) , int(e[2])) for e in entries]
        block = []
        ientries = iter(entries)
        block.append(next(ientries))

        for chr, coord, coverage in ientries:
            if coord - block[-1][1] > 1000:
                display_block(block)
                block = []
            block.append([chr, coord, coverage])

        if len(block):
            display_block(block)

这将显示以下输出：

chr1    236968289   236968374   2.0     85
chr1    237005709   237005709   2.0     0
chr14   22086843    22086851    2.5     8
chr2    5078129     5078129     2.0     0
chr2    5341758     5342443     2.0     685

循环的每次迭代for都会为您提供当前chr和所有匹配的行chr。该脚本遍历每一行并将 and 转换coord为coverage整数。然后它将列表变成一个迭代器。存储第一个匹配行block，然后迭代任何剩余的行。每次距离大于1000块时都会显示并重新启动。然后还会显示最后的任何剩余条目。

使用 Python 2.7.6 测试

于 2016-02-11T14:33:52.453 回答