用Python量化对引用的读取

Question

输入文件（test.sam）：

SN398:470:C8RD3ACXX:7:1111:19077:53994  16  chrI    65374   255 51M *   0   0   TGAGAAATTCTTGAACATTCGTCTGTATTGATAAATAAAACTAGTATACAG IJJJJJJJJJJJJJIJJJIJJJJJJHJJJJJJJJJJJJHHHHHFFFFDB@B AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51 YT:Z:UU NH:i:1

genes.bed文件是参考文件：

chrI    130798  131983  YAL012W 0   +   130798  131983  0   1   1185,   0,
chrI    334 649 YAL069W 0   +   334 649 0   1   315,    0,
chrI    537 792 YAL068W-A   0   +   537 792 0   1   255,    0,
chrI    1806    2169    YAL068C 0   -   1806    2169    0   1   363,    0,
chrI    2479    2707    YAL067W-A   0   +   2479    2707    0   1   228,    0,
chrI    7234    9016    YAL067C 0   -   7234    9016    0   1   1782,   0,
chrI    10090   10399   YAL066W 0   +   10090   10399   0   1   309,    0,
chrI    11564   11951   YAL065C 0   -   11564   11951   0   1   387,    0,
chrI    12045   12426   YAL064W-B   0   +   12045   12426   0   1   381,    0,

脚本如下-它查看两个文件之间的“chr”是否匹配，如果test.sam的第四列（称为genomic_location）位于genes.bed文件的第二列和第三列内，则它将打印genes.bed的第四列并将其计为“1”

#!/usr/bin/env python
import sys

samfile=open('test.sam')  #sorted sam file
bedfile=open('genes.bed') #reference genome
sys.stdout=open('merged.txt', 'w')

lookup = {}
for line in bedfile:
   fields = line.strip().split()
   chrm   = fields[0]
   st     = int(fields[1])
   end    = int(fields[2])
   name   = fields[3]
   if chrm not in lookup:
       lookup[chrm] = {}
   for i in range(st,end):
       if i not in lookup[chrm]:
           lookup[chrm][i] = [name]
       else:
           lookup[chrm][i].append(name)

gene_counts = {}
for line in samfile:
   reads = line.split()
   qname = reads[0]
   flag  = reads[1] # be 0 or 16
   rname=reads[2]
   genomic_location = int(reads[3])
   mapq  = int(reads[4])
   if rname in lookup:
       if genomic_location in lookup[rname]:
           for gene in lookup[rname][genomic_location]:
               if gene not in gene_counts:
                   gene_counts[gene]  = 0
           else:
               gene_counts[gene] += 1

print gene_counts

我需要改变它，当flag（输入文件test.sam中的第二列）为16时，然后从输入文件（test.sam）中的第四列减去51，然后对其进行处理，以查看新生成的整数是否在genes.bed文件的st和end范围内

你认为最好的方法是什么？我需要在脚本中实现这一点，而不是创建一个新的输入文件（test.sam），如果第二列是16，它将首先更改第四列

我想做这个Python。谢谢你的帮助，如果有什么不清楚的地方，请告诉我

Answer 1

也许这里我忽略了一些隐藏的复杂性，但是“when flag（输入文件test.sam中的第二列）是16，然后从inputfile中的第四列减去51”最明显的python实现是：

if flag == 16:
     genomic_location = int(reads[3]) - 51