【生信笔记】python实现DNA反向互补序列的6种方法

1 写在前面的絮絮叨叨

一个练习：用python尽量多种方法来实现DNA反向互补序列。

DNA反向互补序列简单来说就是把序列先倒序读取一次，然后将ATGC对应转换为TACG。

所以把这个过程拆分为反向序列函数和互补序列函数。

python自带的倒序读取方法已经很原始很优秀了，留给我们发挥空间不大，所以重点来看求互补序列可以用哪些骚操作。

1. 首先思路肯定是建立字典再进行字符转换啦【方法1】。

2. translate()是python自带的函数，其实也是建立映射表来实现字符转换，效率很高【方法2】。

3. 最最最原始的方法是用多个if分支来逐个字符判断和替换【方法3】。

4. 也看到有其它博客提到用replace()来实现（Python实现DNA序列互补、反向、反向互补），但是用replace()时要注意对序列的大小写做一点处理，下文再细说【方法4】。

5. 还看到一种用C语言思路实现的方法（生信（二）反向互补序列），是利用了数组的下标对应字符的ASCII码，数组存储互补序列【方法5】。

6. 用replace()构造带{}的替换域，然后利用format对替换域进行替换，妙啊（python编程，获取一段序列的反向互补序列，需要多种方法）【方法5】

如果你还有其它新思路，欢迎留言。
下面马上开始实践吧。

2 反向序列函数

def DNA_reverse(sequence):
    return sequence[::-1]  # 求反向序列

3 互补序列函数

互补序列方法1：用字典dictionary

字典方法很好理解啦，不多说了

# 互补序列方法1：用字典dictionary
def DNA_complement1(sequence):
    # 构建互补字典
    comp_dict = {
        "A":"T",
        "T":"A",
        "G":"C",
        "C":"G",
        "a":"t",
        "t":"a",
        "g":"c",
        "c":"g",
    }
    #求互补序列
    sequence_list = list(sequence)
    sequence_list = [comp_dict[base] for base in sequence_list]
    string = ''.join(sequence_list)
    return string

互补序列方法2：python3 translate()方法

python3自带的str.maketrans(intab, outtab）函数，建立从intab到outtab的映射表。
再用translate(trantab)根据映射表来进行字符转换。
translate()用法参考：https://docs.python.org/zh-cn/3.7/library/stdtypes.html?highlight=translate#str.translate

	# 互补序列方法2：python3 translate()方法
def DNA_complement2(sequence):
    trantab = str.maketrans('ACGTacgtRYMKrymkVBHDvbhd', 'TGCAtgcaYRKMyrkmBVDHbvdh')     # trantab = str.maketrans(intab, outtab)   # 制作翻译表
    string = sequence.translate(trantab)     # str.translate(trantab)  # 转换字符
    return string

互补序列方法3：最原始方法，用多个if分支

# 互补序列方法3：最原始方法，用多个if分支
def DNA_complement3(sequence):
    sequence_list = list(sequence)
    str = []
    for base in sequence_list:
        if base == "A":
            str.append('T')
        elif base == "T":
            str.append('A')
        elif base == "G":
            str.append('C')
        elif base == "C":
            str.append('G')
        elif base == "a":
            str.append('t')
        elif base == "t":
            str.append('a')
        elif base == "g":
            str.append('c')
        elif base == "c":
            str.append('g')
    string = ''.join(str)
    return string

互补序列方法4：对字符串调用replace()

这个方法要注意2个小问题。
1.DNA序列是有大小写的，要先用列表记录原始序列的大小写情况。

seq = "ATGATATAGtatatatgCAAGAGg"  # 原始序列样例

2.replace()是对字符串进行操作，如果先将字符串A转换成T，再将T转换成A的时候就会出问题，此时原始的T和由A转换的T都会被转换掉。解决思路是，先将原始序列都用upper()转换成大写，A->t，T->a，C->g，G->c，就不会出现问题。最后再将原始序列大写的位置转换回大写。

# 互补序列方法4：对字符串调用replace()
def DNA_complement4(sequence):
    sequence_list = list(sequence)
    upper_index=[]

    # 列表upper_index记录哪些位置是小写
    for i in range(len(sequence_list)):
        if sequence_list[i].isupper():
            upper_index.append(i)

    # sequence全部转换为大写后求小写的互补序列
    str = sequence.upper()
    str = str.replace('A', 't')
    str = str.replace('T', 'a')
    str = str.replace('C', 'g')
    str = str.replace('G', 'c')

    # 将原大写序列设回大写
    str_list = list(str)
    for i in upper_index:
        str_list[i] = str_list[i].upper()

    string =''.join(str_list)
    return string

互补序列方法5： ASCII码作为列表下标

# 互补序列方法5： ASCII码作为列表下标 读取对应互补序列
def DNA_complement5(sequence):
    # 建立列表存ASCII码(互补序列方法6)
    comp_tab = []
    for i in range(0, 130):
        comp_tab.append(i)
        trantab = str.maketrans('ACGTacgtRYMKrymkVBHDvbhd', 'TGCAtgcaYRKMyrkmBVDHbvdh')
    for base in trantab:
        comp_tab[base] = chr(trantab[base])

    # comp_tab列表下标是ASCII码：例如'A'的ASCII码是65，comp_tab[65] = 'T'
#    comp_tab = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
#    31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
#     60, 61, 62, 63, 64, 'T', 'V', 'G', 'H', 69, 70, 'C', 'D', 73, 74, 'M', 76, 'K', 78, 79, 80, 81, 'Y', 83, 'A', 85,
#     'B', 87, 88, 'R', 90, 91, 92, 93, 94, 95, 96, 't', 'v', 'g', 'h', 101, 102, 'c', 'd', 105, 106, 'm', 108, 'k', 110,
#     111, 112, 113, 'y', 115, 'a', 117, 'b', 119, 120, 'r', 122, 123, 124, 125, 126, 127, 128, 129]

    string = []
    for j in range(len(sequence)):
        string += comp_tab[ord(sequence[j])]	#求A
    string = ''.join(string)
    return string

互补序列方法6：replace()构造替换域+format()实现替换

# 互补序列方法6：replace()构造替换域+format()实现替换
def DNA_complement6(sequence):
    sequence = sequence.replace('A', '{A}').replace('T', '{T}').replace('C', '{C}').replace('G', '{G}').replace('a', '{a}').replace('t', '{t}').replace('g', '{g}').replace('c', '{c}')
    string = sequence.format(A='T', T='A', C='G', G='C', a='t', t='a', c='g', g='c')
    return string

.replace()构造替换域的过程：

4 测试用例和结果

# seq = "ATGATATAGtatatatgCAAGAGg"    # 临时测试集

# 生成随机序列作为测试集
seq=[]
seq_len = 100000
for i in range(seq_len):
    seq.append(random.choice('ATGCatgc'))
seq = ''.join(seq)
print("原始DNA序列:    seq_len=", seq_len ,"\n", seq, 'len=\n')

# 求反向序列
seq = DNA_reverse(seq)

# 求互补序列
print("DNA反向互补序列:")
for i in [1, 5]:
    begin_time = time()  	# 计时，比较不同方法的性能
    # 依次调用我们写好的6个互补序列函数DNA_complement1/2/3/4/5/6，求结果
    print(eval('DNA_complement'+str(i))(seq))   # eval() 函数用来执行一个字符串表达式，并返回表达式的值。
    end_time = time()
    run_time = end_time - begin_time  
    print('方法', i, '运行时间：', run_time, '\n')

运行结果：

总结

1. 当序列较短时：运行时间差别不大

2. 当序列很长很长时：运行时间 方法2 < 方法1 < 方法6 < 方法3 = 方法4 = 方法5

让这篇文章到你的收藏夹吃灰去吧！
如果文章对你有帮助，请留言或者点个赞，这将是对我极大的支持，谢谢。
欢迎评论欢迎评论欢迎评论~