如何最有效地从 NarrowPeak (BED6+4) 格式文件中检索数据？

Question

我正在从事一个生物信息学项目，该项目涉及非常大的 NarrowPeak 格式文件，如下所示：

（列为 'chrom ,chromStart, chromEnd, name, score ,strand, signalValue ,pValue, qValue ,peak'）


chr1    713835  714424  chr1.1  1000    .   0.1621  10.6    -1  253
chr1    752775  753050  chr1.2  567     .   0.0365  2.09    -1  124
chr1    753373  753448  chr1.3  511     .   0.0243  1.31    -1  27
chr1    762057  763210  chr1.4  1000    .   0.1743  11.5    -1  846
chr1    793358  793612  chr1.5  574     .   0.0379  2.18    -1  121
chr1    805101  805538  chr1.6  783     .   0.0834  5.23    -1  200
chr1    839626  840461  chr1.7  1000    .   0.2079  13.8    -1  510
chr1    842125  842534  chr1.8  641     .   0.0526  3.15    -1  199
chr2    846334  846381  chr1.9  510     .   0.0241   1.3    -1  17
chr2    846545  846937  chr1.10 562     .   0.0355  2.03    -1  198
chr2    848219  848588  chr1.11 605     .   0.0448  2.64    -1  179
chr2    851784  852887  chr1.12 734     .   0.0728  4.51    -1  323

我有一种方法可以为文件编制索引（在 samtools 中使用 tabix）并通过输入字符编号和范围并获取正确的行来查询行（我使用它是因为我的队友说这是最快的方法） .

例如，如果我的查询是 chr1:713835-714424，我将得到结果：

chr1    713835  714424  chr1.1  1000    .   0.1621  10.6    -1  253

我希望能够输入chr num和range，单独得到253分

我通过将结果转换为数据框（pandas）来做到这一点。

def turn_query_results_into_df(query_results):
    replaced=(query_results.replace('\t',','))
    stripped = (replaced.strip())
    lines = (replaced.split("\n") )
    data=[]
    for line in lines:
        data.append((line.split(',')))
    mydf1=pd.DataFrame(data,columns=('chrom chromStart chromEnd name score strand signalValue pValue qValue peak'.split()))
    mydf1=mydf1.drop(mydf1.tail(1).index)
    mydf1=mydf1.astype({'chromStart':'int32','chromEnd':'int32','score':'int32','signalValue':'float64','pValue':'float64','peak':'float64'})
    return mydf1

print(peak_val = df1['peak'].values)

有没有更好的方法？

如果有samtools的使用方法，那最好不过，但我很难理解。

谢谢

Answer 1

Pandas 是严重的矫枉过正。如果您也使用 tabix 进行查询，则使用如下命令序列：

$ tabix input.bed.gz # index the input
$ tabix input.bed.gz chr1:713835-714424 # query the input
chr1    713835  714424  chr1.1  1000    .   0.1621  10.6    -1  253

您可以将 tabix 查询的输出通过管道传输到 Unix cut 实用程序到 select 只有第 10 列：

$ tabix input.bed.gz chr1:713835-714424 | cut -f 10
253

如何最有效地从 NarrowPeak (BED6+4) 格式文件中检索数据？

How to most efficiently retrieve data from NarrowPeak (BED6+4) format files?

bioinformatics

pandas

samtools