使用Python计算DNA链的反向互补序列

Question

我有一段DNA序列，想用Python来获取它的反向互补序列。这段序列在一个CSV文件的某一列里，我想把反向互补序列写到同一个文件的另一列。问题是，有些单元格里包含的不是A、T、G和C这几种碱基。我用下面这段代码成功得到了反向互补序列：

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} 
    bases = list(seq) 
    bases = [complement[base] for base in bases] 
    return ''.join(bases)
    def reverse_complement(s):
        return complement(s[::-1])

    print "Reverse Complement:"
    print(reverse_complement("TCGGGCCC"))

但是，当我尝试用下面的代码查找那些不在互补字典里的项目时，我只得到了最后一个碱基的互补，而没有遍历所有的碱基。我想知道该怎么解决这个问题。

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} 
    bases = list(seq) 
    for element in bases:
        if element not in complement:
            print element  
        letters = [complement[base] for base in element] 
        return ''.join(letters)
def reverse_complement(seq):
    return complement(seq[::-1])

print "Reverse Complement:"
print(reverse_complement("TCGGGCCCCX"))

Answer 1

反向互补序列最快的一行代码是下面这个：

def rev_compl(st):
    nn = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(nn[n] for n in reversed(st))

Answer 2

import string
old_chars = "ACGT"
replace_chars = "TGCA"
tab = string.maketrans(old_chars,replace_chars)
print "AAAACCCGGT".translate(tab)[::-1]

这段代码会给你一个反向的互补序列，结果是 ACCGGGTTTT。

Answer 3

一般来说，生成器表达式比原来的代码更简单，而且不会创建额外的列表对象。如果需要插入多个字符，可以参考其他的答案。

complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
seq = "TCGGGCCC"
reverse_complement = "".join(complement.get(base, base) for base in reversed(seq))

Answer 4

其他的回答都很好，但如果你打算处理真实的DNA序列，我建议你使用Biopython。因为如果你遇到像“-”、“*”这样的字符，或者一些不确定的情况怎么办？如果你想对你的序列进行进一步的操作，你想为每种文件格式都写一个解析器吗？

你问的代码其实很简单：

from Bio.Seq import Seq

seq = Seq("TCGGGCCC")

print seq.reverse_complement()
# GGGCCCGA

如果你想进行其他转换：

print seq.complement()
print seq.transcribe()
print seq.translate()

输出结果

AGCCCGGG
UCGGGCCC
SG

而且如果你碰到奇怪的字符，也不需要一直往你的程序里加代码。Biopython会帮你处理这些问题：

seq = Seq("TCGGGCCCX")
print seq.reverse_complement()
# XGGGCCCGA

Answer 5

字典中的 get 方法可以让你在查找一个键时，如果这个键不在字典里，可以指定一个默认值。作为一个预处理步骤，我建议你把所有不是 'ATGC' 的碱基用单个字母（或者标点符号、数字，反正是不会出现在你的序列中的东西）来表示，然后把序列反转，再把这些单个字母替换回它们原来的样子。另一种方法是先反转序列，然后再把像 sni 这样的内容替换成 ins。

alt_map = {'ins':'0'}
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'} 

def reverse_complement(seq):    
    for k,v in alt_map.iteritems():
        seq = seq.replace(k,v)
    bases = list(seq) 
    bases = reversed([complement.get(base,base) for base in bases])
    bases = ''.join(bases)
    for k,v in alt_map.iteritems():
        bases = bases.replace(v,k)
    return bases

>>> seq = "TCGGinsGCCC"
>>> print "Reverse Complement:"
>>> print(reverse_complement(seq))
GGGCinsCCGA