![如何获取与另一个文件中的标题行相对应的 fasta 序列](https://linux22.com/image/153304/%E5%A6%82%E4%BD%95%E8%8E%B7%E5%8F%96%E4%B8%8E%E5%8F%A6%E4%B8%80%E4%B8%AA%E6%96%87%E4%BB%B6%E4%B8%AD%E7%9A%84%E6%A0%87%E9%A2%98%E8%A1%8C%E7%9B%B8%E5%AF%B9%E5%BA%94%E7%9A%84%20fasta%20%E5%BA%8F%E5%88%97.png)
我有一个包含标题行的文件 ( file1
),另一个文件是 fasta 格式的序列 ( file2
)。如果file1
中存在来自的标题行,我想 grep fasta 序列file2
。
例子:
file1
:>sp|B7UM99|TIR_ECO27 >sp|P06616|ERA_ECOLI
file2
:>sp|B7UM99|TIR_ECO27 MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA RGGTGHLISSTGALGSRSLFSPLRNSMADS VDSRDIPGLPTNPSRLAAATSETCLLGGFE VLHDKGPLDILNTQIGPSAFRVEVQADGTH ...... >sp|P06616|ERA_ECOLI MSIDKSYCGFIAIVGRPNVGKSTLLNKLL GQKISITSRKAQTTRHRIVGIHTEGAYQAIY VDTPGLHMEEKRAINRLMNKAASSSIGDVE LVIFVVEGTRWTPDDEMVLNKLREGKAPVI ............ >sp|P0AD68|HUMAN MKAAAKTQKPKRQEEHANFISWRFALLCGC ILLALAFLLGRVAWLQVISPDMLVKEGDMR SLRVQQVSTSRGMITDRSGRPLAVSVPVKA IWADPKEVHDAGGISVGDRWKALANALNIP .............
- 所需输出
>sp|B7UM99|TIR_ECO27 MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA RGGTGHLISSTGALGSRSLFSPLRNSMADS VDSRDIPGLPTNPSRLAAATSETCLLGGFE VLHDKGPLDILNTQIGPSAFRVEVQADGTH ...... >sp|P06616|ERA_ECOLI MSIDKSYCGFIAIVGRPNVGKSTLLNKLL GQKISITSRKAQTTRHRIVGIHTEGAYQAIY VDTPGLHMEEKRAINRLMNKAASSSIGDVE LVIFVVEGTRWTPDDEMVLNKLREGKAPVI ............
答案1
给定一个 Fasta 文件,其序列行为等长,
$ cat file.fa
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLLG
QKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
>sp|P0AD68|HUMAN
MKAAAKTQKPKRQEEHANFISWRFALLCGC
ILLALAFLLGRVAWLQVISPDMLVKEGDMR
SLRVQQVSTSRGMITDRSGRPLAVSVPVKA
IWADPKEVHDAGGISVGDRWKALANALNIP
.............
和一个包含序列名称的查询文件,
$ cat query
sp|B7UM99|TIR_ECO27
sp|P06616|ERA_ECOLI
然后samtools
可以用作
$ samtools faidx file.fa -r query
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAARGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFEVLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLLGQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVELVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
答案2
有一个名为fastagrep
它的实用程序似乎可以做你想做的事。您的数据位于文件 data1 和 data2 中:
# Utility functions: print-as-echo, print-line-with-visual-space.
pe() { for _i;do printf "%s" "$_i";done; printf "\n"; }
pl() { pe;pe "-----" ;pe "$*"; }
pl " Input data file $FILE1:"
head -20 $FILE1
pl " Input data file $FILE2:"
head -20 $FILE2
pl " Expected output:"
cat $E
pl " Results:"
fastagrep -f $FILE1 $FILE2
生产:
-----
Input data file data1:
>sp|B7UM99|TIR_ECO27
>sp|P06616|ERA_ECOLI
-----
Input data file data2:
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
>sp|P0AD68|HUMAN
MKAAAKTQKPKRQEEHANFISWRFALLCGC
ILLALAFLLGRVAWLQVISPDMLVKEGDMR
SLRVQQVSTSRGMITDRSGRPLAVSVPVKA
IWADPKEVHDAGGISVGDRWKALANALNIP
-----
Expected output:
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
-----
Results:
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
这是在这样的系统上:
OS, ker|rel, machine: Linux, 3.16.0-7-amd64, x86_64
Distribution : Debian 8.11 (jessie)
perl 5.20.2
fastagrep 的详细信息:
fastagrep extract sequences from a multi-FASTA file by regex. (what)
Path : ~/bin/fastagrep
Version : fastagrep version 0.3
Length : 392 lines
Type : Perl script, ASCII text executable
Shebang : #!/usr/bin/env perl
Home : https://github.com/rec3141/rec-genome-tools (doc)
Modules : (for perl codes)
strict 1.08
warnings 1.23
Getopt::Std 1.10
IO::File 1.16
Data::Dumper 2.151_01
最美好的祝愿...干杯,drl
答案3
这个 awk 命令可能会做你想做的事:
$ cat file1 | xargs -I{} awk -v tag={} '$0==tag{p=1; print; next} /sp/{p=0}p' file2
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
该过程非常简单,几乎不需要任何解释。
答案4
使用乐(以前称为 Perl_6)
~$ raku -e 'my @query = "/path/to/fastaIDs.txt".IO.lines;
.grep(/@query/).print for lines.join("\n").split(/^^ <?before \> >/, :skip-empty);' file.fa
#OR:
~$ raku -e 'my @query = "/path/to/fastaIDs.txt".IO.lines;
.grep(/@query/).print for lines.join("\n").comb(/^^ \> <-[>]>+ /);' file.fa
#OR (more simply):
~$ raku -e 'my @query = "/path/to/fastaIDs.txt".IO.lines;
.print for lines.join("\n").comb(/^^ \> @query <-[>]>+ /);' file.fa
Raku 是 Perl 家族的一种编程语言。 Raku 是生物信息文本处理的不错选择,因为它实现了非常强大的正则表达式/语法引擎。
在顶部,从右到左阅读每个代码示例,可以使用(第一个示例)或(第二个示例)fasta.fa
将文件分解为单独的记录。文件的行被读入,在换行符上连接,并且:split
comb
fasta.fa
\n
split
在行首>
右角之前,使用该:skip-empty
选项消除空白元素。注:<?before … >
Raku 中是积极的前瞻。这些记录将grep
根据所需的 FastaID 进行 ped。
或者
comb
-through 查找行首为>
右角、后跟一个或多个非右角字符的元素。 Note<-[…]>
是 Raku 中的负字符类(而<+[…]>
是正字符类)。这些记录被grep
ped 为所需的 FastaID。在第三个示例中,连续的comb
步骤grep
被简化为单个comb
步骤。
输入示例 ( fastaIDs.txt
):
sp|B7UM99|TIR_ECO27
sp|P06616|ERA_ECOLI
输入示例 ( file.fa
):
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
>sp|P0AD68|HUMAN
MKAAAKTQKPKRQEEHANFISWRFALLCGC
ILLALAFLLGRVAWLQVISPDMLVKEGDMR
SLRVQQVSTSRGMITDRSGRPLAVSVPVKA
IWADPKEVHDAGGISVGDRWKALANALNIP
............
示例输出:
>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
注意:本质上,我们从下面的 FASTA 记录分隔代码开始,并添加路径/grep 代码来过滤感兴趣的记录:
~$ raku -e '.print for lines.join("\n").split(/^^ <?before \> >/, :skip-empty);' file.fa
#OR:
~$ raku -e '.print for lines.join("\n").comb(/^^ \> <-[>]>+ /);' file.fa