背景

最近參加了個生信的面試，記錄一下有意思的面試題。

題目描述

要求從提供的*.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的堿基”，代碼中實現時不僅實現了這個功能，還通過可變參數小技巧，實現可以返回多段堿基序列

• 推文多平臺同步發布，公眾號內容食用更佳
• 更多內容，請關注微信公眾號“生信礦工”