как найти последовательности fasta, если строка заголовка совпадает с другим файлом

У меня есть файл строк заголовков (файл 1), а другой файл — это последовательности в формате fasta (файл 2). Мне нужны последовательности grep fasta, если строка заголовка из файла 1 совпадает с файлом 2. Пример: Файл 1:

>sp|B7UM99|TIR_ECO27
>sp|P06616|ERA_ECOLI

Файл 2:

>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
............
0
09.05.2019, 23:54
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       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

С наилучшими пожеланиями... ура, дрл

0
28.01.2020, 03:42

Дан файл 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
............
0
28.01.2020, 03:42

Теги

Похожие вопросы