網站首頁 編程語言 正文
引言
在基因組分析中，我們經常會有這么一個需求，就是在一個fasta文件中提取一些序列出來。有時這些序列是一段完整的序列，而有時僅僅為原fasta文件中某段序列的一部分。特別是當數據量很多時，使用肉眼去挑選序列會很吃力，那么這時我們就可以通過簡單的編程去實現了。
例如此處在網盤附件中給定了某物種的全基因組序列（0-refer/ Bacillus_subtilis.str168.fasta），及其基因組gff注釋文件（0-refer/ Bacillus_subtilis.str168.gff）。
假設在這里我們對該物種進行研究，通過gff注釋文件中的基因功能描述字段，加上對相關資料的查閱等，定位到了一些特定的基因。
接下來我們期望基于gff文件中對這些基因位置的描述，在全基因組序列fasta文件中將這些基因找到并提取出來，得到一個新的fasta文件，新文件中只包含目的基因序列。
請使用python3編寫一個可以實現該功能的腳本。
示例
一個示例腳本如下（可參見網盤附件“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 [list]\npython3 seq_select.py --input [input.fasta] --output [output.fasta] --list [list]')
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 = '[list]', 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
原文鏈接：https://www.jianshu.com/p/4c19389b2f43
相關推薦
- 2022-07-15 關于在Redis中使用Pipelining加速查詢的問題_Redis
- 2023-02-15 Python進行ffmpeg推流和拉流rtsp、rtmp實例詳解_python
- 2022-09-13 C++中的偽隨機數_C 語言
- 2022-02-15 H5移動端大轉盤抽獎插件, 簡單、易用、無依賴
- 2022-06-26 C++中declspec(dllexport)和declspec(dllimport)?的用法介紹_
- 2022-06-28 python神經網絡Keras構建CNN網絡訓練_python
- 2022-11-22 python正則表達式中匹配次數與貪心問題詳解(+?？*)_python
- 2022-09-14 python重寫方法和重寫特殊構造方法_python
- 最近更新
-
- window11 系統安裝 yarn
- 超詳細win安裝深度學習環境2025年最新版（
- Linux 中運行的top命令 怎么退出？
- MySQL 中decimal 的用法？ 存儲小
- get 、set 、toString 方法的使
- @Resource和 @Autowired注解
- Java基礎操作-- 運算符，流程控制 Flo
- 1. Int 和Integer 的區別，Jav
- spring @retryable不生效的一種
- Spring Security之認證信息的處理
- Spring Security之認證過濾器
- Spring Security概述快速入門
- Spring Security之配置體系
- 【SpringBoot】SpringCache
- Spring Security之基于方法配置權
- redisson分布式鎖中waittime的設
- maven:解決release錯誤：Artif
- restTemplate使用總結
- Spring Security之安全異常處理
- MybatisPlus優雅實現加密？
- Spring ioc容器與Bean的生命周期。
- 【探索SpringCloud】服務發現-Nac
- Spring Security之基于HttpR
- Redis 底層數據結構-簡單動態字符串（SD
- arthas操作spring被代理目標對象命令
- Spring中的單例模式應用詳解
- 聊聊消息隊列，發送消息的4種方式
- bootspring第三方資源配置管理
- GIT同步修改后的遠程分支
提供CDN加速