背景
最近参加了个生信的面试，记录一下有意思的面试题。
题目描述
要求从提供的*.fasta文件出发：
-
获得序列的反向互补序列，并统计信息：序列条数，碱基总数，N50，N90，GC 含量。 -
提取每条序列上 32bp-332bp、780bp-992bp 的序列。 -
统计单碱基重复 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
-
-
-
将序列拼接起来，中间添加 300 个 N 从而形成一条序列，保存的结果文件中每行 60 个碱基。
![[Python|生信]从Fasta文件出发获取序列的基本信息-编程知识网](https://img.internicdomainnames.com/2023/11/20231111072920299.png "[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)
其他
-
个人觉得在写代码过程中多写注释，避免以后回头再看不知道怎么用 -
多写函数，只需关注输入输出即可 -
不要为了解决而解决问题，尽量让代码更具泛化能力。比如：原题目中要求“返回32bp-332bp、780bp-992bp的碱基”，代码中实现时不仅实现了这个功能，还通过可变参数小技巧，实现可以返回多段碱基序列
-
推文多平台同步发布，公众号内容食用更佳 -
更多内容，请关注微信公众号“生信矿工” ![[Python|生信]从Fasta文件出发获取序列的基本信息-编程知识网](https://img.internicdomainnames.com/2023/11/20231111152921706.jpeg "[Python|生信]从Fasta文件出发获取序列的基本信息-编程知识网")
本文由 mdnice 多平台发布
支付宝扫一扫
微信扫一扫
.png)
.png)
.png)
.png)
