背景

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

题目描述

要求从提供的*.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 个碱基。
[Python|生信]从Fasta文件出发获取序列的基本信息-编程之家

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的碱基”，代码中实现时不仅实现了这个功能，还通过可变参数小技巧，实现可以返回多段碱基序列
  • 推文多平台同步发布，公众号内容食用更佳
  • 更多内容，请关注微信公众号“生信矿工” [Python|生信]从Fasta文件出发获取序列的基本信息-编程之家

本文由 mdnice 多平台发布