如何比较两个文件，使用第一个文件中的 ID 提取特定序列并打印输出

Question 1

您可以编写脚本来直接执行此操作，但我建议您使用下面的两个小脚本。我花了很多年时间处理这类数据，它们非常有用。如果您打算进行大量 FASTA 文件操作，我建议您使用它们。

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

将上述脚本保存在您的$HOME/bin目录中（如果不存在，请创建它，然后注销并重新登录以将其添加到您的$PATH）。然后使用使脚本可执行chmod a+x ~/bin/TblToFasta ~/bin/FastaToTbl。

这FastaToTbl会将 FASTA 序列文件转换为 tbl 格式（序列 ID，后跟制表符和序列，全部在一行中）。一旦它们变成 tbl 格式，并且 ID 和序列在同一行上，您就可以使用它来grep搜索您的 ID。然后，您将输出通过TblToFastatoi 转换回 FASTA：

$ FastaToTbl file2 | grep -wFf file1 | TblToFasta
>AAX50250.1 cds:annotated chromosome:ASM1212v1:Chromosome:1:1770:1 gene:CTA_0001.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical protein 
ATGAGCATCAGGGGAGTAGGAGGCAACGGGAATAGTCGAATCCCTTCTCATAATGGGGAT
GGATCGAATCGCAGAAGTCAAAATACGAAGAATAAAGTTGAAGATCGAGTTCGTTCTCTA
TATTCATCTCGTAGTAACGAAAATAGAGAATCTCCTTATGCAGTAGTAGACGTCAGCTCT
>AAX50251.1 cds:annotated chromosome:ASM1212v1:Chromosome:1915:2187:-1 gene:CTA_0002.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:hypothetical membrane associated protein 
ATGCTTTGTAAAGTTTGTAGAGGATTATCTTCTCTTATTGTTGTTCTCGGAGCTATAAAC
ACTGGAATTTTAGGAGTAACAGGGTATAAGGTAAACCTACTTACTCACCTGCTTGGTGAA
GGAACCATGTGGACACAAGCAGCTTATGTAGTAACGGGAATCGCTGGGGTTATGGTCTGC
>AAX50252.1 cds:annotated chromosome:ASM1212v1:Chromosome:2388:2690:1 gene:CTA_0003.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:gatC description:glutamyl-tRNA(Gln) amidotransferase subunit C 
ATGACAGAGTCATATGTAAACAAAGAAGAAATCATCTCTTTAGCAAAGAATGCTGCATTG
GAGTTGGAAGATGCCCACGTGGAAGAGTTCGTAACATCTATGAATGACGTCATTGCTTTA
ATGCAGGAAGTAATCGCGATAGATATTTCGGATATCATTCTTGAAGCTACAGTGCATCAT
TTCGTTGGTCCAGAGGATCTTAGAGAAGACATGGTGACTTCGGATTTTACTCAAGAAGAA
TAG

使用的选项grep是：

   -F, --fixed-strings
          Interpret PATTERN as a list of fixed strings  (instead  of  regular
          expressions), separated by newlines, any of which is to be matched.

   -f FILE, --file=FILE
          Obtain  patterns  from  FILE, one per line.  If this option is used
          multiple times or is combined with the -e (--regexp) option, search
          for all patterns given.  The empty file contains zero patterns, and
          therefore matches nothing.
   -w, --word-regexp
          Select only those lines containing matches that form  whole  words.
          The  test  is  that  the  matching  substring must either be at the
          beginning of the  line,  or  preceded  by  a  non-word  constituent
          character.   Similarly, it must be either at the end of the line or
          followed by a  non-word  constituent  character.   Word-constituent
          characters are letters, digits, and the underscore.

所以，-f让我们给它一个文件来读取搜索 [模式，告诉-F它将模式视为字符串而不是正则表达式（否则，如果您的 ID 包含.它们经常包含的内容，那将被视为“匹配任何字符”）并且确保-w我们只匹配整个 ID，所以foo不会匹配foobar。

另外，这里有一个快速而肮脏的 perl 单行程序来做你想做的事情：

perl -e 'open(F,"File1.txt");while(<F>){/(\S+)/; $k{$1}++}; while(<>){if(/>\s*(\S+?)(\.| )/){if($k{$1}){$k=1}else{$k=0}; } print if $k==1;}' File2.fa

或者，稍微不那么脏一点：

perl -e 'open(F,"File1.txt");
         while(<F>){
            /(\S+)/; 
            $k{$1}++
         } 
         while(<>){
            if(/>\s*(\S+?)(\.| )/){
                if($k{$1}){$k=1}
                else{$k=0}
            } 
            print if $k==1
        }' File2.fa

Answer