我尝试遍历一个文件中的每个特性(每行1个),并根据第二个文件中该行的一列查找所有匹配的特性。我有这个解决方案,在小文件上做我想要的,但在大文件上非常慢(我的文件有20000000行)。Here's a sample of the two input files.
我的(慢)代码:
FEATUREFILE = 'S2_STARRseq_rep1_vsControl_peaks.bed'
CONSERVATIONFILEDIR = './conservation/'
with open(str(FEATUREFILE),'r') as peakFile, open('featureConservation.td',"w+") as outfile:
for line in peakFile.readlines():
chrom = line.split('\t')[0]
startPos = int(line.split('\t')[1])
endPos = int(line.split('\t')[2])
peakName = line.split('\t')[3]
enrichVal = float(line.split('\t')[4])
#Reject negative peak starts, if they exist (sometimes this can happen w/ MACS)
if startPos > 0:
with open(str(CONSERVATIONFILEDIR) + str(chrom)+'.bed','r') as conservationFile:
cumulConserv = 0.
n = 0
for conservLine in conservationFile.readlines():
position = int(conservLine.split('\t')[1])
conservScore = float(conservLine.split('\t')[3])
if position >= startPos and position <= endPos:
cumulConserv += conservScore
n+=1
featureConservation = cumulConserv/(n)
outfile.write(str(chrom) + '\t' + str(startPos) + '\t' + str(endPos) + '\t' + str(peakName) + '\t' + str(enrichVal) + '\t' + str(featureConservation) + '\n')
首先,每次从
peakFile
读取一行时,您都会遍历所有的conservationFile
,因此在if语句中的n+=1
之后插入一个break
,这应该会有所帮助。假设只有一个匹配。你知道吗另一种选择是尝试使用mmap,这可能有助于缓冲
Bedtools就是为此而设计的,特别是
intersect
函数:http://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
对于我来说,最好的解决方案似乎是为熊猫重写上面的代码。以下是在一些非常大的文件上对我有效的方法:
相关问题 更多 >
编程相关推荐