2

我目前在一个文件中有一个基因列表。每条线都有一条带有它的信息的染色体。这样的条目显示为：

NM_198212 chr7 + 115926679 115935830 115927071 11593344 2 115926679 ,'115933260', 115927221 ,'115935830',

染色体的序列从碱基115926679开始，一直到（但不包括）碱基115935830

如果我们想要拼接序列，我们使用外显子。第一个从 115926679延伸到155927221，第二个从 '115933260' 到 '115935830'

但是，我在互补序列上遇到了一个问题，例如：

NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1 245941755，'245942680'

由于第 3 列是“-”，因此这些坐标是指反义链（链的补码）。第一个碱基（粗体）与有义链上的最后一个碱基（斜体）匹配。由于文件只有有义支架，我需要尝试将反义链上的坐标转换为有义链，挑选出正确的序列，然后对其进行反向补码。

也就是说，我只编程了大约半年，并且不知道如何开始这样做。

我写了一个正则表达式：

'(NM_\d+)\s+(chr\d+)([(\+)|(-)])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+),(\d+),s+(\d+),(\d+),'

但现在我不确定如何启动这个功能......如果有人可以帮助我开始这方面的工作，也许让我知道如何做到这一点，我将非常感激。

好的：假设这是 25 号染色体：

啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊！

（每个字符有 10 个）。

现在：如果我正在寻找未剪接的基因：chr25 + 10 20

然后基因从位置 10（从 0 开始）开始，向上但不包括位置 20。所以它：

中交中交中交

这很简单。它非常适合 python 字符串切片。

如果我给你，它会更令人困惑：

chr25 - 10 20

你所拥有的是积极的一面。但是这个基因在负（互补）链上。记住染色体看起来像双链：

啊啊
啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊！

我们正在寻找底部链上的基因。这意味着我们从右边开始从 0 开始计数。上链从左开始编号，下链从右开始。所以我在这里想要的是 AAAAAAAAAA。

问题是我只给你最上面的部分。我不会给你底线。（你可以从顶层生成自己——但考虑到它有多大，我建议不要这样做。）

所以你需要转换坐标。在底链上，碱基 0（最右边的 C）与顶链上的碱基 39 相对。以 1 为底数与 38 为底数。以 2 为底数以 37 为底数。（重要一点：注意每次将这两个数字相加时会发生什么。）因此，以 10 为底数以 29 为底数，以 19 为底数以 20 为底数。

所以：如果我想在底部链上找到基数 10-20，我可以查看顶部的基数 20-29（然后反向补码）。

我需要弄清楚如何将底部链上的坐标转换为底部链上的等效坐标。是的：这很混乱

我试过weronika的原始答案：

fields = line.split(' \t')
geneID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
    start,end = -(start + 1), -(end + 1) # this part was changed from original answer.

这是在正确的轨道上，但这还不够。这将把 10 和 20 变成 20 和 10。

而且我知道我可以通过这样做来反向补充字符串：

r = s[::-1]
bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
l = list(r)
o = [bc[base] for base in l]
     return ''.join(o)

已编辑！这看起来对吗？！

fp2 = open('chr22.fa', 'r')
fp = open('chr22.fa', 'r')
for line in fp2:
    newstring = ''
    z = line.strip()
    newstring += z
for line in fp:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]
    start = int(fields[3])
    end = int(fields[4])
    bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'a': 't', 't': 'a', 'c':'g', 'g':'c', 'N':'N', 'n':'n'}
    l = list(newstring)        
    if strand == '+':
        geneseq = ''.join([bc[base] for base in l[start:end]]) 
    if strand == '-':
        newstart, newend = -(start + 1), -(end + 1)
        genseq = ''.join([bc[base] for base in l[newstart:newend]])
4

4 回答 4

0

我不明白域问题，但看起来你试图在一个正则表达式中塞进太多东西。尝试将其分解为更简单的子问题，如下所示（在伪代码中）：

if third column is a '+'
    parseRegularSequence()
else
    parseComplementarySequence()
于 2012-04-21T21:51:22.620 回答
0

正如我在对该问题的评论中所指出的，这似乎是一种非常奇怪的文件格式，因此我最初感到困惑。

注意：如果这是标准生物学文件格式之一，您最好使用Biopython或类似的东西来解析它。

如果您想进行自己的解析，正则表达式似乎仍然是错误的方法 - 使用简单的空格/制表符分隔文件难以阅读且不必要。我假设你已经解析了你的染色体 fasta 文件，并将所有染色体的序列作为chrom_sequencesname:seq 字典，并且你有一个reverse_complement函数（这两件事都可以很容易地手动实现，但用 biopython 可能更好） .

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start, end = [int(x) for x in fields[3:5])
this_chromosome_seq = chrom_sequences[chr]
# if strand is +, just take the sequence based on the start-end position
if strand == '+':
  # be careful about off-by-one errors here - some formats give you a 1-based position, 
  #  other formats make it 0-based, and they can also be end-inclusive or end-exclusive
  gene_sequence = this_chromosome_seq[start:end]
# if your coordinates really are given as antisense strand coordinates when strand is -,
#  you just need to subtract them from the chromosome length to get sense-strand coordinates,
#  (switching start and end so they're still in smaller-to-larger order),
#  and then reverse-complement the resulting sequence. 
if strand == '-':
  chrom_length = len(this_chromosome_seq)
  # again, be careful about off-by-one errors here!
  new_start,new_end = chrom_length-end, chrom_length-start
  gene_sequence = reverse_complement(this_chromosome_seq[new_start:new_end])

原始答案，实际上并没有按照要求做：

如果您只想获取开始/结束位置，请执行以下操作：

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
  start,end = end,start

然后，您必须解析您的 fasta 文件以实际获取这些起点坐标的序列，并在strand=='-'. 同样，我认为 Biopython 可以为您完成大部分工作。

另一个注意事项 - Biostar是一个很好的 StackExchange 类型的站点，专门用于生物信息学，您可能会在那里得到更好的答案。

于 2012-04-21T22:42:04.927 回答
0

如果用负数对字符串进行切片，Python 会从末尾倒数。

complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
s = "AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG"
"".join(complement[c] for c in s[-20:-10])

编辑：

你的编辑看起来对我来说是正确的，是的。我非常不擅长检查栅栏错误，但无论如何你比我更能看到这些错误！

我已经将您的代码重写为更加 Pythonic，但没有改变任何实质性内容。

bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N':'N'}

f = open('chr22.fa')
l = ''.join(line.strip() for line in f)
f.seek(0)

for line in f:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]

    start = int(fields[3])
    end = int(fields[4])

    if strand == '-':
        start, end = -(start + 1), -(end + 1)

    geneseq = ''.join(bc[base.upper()] for base in l[start:end])
于 2012-04-22T16:16:20.543 回答
0

我想（尤其是因为您的文件很大）直接从文件缓冲区读取和写入会容易得多。

假设您已经解析了头文件。您解析的行如下所示：

line = "NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1"

然后，您确定开始/结束位置是什么（反义坐标）：

name, chromosome, direction, start, end = line[:5]

然后，执行以下操作：

#Open up the file `chr1.txt`.
f = open("chr1.txt", "r")

#Determine the read length.
read_len = end - start

#Seek to the appropriate position (notice the second argument, 2 -- this means
#seek from the end of the file)
f.seek(-end, 2)

#Read the data
sense_seq = f.read(read_len)

之后，只需将序列转换为反义即可。

简单的 ：）

于 2012-04-23T01:15:21.523 回答