如何最有效地从NarrowPeakBED6 + 4格式的文件中检索数据？

问题描述

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

（列为“ chrom，chr​​omStart，chromEnd，名称，得分，链，signalValue，pValue，qValue，峰值”）

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）并通过输入chr编号和范围并获取正确的行来查询行（我正在使用它，因为我的队友说这是最快的方法）。 / p>

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

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

我希望能够输入字符数和范围并单独获得253分。

我通过将结果转换为一个数据框（带有大熊猫）来实现的。

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的方法，那将是最好的方法，但是我很难理解它。

谢谢

解决方法

熊猫是严重的过度杀伤力。如果您也使用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实用程序以仅选择第十列：

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