python - 使用 python 计算 gcskews 的问题

Question

每个人。

我在 python 中计算 gcskews 时遇到了一些问题。我的两个主要输入是 fasta 文件和 bed 文件。床文件有 gn(0)、gene_type(1)、基因名称(2)、染色体(3)、strand(4)、num(5)、start(6) 列。（这些数字是 python 中的索引号。 ) 然后我尝试使用一些函数，这些函数可以从每个基因的起始位点计算有义和反义链的 gcskew。窗口是 100bp，这些是功能。

import re
import sys
import os

# opening bed file
content= []
with open("gene_info.full.tsv") as new : 
 for line in new :
     content.append(line.strip().split())
content = content[1:]       

def fasta2dict(fil):
     dic = {}
     scaf = ''
     seq = []
     for line in open(fil):
         if line.startswith(">") and scaf == '':
             scaf = line.split(' ')[0].lstrip(">").replace("\n", "")
         elif line.startswith(">") and scaf != '':
             dic[scaf] = ''.join(seq)
             scaf = line.split(' ')[0].lstrip(">").replace("\n", "")
             seq = []
         else:
             seq.append(line.rstrip())
     dic[scaf] = ''.join(seq)
     return dic

dic_file = fasta2dict("full.fa")

# functions for gc skew
def GC_skew_up(strand, loc, seq, window = 100) : # need -1 for index
  values_up = []
  loc = loc - 1
  if strand == "+" :
    sp_up = seq[loc - window : loc]
    g_up = sp_up.count('G') + sp_up.count('g') 
    c_up = sp_up.count('C') + sp_up.count('c') 
    try :
        skew_up = (g_up - c_up) / float(g_up + c_up)
    except ZeroDivisionError:
        skew_up = 0.0
    values_up.append(skew_up)
  elif strand == "-" :
    sp_up = seq[loc : loc + window]
    g_up = sp_up.count('G') + sp_up.count('g') 
    c_up = sp_up.count('C') + sp_up.count('c')
    try :
        skew_up = (c_up - g_up) / float(g_up + c_up)
    except ZeroDivisionError:
        skew_up = 0.0
    values_up.append(skew_up)

  return values_up  

def GC_skew_dw(strand, loc, seq, window = 100) :
  values_dw = []
  loc = loc - 1
  if strand == "+" :
    sp_dw = seq[loc : loc + window]
    g_dw = sp_dw.count('G') + sp_dw.count('g') 
    c_dw = sp_dw.count('C') + sp_dw.count('c') 
    try :
        skew_dw = (g_dw - c_dw) / float(g_dw + c_dw)
    except ZeroDivisionError:
        skew_dw = 0.0
    values_dw.append(skew_dw)
  elif strand == "-" :
    sp_dw = seq[loc - window : loc]
    g_dw = sp_dw.count('G') + sp_dw.count('g') 
    c_dw = sp_dw.count('C') + sp_dw.count('c')
    try :
        skew_dw = (c_dw - g_dw) / float(g_dw + c_dw)
    except ZeroDivisionError:
        skew_dw = 0.0
    values_dw.append(skew_dw)   

  return values_dw  

正如我所说，我想从基因的起始位点计算 100bp 链的 gcskews。
因此，我编写了从 bed 文件中获取染色体名称并从 Fasta 文件中获取序列数据的代码。
然后根据基因名称和链信息，我预计代码会找到正确的起始位点并计算 100bp 窗口的 gcskew。
但是，当我运行此代码时， - strand 的 gcskew 是错误的，但 + strand 是正确的。（我得到了正确的 gcskew 数据并使用了它。）
gcskew 与正确的数据不同，但我不知道是什么问题。
谁能告诉我这段代码有什么问题？
提前致谢！

window = 100
gname = []
up = []
dw = []

for match in content :
seq_chr = dic_file[str(match[3])]
if match[4] == "+" :
    strand = match[4]
    new = int(match[6])
    sen_up = GC_skew_up(strand, new, seq_chr, window = 100)
    sen_dw = GC_skew_dw(strand, new, seq_chr, window = 100)
    gname.append(match[2])
    up.append(str(sen_up[0]))
    dw.append(str(sen_dw[0]))

if match[4] == "-" :
    strand = match[4]
    new = int(match[6])
    an_up = GC_skew_up(strand, new, seq_chr, window = 100)
    an_dw = GC_skew_dw(strand, new, seq_chr, window = 100)
    gname.append(match[2])
    up.append(str(an_up[0]))
    dw.append(str(an_dw[0]))

tot = zip(gname, up, dw)

标签： pythonfunctionbioinformatics