fastaファイルからサブセットを抽出する

Question 1

すべてのシーケンスが1行に収まる単純な例では、次のものを使用できますgrep（grepそのオプションをサポートしていない場合は--no-group-separator出力を渡すgrep -v -- '--'）。

$ grep -wEA1 --no-group-separator 'chr1|chr2|chr21|chrX' file.fa 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

複数の行のシーケンスがあると仮定すると（染色体を扱う場合はそうである可能性が高い）、まずそれを行にリンクする必要があります。これは状況がかなり複雑になった。awk1行を使用できます。

$ awk -vRS=">" 'BEGIN{t["chr1"]=1;t["chr2"]=1;t["chr21"]=1;t["chrX"]=1}
                {if($1 in t){printf ">%s",$0}}' file.fa
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

上記のスクリプトはレコード区切り文字を>（vRS=">"）に設定します。これは、「線」が定義さ>~れていることを意味します\n。次に、ブロックは、BEGIN各ターゲットシーケンスＩＤがキーである配列を設定する。残りは各「ライン」（シーケンス）をチェックし、最初のフィールドが配列（）にある場合は、前に現在の$i in t「ライン」（）を印刷します。$0>

このようなタスクを頻繁に実行すると、配列の作成がすぐに迷惑になることがあります。個人的には、FASTAをtbl形式（sequence_name<TAB>sequence1行に1つのシーケンス）に変換し、その逆に変換するために、次の2つのスクリプト（以前のラボパートナーから継承した）を使用します。

FastaToTbl

#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else 
            printf "%s\t", substr($0,2,length($0)-1)
        else 
                printf "%s", $0
}END{printf "\n"}

TblToFasta

#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld
     if (fld == NF-1){
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls){
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

ファイルに保存し$PATHて実行可能にすると、ターゲットgrepシーケンスに簡単に使用できます（上記とは異なり、複数行シーケンスに対して機能します）。

$ FastaToTbl file.fa | grep -wE 'chr1|chr2|chr21|chrX' | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

grep検索対象のファイルを渡すことができるので、拡張する方が簡単です。

$ cat ids.txt 
chr1
chr2
chr21
chrX

$ FastaToTbl file.fa | grep -wFf ids.txt | TblToFasta
>chr1 
ACGGTGTAGTCG
>chr2 
ACGTGTATAGCT
>chr21 
ACGTTGATGAAA
>chrX 
GTACGGGGGTGG

最後に、大規模なシーケンスを処理する場合は、これを実行できるさまざまなツールの1つを使用することを検討できます。長期的に見ると、より速く効率的になります。たとえばfastafetch、exonerate母音。ほとんどのLinuxディストリビューションのリポジトリで利用できます。sudo apt-get install exonerateたとえば、Debianベースのシステムでは。インストール後、次のことができます。

## Create the index    
fastaindex -f file.fa -i file.in
## Loop and retrieve your sequences
for seq in chr1 chr2 chr21 chrX; do
     fastafetch -f file.fa -i file.in -q "$seq" 
done 
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

または私のものを使用することもできます。retrieveseqs.pl、ここには他の素晴らしい機能もあります。

$ retrieveseqs.pl -h

    retrieveseqs.pl will take one or more lists of ids and extract their sequences from 
    multi FASTA file

    USAGE : retrieveseqs.pl [-viofsn] <FASTA sequence file> <desired IDs, one per line>

    -v : verbose output, print a progress indicator (a "." for every 1000 sequences processed)
    -V : as above but a "!" for every desired sequence found.
    -f : fast, takes first characters of name "(/^([^\s]*)/)" given until the first space as the search string
         make SURE that those chars are UNIQUE.
    -i : use when the ids in the id file are EXACTLY identical
         to those in the FASTA file
    -h : Show this help and exit.
    -o : will create one fasta file for each of the id files
    -s : will create one fasta file per id
    -n : means that the last arguments (after the sequence file)
         passed are a QUOTED list of the names desired.
    -u : assumes uniprot format ids (separated by |)

あなたの場合は、次のことを行います。

$ retrieveseqs.pl -fn file.fa "chr1 chr2 chr21 chrX"
[7 (4/4 found]
>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

これは私が仕事のために書いたものであり、正しく文書化またはサポートされていません。それにもかかわらず、私は何年も幸せに使用しています。

Answer