如何用grep命令实现三蛋白文件的独有序列提取（存在-缺失分析）

前提假设 #

假设你的蛋白质序列文件是每行存储一条序列（无额外注释行）。如果是FASTA格式（带>开头的序列ID），可参考后面的预处理部分。

提取单个文件独有的序列 #

假设三个文件分别为file1.txt、file2.txt、file3.txt，执行以下命令提取各文件独有的序列：

• 仅在file1.txt中存在的序列：

  grep -vxFf file2.txt file1.txt | grep -vxFf file3.txt
  

  参数说明：

  • -v：反向匹配（输出不匹配的行）
  • -x：整行精确匹配（避免部分匹配短序列）
  • -F：将模式视为固定字符串（禁用正则，避免序列中的特殊字符干扰）
  • -f：从指定文件读取匹配模式

• 仅在file2.txt中存在的序列：

  grep -vxFf file1.txt file2.txt | grep -vxFf file3.txt
  

• 仅在file3.txt中存在的序列：

  grep -vxFf file1.txt file3.txt | grep -vxFf file2.txt
  

处理FASTA格式文件 #

如果你的文件是FASTA格式（序列ID行以>开头，序列可能跨多行），需要先将每个序列的ID和对应序列合并为单行，方便后续处理：

1. 将FASTA文件转换为制表符分隔的TSV文件（ID列 + 序列列）：

   awk '/^>/ {if (seq) print id "\t" seq; id=$0; seq=""; next} {seq = seq $0} END {print id "\t" seq}' file1.fasta > file1.tsv
   awk '/^>/ {if (seq) print id "\t" seq; id=$0; seq=""; next} {seq = seq $0} END {print id "\t" seq}' file2.fasta > file2.tsv
   awk '/^>/ {if (seq) print id "\t" seq; id=$0; seq=""; next} {seq = seq $0} END {print id "\t" seq}' file3.fasta > file3.tsv
   

2. 提取各文件的序列列作为匹配参考：

   cut -f2 file2.tsv > file2_seqs.txt
   cut -f2 file3.tsv > file3_seqs.txt
   

3. 提取仅file1.tsv中存在的序列（含ID）：

   grep -vxFf file2_seqs.txt file1.tsv | grep -vxFf file3_seqs.txt
   

内容的提问来源于stack exchange，提问作者idk_102