0

我正在从事一个生物信息学项目，该项目涉及非常大的 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）并通过输入 chr 编号和范围并获取正确的行来查询一行（我正在使用它，因为我的队友说这是最快的方法）。

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

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

我希望能够输入 chr num 和 range 并单独获得分数 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 的方法，那将是最好的，但我很难理解它。

谢谢

4

1 回答 1

0

熊猫是严重的矫枉过正。如果您也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查询的输出通过管道传输到 Unixcut实用程序以仅选择第 10 列：

$ tabix input.bed.gz chr1:713835-714424 | cut -f 10
253
于 2020-10-12T17:08:42.713 回答