我有一個多 fasta 序列檔案：test.fasta
>Ara_001
MGIKGLTKLLADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGE
VTSHLQGMFNRTIRLLEAGIKPVYVFDGKPPELKRQELAKRYSKRADATADLTGAIEAGN
>Ara_002
MGIKGLTKLLADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGE
VTSHLQGMFNRTIRLLEAGIKPVYVFDGKPPELKRQELAKRYSKRADATADLTGAIEAGN
>Ara_003
MGIKGLTKLLAEHAPRAAAQRRVEDYRGRVIAIDASLSIYQFLVVVGRKGTEVLTNEAEG
LTVDCYARFVFDGEPPDLKKRELAKRSLRRDDASEDLNRAIEVGDEDSIEKFSKRTVKIT
我有另一個帶有范圍的串列檔案：range.txt
Ara_001 3 60
Ara_002 10 80
Ara_003 20 50
我想提取定義的區域。
我的預期輸出將是：
>Ara_001
KGLTKLLADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGE
VT
>Ara_002
ADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGE
VTSHLQGMFNRTIRLLEAGIKPVYVFDGKP
>Ara_003
RRVEDYRGRVIAIDASLSIYQFLVVVGRKG
我試過了：
#!/bin/bash
lines=$(awk 'END {print NR}' range.txt)
for ((a=1; a<= $lines ; a ))
do
number=$(awk -v lines=$a 'NR == lines' range.txt)
grep -v ">" test.fasta | awk -v lines=$a 'NR == lines' | cut -c$number
done;
uj5u.com熱心網友回復：
不要重新發明輪子。使用為此目的撰寫并廣泛使用的標準生物資訊學工具。在您的示例中，使用bedtools getfasta. 將您的區域檔案重新格式化為 3 列格式，然后：
bedtools getfasta -fi test.fasta -bed range.bed
安裝bedtools套件，例如，使用conda，特別是miniconda，像這樣：
conda create --name bedtools bedtools
uj5u.com熱心網友回復：
使用生物蟒：
# read ranges as dictionary
with open('range.txt') as f:
ranges = {ID: (int(start), int(stop)) for ID, start, stop in map(lambda s: s.strip().split(), f)}
# {'Ara_001': (3, 60), 'Ara_002': (10, 80), 'Ara_003': (20, 50)}
# load input fasta and slice
from Bio import SeqIO
with open ('test.fasta') as handle:
out = [r[slice(*ranges[r.id])] for r in SeqIO.parse(handle, 'fasta')]
# export sliced sequences
with open('output.fasta', 'w') as handle:
SeqIO.write(out, handle, 'fasta')
輸出檔案：
>Ara_001
KGLTKLLADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGE
>Ara_002
ADNAPSCMKEQKFESYFGRKIAVDASMSIYQFLIVVGRTGTEMLTNEAGEVTSHLQGMFN
RTIRLLEAGI
>Ara_003
RRVEDYRGRVIAIDASLSIYQFLVVVGRKG
注意。使用這個快速代碼，range.txt 中的每個序列 id 都必須有一個條目，但是很容易修改它以在沒有它的情況下使用默認行為。
uj5u.com熱心網友回復：
當我寫下最初的答案時，我誤解了這個問題。您的情況不太依賴于任何編程語言。您似乎需要 Timur Shtatland 的答案中的實用程式。安裝該實用程式或使用處理您的兩個檔案的 mozway 的 python 代碼片段。
我讀了你的問題的格式，好像你有 4 個單獨的檔案，而不是兩個。
uj5u.com熱心網友回復：
使用我正在開發的軟體包Biotite ，可以通過以下方式完成：
import biotite.sequence.io.fasta as fasta
input_file = fasta.FastaFile.read("test.fasta")
output_file = fasta.FastaFile()
with open("range.txt") as file:
for line in file.read().splitlines():
seq_id, start, stop = line.split()
start = int(start)
stop = int(stop)
output_file[seq_id] = input_file[seq_id][start : stop]
output_file.write("path/to/output.fasta")
轉載請註明出處，本文鏈接：https://www.uj5u.com/caozuo/506245.html