背景

最近参加了个生信的面试，记录一下有意思的面试题。

题目描述

要求从提供的*.fasta文件出发：

1. 获得序列的反向互补序列，并统计信息：序列条数，碱基总数，N50，N90，GC 含量。
   
2. 提取每条序列上 32bp-332bp、780bp-992bp 的序列。
   
3. 统计单碱基重复 4 次及以上的序列在每条序列上出现的次数。如 AAAAA 或者 TTTT 等

   • 如 “AAATTTTTTTCCCCAAAAAAA”，结果如下：

     • the 4th character T repeat 7 times
       
     • the 11th character C repeat 4 times
       
     • the 15th character A repeat 7 times
       
4. 将序列拼接起来，中间添加 300 个 N 从而形成一条序列，保存的结果文件中每行 60 个碱基。
   

Code

将fasta文件读入为字典

使用方法见代码内注释

• 要点：通过fasta文件中的">"识别序列的描述，存为字典的“key”,把序列信息存为字典的“value”
  

# 功能：读取解压后的*.fasta文件，将内容按键值对形式存为字典# 输入： ## fasta_name：解压后的文件名# 输出：## fa_dict：包含描述和序列的字典def fasta2dict(fasta_name):        with open(fasta_name) as fa:        fa_dict = {}        for line in fa:            # 去除末尾换行符            line = line.replace('\n','')            if line.startswith('>'):                # 去除 > 号                seq_name = line[1:]                fa_dict[seq_name] = ''            else:                # 去除末尾换行符并连接多行序列                fa_dict[seq_name] += line.replace('\n','')    return fa_dict

• 根据字典长度获取序列条数
  

my_fasta_dict = fasta2dict("exercise.fa")print("序列条数为 : %d" %len(my_fasta_dict))

获取反向互补序列

• 要点：先获取互补序列，再获取反向序列（通过切片实现）
  

# 功能：将给定的碱基序列转化为反向互补序列def DNA_rever_complement(sequence):    trantab = str.maketrans('ACGTacgt', 'TGCAtgca')     # trantab = str.maketrans(intab, outtab)   # 制作翻译表    string = sequence.translate(trantab)     # str.translate(trantab)  # 转换字符    string_tmp = list(string)[::-1]    string_fin = ''.join(string_tmp)

    return string_fin

获取多个信息

• 要点：

  • 获取fasta文件反向互补序列
    
  • 计算fasta文件总碱基个数
    
  • N50/N90、GC含量
    
  • 返回指定位置指定长度的碱基
    
  • 将所有碱基以300连续“N”作为间隔得到一整条序列
    

# 功能： 获取fasta文件反向互补序列、计算fasta文件总碱基个数、N50/N90、GC含量、返回指定位置指定长度的碱基、将所有碱基以300连续“N”作为间隔得到一整条序列# 输入：## fasta_dict：值为序列信息的字典## N_percent: N50时为50, N90时为90, 也可以为0~100内的其他整数## *cut_local：可变长度参数，指定需要提取的碱基位置信息，数据为偶数且为整数。注：此参数可以缺省，表示不进行”返回指定位置指定长度的碱基“的计算# 输出：总碱基长度，N50或N90以及任意Nn，GC含量，截取指定长度碱基存入字典，将所有碱基以300连续“N”作为间隔得到一整条序列## rever_comple_seq：获取fasta文件反向互补序列，输出为字典## all_dna_len: fasta文件中总碱基数## N_len：N50或N90长度, 输入参数（即，N_percent）需要对应 ## GC_percent：GC含量## fasta_dict_cut：获得指定位置，指定长度的碱基，输出为字典## all_dna_fin：将所有碱基以300连续“N”作为间隔得到一整条序列，输出为列表def return_Nn(fasta_dict, N_percent:int, *cut_local:int):    import sys

    rever_comple_seq = {}    fasta_dict_counts = {}    fasta_dict_cut = {}    all_dna = []    all_dna_len = 0    len_tmp = 0    nGC = 0

    # 将原来字典中碱基序列换成序列长度    ## key为序列描述信息，value为碱基信息    for key,value in fasta_dict.items():        # GC碱基数        nGC_tmp = value.count("g") + value.count("G") + value.count("c") + value.count("C")        nGC += nGC_tmp         # 在每条序列后面接300个连续“N”        all_dna_tmp = value + 'N'*300        all_dna += all_dna_tmp

        # 调用上述自定义的“DNA_rever_complement”函数， 获取反向互补序列，并存为字典        rever_comple_seq[key] = DNA_rever_complement(value)        fasta_dict_counts[key] = len(value)

        # 获取指定长度，指定位置碱基序列        ## 从第三个参数开始为指定碱基的位置参数（即，cut_local参数），可选参数        ## 每两个为一对，该参数个数必须为偶数        cut_dna = []        # 技巧：以步长为2取cut_local参数信息        for i in range(len(cut_local))[::2]:            if i + 1 < len(cut_local):                cut_dna.append(value[cut_local[i]:cut_local[i + 1]])            else:                print("Error:the lengths of cut_local parameters must be even!!!")                 sys.exit()        fasta_dict_cut[key] = cut_dna

        # *.fasta文件总碱基数        all_dna_len += len(value)    # 所有碱基整合为一条序列    all_dna_fin = ''.join(all_dna)    # GC含量    GC_percent = nGC/all_dna_len    #print("GC_percent:%f" %GC_percent)    # N50或N90    cutoff = int(all_dna_len*(N_percent/100))    # 将得到的序列长度字典按序列长度降序排序，用于计算N50,N90    fasta_dict_counts_order = dict(sorted(fasta_dict_counts.items(), key=lambda k: k[1], reverse = True)) # 按碱基长度排序得到新字典    # 注：lambda为匿名函数    for key,value in fasta_dict_counts_order.items():        len_tmp += value        if len_tmp >= cutoff:            #print("N%d为：%d" %(N_percent, value))            N_len = value            break

    return rever_comple_seq, all_dna_len, N_len, GC_percent, fasta_dict_cut, all_dna_fin

sum_1 = return_Nn(my_fasta_dict, 90, 0, 30, 10, 40)

将列表内长字符串以指定长度输出至文件

# 功能：将列表内长字符串以指定长度输出至文件# 参数：## input_file: 待输入的字符串，可以是列表或'str'## row_length: 每行指定输出的长度## output_file: 输出文件名称

def outPutStrLength(input_file, row_length, output_file):    import math    startpoint = 0    input_file = input_file.replace('\n', '') # 替换掉字符串中换行符    seg = math.ceil(len(input_file)/row_length) # 每行输出row_length 个向上取整有seg行，    with open(output_file, 'w') as f:        for i in range(seg):            startpoint = row_length * i  #每行的索引点

            f.write(input_file[startpoint : startpoint + row_length] + "\n") # 索引字符串            # print(tempStr[startpoint : startpoint + row_length]) 

outPutStrLength(sum_1[-1], 60, "fasta_test.fa")

统计单碱基重复 4 次及以上的序列在每条序列上出现的次数

• 要点：给每条序列添加一个符号位
  
• 如：
  
• AAATCGGGGTCCCC
  
• 01200012300123
  

# 功能：统计单碱基重复 4 次及以上的序列在每条序列上出现的次数，并以“The 4th character A repeat 6 times”形式输出# 输入：## fasta_dict：字典，key为序列说明，value为序列信息## repeat_times_least: 整型，输出的碱基最少重复次数# 输出：## 格式化输出

def repeat_compute(fasta_dict, repeat_times_least:int):    for key,value in fasta_dict.items():        print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")        print(key)        test_str = list(value)

        flag_list = [0]        flag = 0        Zero_flag_index = 0        Zero_flag_index_list = []        for i in range(len(test_str) - 1):            Zero_flag_index += 1            if test_str[i] == test_str[i + 1]:                flag += 1                flag_list.append(flag)            else:                flag = 0                flag_list.append(flag)                Zero_flag_index_list.append(Zero_flag_index)

        need_index_tmp = [Zero_flag_index_list[i] for i in range(len(Zero_flag_index_list)) if flag_list[Zero_flag_index_list[i] - 1] >= repeat_times_least -1]        if len(flag_list) - Zero_flag_index_list[-1] >= repeat_times_least:            need_index = [i-1 for i in need_index_tmp]            need_index.append(len(flag_list) - 1)        else:            need_index = [i-1 for i in need_index_tmp]

        for i in need_index:            print("The %dth character %s repeat %d times" %(i - flag_list[i] + 1, test_str[i], flag_list[i] + 1))

repeat_compute(my_fasta_dict, 4)

其他

1. 个人觉得在写代码过程中多写注释，避免以后回头再看不知道怎么用
   
2. 多写函数，只需关注输入输出即可
   
3. 不要为了解决而解决问题，尽量让代码更具泛化能力。比如：原题目中要求“返回32bp-332bp、780bp-992bp的碱基”，代码中实现时不仅实现了这个功能，还通过可变参数小技巧，实现可以返回多段碱基序列
   

• 推文多平台同步发布，公众号内容食用更佳
  
• 更多内容，请关注微信公众号“生信矿工”