如何从带有其他文件列表的fasta文件中提取用户定义区域答案

【问题标题】：How to extract user defined region from an fasta file with a list in other file如何从带有其他文件列表的fasta文件中提取用户定义区域
【发布时间】：2022-08-24 20:25:12
【问题描述】：

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

请更详细地更新问题...$number 来自哪里？ range.txt 中的 2 个数字指的是什么 - 起始位置和结束位置 - 要提取的字符串的起始位置和长度 - 其他；以及这两个数字如何跨行应用 fasta 文件？
还可以考虑查看How do I format my posts，然后使用正确的格式更新您的问题；查看您的问题历史记录，您可能还想查看What should I do when someone answers my question，然后考虑查看您的问题历史记录

标签： python bash perl bioinformatics fasta

【解决方案1】：

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

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

安装bedtools套件，例如，使用conda，特别是miniconda，像这样：

conda create --name bedtools bedtools

【讨论】：

【解决方案2】：

使用生物蟒：

# 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 都必须有一个条目，但是很容易修改它以在没有它的情况下使用默认行为。

【讨论】：

嗨 mozway，我收到以下错误：./extract_region.py：第 2 行：意外标记附近的语法错误 (' ./extract_region.py: line 2: with open('region.txt') as f:'
什么是完整的错误回溯，您是否使用了提供的相同文件？

【解决方案3】：

当我写下最初的答案时，我误解了这个问题。您的情况不太依赖于任何编程语言。您似乎需要 Timur Shtatland 的答案中的实用程序。要么安装该实用程序，要么使用 mozway 的 python 代码 sn-p 来处理两个都你的文件。

我读了你的问题的格式，好像你有 4 个单独的文件，而不是两个。

【讨论】：

亲爱的 PixelBlurb，感谢您的快速帮助。使用 ./extract.py test.fasta，我得到文件“extract.py”，第 10 行，在 <module> tmp = open(list_[0], 'r') IOError: [Errno 2] No such file 或目录：'>Ara_001'

【解决方案4】：

使用我正在开发的包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")

【讨论】：