perl 如何从fasta文件中提取用户定义的区域,并在其他文件中列出

3okqufwl  于 2022-11-15  发布在  Perl
关注(0)|答案(4)|浏览(172)

我有一个多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;
tuwxkamq

tuwxkamq1#

不要重复发明。使用为此目的编写并广泛使用的标准生物信息学工具。在您的示例中,使用bedtools getfasta。将区域文件重新格式化为3列床格式,然后:

bedtools getfasta -fi test.fasta -bed range.bed

安装bedtools套件,例如,使用conda,特别是miniconda,如下所示:

conda create --name bedtools bedtools
toe95027

toe950272#

使用biopython:

# 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都必须有一个条目,但是,如果没有它,修改它以使用默认行为非常容易。*
j8ag8udp

j8ag8udp3#

当我写我的第一个答案时,我误解了这个问题。你的情况并不太依赖于任何编程语言。你似乎需要Timur Shtatland的答案中的实用程序。要么安装该实用程序,要么使用mozway的python代码片段来处理 * 两个 * 你的文件。
我读你的问题的格式,如果你有4个单独的文件,而不是两个。

lrpiutwd

lrpiutwd4#

对于我正在开发的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")

相关问题