☰

他のファイルのヘッダー行に対応するfastaシーケンスを取得する方法

tag-icon • tag-icon text-processing bioinformatics
他のファイルのヘッダー行に対応するfastaシーケンスを取得する方法

file1ヘッダー行を含むファイル（）とfasta形式のシーケンスである他のファイル（）がありますfile2。file1のヘッダ行がある場合、fastaシーケンスをgrepしたいと思います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は非常に強力な正規表現/文法エンジンを実装しているため、バイオインフォマティクスのテキスト処理に適した選択です。

上部の各コード例を右から左に読み込んだら、(最初の例)または(2番目の例)を使用してfasta.faファイルを別々のレコードに分割できます.ファイルの行を読み、改行文字で連結し、次の操作を行います。splitcombfasta.fa\n

  • split行の右端の前にある空白要素を削除するには、>このオプションを使用します。:skip-empty注：<?before … >Rakuの肯定的な予測。これらの履歴は、grep必須のFastaIDに基づいて修正されます。

または

  • comb>-throughは右端で始まり、その後に右以外の文字が1つ以上来る要素を探します。 Note<-[…]>は Raku の否定的な文字クラス (<+[…]>しかし . これらのレコードはgrep必須 FastaID にリンクされています。 3 番目の例では、連続した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

関連情報