☰

如何获取与另一个文件中的标题行相对应的 fasta 序列

tag-icon • tag-icon text-processing bioinformatics
如何获取与另一个文件中的标题行相对应的 fasta 序列

我有一个包含标题行的文件 ( 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将文件分解为单独的记录。文件的行被读入，在换行符上连接，并且：splitcombfasta.fa\n

  • split在行首>右角之前，使用该:skip-empty选项消除空白元素。注：<?before … >Raku 中是积极的前瞻。这些记录将grep根据所需的 FastaID 进行 ped。

或者

  • comb-through 查找行首为>右角、后跟一个或多个非右角字符的元素。 Note<-[…]>是 Raku 中的负字符类（而<+[…]>是正字符类）。这些记录被grepped 为所需的 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

https://docs.raku.org/routine/grep
https://raku.org

相关内容