此篇文章關鍵給大家介紹了應用Python腳本制作獲取基因組測序指定位置編碼序列的實例詳細說明，感興趣的小伙伴值得借鑒參考一下，也希望能有一定的幫助，祝愿大家多多的發展，盡早漲薪

　　前言

　　在基因組分析中，大家常常會有這樣一個要求，便是在一個fasta文件中獲取某些編碼序列出去。有時候這種編碼序列注定是完備的編碼序列，而有時候只是為原fasta文件中某一段編碼序列中的一部分。尤其是當信息量許多時，應用人眼來挑選編碼序列會很費勁，因此這個時候大家就可以用編程邏輯去完成了。

　　比如這里在百度云盤配件中給出了某種群全基因組序列（0-refer/Bacillus_subtilis.str168.fasta），及其基因組測序gff注解文件（0-refer/Bacillus_subtilis.str168.gff）。

　　假定在這兒對于該種群進行分析，根據gff注解文檔中的基因功能敘述字段名，再加上對相關信息的查閱等，精準定位到了一部分特定的基因。

　　下面我們期待基于gff文件中對這種遺傳基因部位的描寫，在這一基因組序列fasta文件里將這種遺傳基因尋找并分離出來，獲得了一個新的fasta文件，新文件上只包括目地基因序列。

　　請使用python3撰寫1個能夠實現這個功能的腳本制作。

　　示例

　　一個示例腳本如下（可參見網盤附件“seq_select1.py”）。

　　為了實現以上目的，我們首先需要準備一個txt文件（以下稱其為list文件，示例list.txt可參見網盤附件），基于gff文件中所記錄的基因位置信息，填入類似以下的內容（列與列之間以tab分隔）。

　　#下列內容保存到list.txt

　　gene46 NC_000964.3 42917 43660+
　　NP_387934.1 NC_000964.3 59504 60070+
　　yfmC NC_000964.3 825787 826734-
　　cds821 NC_000964.3 885844 886173-

　　第1列，給所要獲取的新序列命個名稱；

　　第2列，所要獲取的序列所在原序列ID；

　　第3列，所要獲取的序列在原序列中的起始位置；

　　第4列，所要獲取的序列在原序列中的終止位置；

　　第5列，所要獲取的序列位于原序列的正鏈（+）或負鏈（-）。

　　之后根據輸入文件，即輸入fasta文件及記錄所要獲取序列位置的list文件中的內容，編輯py腳本。

　　打開fasta文件“Bacillus_subtilis.scaffolds.fasta”，使用循環逐行讀取其中的序列id及堿基序列，并將每條序列的所有堿基合并為一個字符串；將序列id及該序列合并后的堿基序列以字典的形式存儲（字典樣式{'id':'base'}）。

　　打開list文件“list.txt”，讀取其中的內容，存儲到字典中。字典的鍵為list文件中的第1列內容；字典的值為list文件中第2-5列的內容，并按tab分割得到一個列表，包含4個字符分別代表list文件中第2-5列的信息）。

　　最后根據讀取的list文件中序列位置信息，在讀取的基因組中截取目的基因序列。由于某些基因序列可能位于基因組負鏈中，需取其反向互補序列，故首先定義一個函數rev()，用于在后續調用得到反向互補序列。在輸出序列名稱時，還可選是否將該序列的位置信息一并輸出（name_detail=True/False）。

　　<pre class="r"style="overflow-wrap:break-word;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;orphans:2;text-align:start;text-indent:0px;text-transform:none;widows:2;word-spacing:0px;-webkit-text-stroke-width:0px;text-decoration-style:initial;text-decoration-color:initial;margin-top:0px;margin-bottom:10px;padding:9.5px;border-radius:4px;background-color:rgb(245,245,245);box-sizing:border-box;overflow:auto;font-size:13px;line-height:1.42857;color:rgb(51,51,51);word-break:break-all;border:1px solid rgb(204,204,204);font-family:"Times New Roman";">#!/usr/bin/env python3

　　#-*-coding:utf-8-*-
　　#初始傳遞命令
　　input_file='Bacillus_subtilis.str168.fasta'
　　list_file='list.txt'
　　output_file='gene.fasta'
　　name_detail=True
　　##讀取文件
　　#讀取基因組序列
　　seq_file={}
　　with open(input_file,'r')as input_fasta:
　　for line in input_fasta:
　　line=line.strip()
　　if line[0]=='>':
　　seq_id=line.split()[0]
　　seq_file[seq_id]=''
　　else:
　　seq_file[seq_id]+=line
　　input_fasta.close()
　　#讀取列表文件
　　list_dict={}
　　with open(list_file,'r')as list_table:
　　for line in list_table:
　　if line.strip():
　　line=line.strip().split('t')
　　list_dict[line[0]]=[line[1],int(line[2])-1,int(line[3]),line[4]]
　　list_table.close()
　　##截取序列并輸出
　　#定義函數，用于截取反向互補
　　def rev(seq):
　　base_trans={'A':'T','C':'G','T':'A','G':'C','N':'N','a':'t','c':'g','t':'a','g':'c','n':'n'}
　　rev_seq=list(reversed(seq))
　　rev_seq_list=[base_trans[k]for k in rev_seq]
　　rev_seq=''.join(rev_seq_list)
　　return(rev_seq)
　　#截取序列并輸出
　　output_fasta=open(output_file,'w')
　　for key,value in list_dict.items():
　　if name_detail:
　　print('>'+key,'['+value[0],value[1]+1,value[2],value[3]+']',file=output_fasta)
　　else:
　　print('>'+key,file=output_fasta)
　　seq=seq_file['>'+value[0]][value[1]:value[2]]
　　if value[3]=='+':
　　print(seq,file=output_fasta)
　　elif value[3]=='-':
　　seq=rev(seq)
　　print(seq,file=output_fasta)
　　output_fasta.close()</pre>

　　編輯該腳本后運行，輸出新的fasta文件“gene.fasta”，其中的序列即為我們所想要得到的目的基因序列。

　　擴展：

　　網盤附件“seq_select.py”為添加了命令傳遞行的python3腳本，可在shell中直接進行目標文件的I/O處理。該腳本可指定輸入fasta序列文件以及記錄有所需提取序列位置的列表文件，輸出的新fasta文件即為提取出的序列。

　　#!/usr/bin/env python3
　　#-*-coding:utf-8-*-
　　#導入模塊，初始傳遞命令、變量等
　　import argparse
　　parser=argparse.ArgumentParser(description='n該腳本用于在基因組特定位置截取序列，需額外輸入記錄有截取序列信息的列表文件',add_help=False,usage='npython3 seq_select.py-i[input.fasta]-o[output.fasta]-l<ul>npython3 seq_select.py--input[input.fasta]--output[output.fasta]--list<ul>')
　　required=parser.add_argument_group('必選項')
　　optional=parser.add_argument_group('可選項')
　　required.add_argument('-i','--input',metavar='[input.fasta]',help='輸入文件，fasta格式',required=True)
　　required.add_argument('-o','--output',metavar='[output.fasta]',help='輸出文件，fasta格式',required=True)
　　required.add_argument('-l','--list',metavar='<ul>',help='記錄“新序列名稱/序列所在原序列ID/序列起始位置/序列終止位置/正鏈（+）或負鏈（-）”的文件，以tab作為分隔',required=True)
　　optional.add_argument('--detail',action='store_true',help='若該參數存在，則在輸出fasta的每條序列id中展示序列在原fasta中的位置信息',required=False)
　　optional.add_argument('-h','--help',action='help',help='幫助信息')
　　args=parser.parse_args()
　　##讀取文件
　　#讀取基因組序列
　　seq_file={}
　　with open(args.input,'r')as input_fasta:
　　for line in input_fasta:
　　line=line.strip()
　　if line[0]=='>':
　　seq_id=line.split()[0]
　　seq_file[seq_id]=''
　　else:
　　seq_file[seq_id]+=line
　　input_fasta.close()
　　#讀取列表文件
　　list_dict={}
　　with open(args.list,'r')as list_file:
　　for line in list_file:
　　if line.strip():
　　line=line.strip().split('t')
　　list_dict[line[0]]=[line[1],int(line[2])-1,int(line[3]),line[4]]
　　list_file.close()
　　##截取序列并輸出
　　#定義函數，用于截取反向互補
　　def rev(seq):
　　base_trans={'A':'T','C':'G','T':'A','G':'C','a':'t','c':'g','t':'a','g':'c'}
　　rev_seq=list(reversed(seq))
　　rev_seq_list=[base_trans[k]for k in rev_seq]
　　rev_seq=''.join(rev_seq_list)
　　return(rev_seq)
　　#截取序列并輸出
　　output_fasta=open(args.output,'w')
　　for key,value in list_dict.items():
　　if args.detail:
　　print('>'+key,'['+value[0],value[1]+1,value[2],value[3]+']',file=output_fasta)
　　else:
　　print('>'+key,file=output_fasta)
　　seq=seq_file['>'+value[0]][value[1]:value[2]]
　　if value[3]=='+':
　　print(seq,file=output_fasta)
　　elif value[3]=='-':
　　seq=rev(seq)
　　print(seq,file=output_fasta)
　　output_fasta.close()

　　適用上述示例中的測試文件，運行該腳本的方式如下。

　　#python3 seq_select.py-h
　　python3 seq_select.py-i Bacillus_subtilis.str168.fasta-l list.txt-o gene.fasta--detail
　　源碼提取鏈接:https://pan.baidu.com/s/1kUhBTmpDonCskwmpNIJPkA?pwd=ih9n

　　提取碼:ih9n