从 fasta 文件中提取子集

Question 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

假设您有多行序列（如果您正在处理染色体，这似乎很可能），您需要首先将它们连接成一行。这使事情变得相当复杂。您可以使用awk单行：

$ 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块设置一个数组，其中每个目标序列 ID 都是一个键。其余部分只是检查每个“行”（序列），如果第一个字段位于数组 ( $i in t) 中，则它会打印当前“行”( $0)，前面带有>.

如果你经常做这种事情，写数组很快就会变得很烦人。就我个人而言，我使用下面的两个脚本（这是我从前实验室伙伴那里继承来的）将 FASTA 转换为 tbl 格式（sequence_name<TAB>sequence每行一个序列）并返回：

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

最后，如果您要处理大型序列，您可能会考虑使用可以为您完成此类操作的各种工具之一。从长远来看，它们将更快、更高效。例如，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