从fasta文件中提取子集

Question

我有一个 fasta 文件，它看起来像这样：

>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chrUn
ACGTGGATATTT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG
>chrUn5
TGATAGCTGTTG

我只想提取chr1, chr2, chr21,chrX以及它们的序列。所以我想要的输出是：

>chr1
ACGGTGTAGTCG
>chr2
ACGTGTATAGCT
>chr21
ACGTTGATGAAA
>chrX
GTACGGGGGTGG

如何在 unix 命令行中执行此操作？

Answer 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) 前面是 a >。

如果你要经常做这种事情，写数组很快就会变得很烦人。就个人而言，我使用以下两个脚本（我从前实验室伙伴那里继承而来）将 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 发行版的存储库中可用。例如，在基于 Debian 的系统上，您可以使用 安装它sudo apt-get install exonerate。安装后，您可以执行以下操作：

## 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

或者，您可以使用 my own 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

请注意，这是我为自己的工作编写的内容，并没有很好的文档记录或支持。尽管如此，多年来我一直在愉快地使用它。