RNA-seq数据与特定基因的关联

Question

我有一个基因列表(作为一个bed文件)和一个全基因组RNA-seq数据集(也存储为一个bed文件)。我目前正在尝试开发一个python脚本，它允许我提取转录起始点上游500bp到下游2000bp的读取计数，即基因的开始，并将这些值存储在一个数组中以备将来使用。

目前，我的脚本如下所示：

feature_genes=np.zeros((6576, 2501))

for lines in feature:
   for i in range(0,6575):
        if line[5]==lines[5] and line[5]=='+' and line[0]==lines[0] and int(lines[1])>=int(line[1])- 500 and int(lines[1])<=int(line[1])+2000:
            feature_genes[i][int(lines[1])-int(line[1])+500]=lines[4] 
        elif line[5]=='-' and line[0]==lines[0] and int(lines[2])+500>=int(line[2]) and int(lines[2])-2000<=int(line[2]) and lines[5]=='-':
            feature_genes[i][-1*(int(lines[2])-int(line[2])-500)]=lines[4]

其中feature代表我在文件中的读数，基因我的基因列表，每一行都包含特定核苷酸(这是链特定信息，不包括任何未观察到读数的碱基对)或基因位置的读数。

注意：.bed文件的格式如下：

Position 

0 chromosome
1 transcription start site
2 transcription termination site
3 feature name
4 read count
5 strand 

有没有人能想出一个有效的方法呢？我的代码永远运行不完(python菜鸟)。

Answer 1

简单的答案是不使用python，而是使用床上工具。有几种方法可以做到这一点，下面是其中一种：

1)将TSS上游扩展x个核苷酸，并将其下游扩展x个核苷酸，这样数学就已经得到了处理。

2)使用带有abam选项的intersectBed输出覆盖感兴趣区域的RNA-Seq读数(或者，如果只想要覆盖深度，可以使用coverageBed)