如何提取Cluster中的Singleton蛋白质序列ID？现有Bash脚本失效求助

我有一个包含蛋白质序列聚类的大型数据集，每个聚类以编号和多行序列条目呈现。部分序列在聚类中多次出现，部分仅出现一次（即Singleton）。需要提取每个聚类中仅出现一次的蛋白质序列ID。

数据集示例 #

>Cluster 0
0       310aa, >ref_ENST00000279791... at 100.00%
1       415aa, >ref_ENST00000641310... *
>Cluster 1
0       310aa, >ENST00000279791.590... at 100.00%
1       310aa, >ENST00000332650.693... at 100.00%
2       413aa, >ENST00000641310.590... *
3       310aa, >ENST00000279791.590... at 99.35%
4       310aa, >ENST00000332650.693... at 99.35%
>Cluster 2
0       399aa, >ENST00000641310.394... *
>Cluster 3
0       311aa, >ENST00000641081.179... at 96.14%
1       395aa, >ENST00000641310.395... *
2       311aa, >ENST00000641581.842... at 96.14%
3       311aa, >ENST00000641668.842... at 96.14%
4       311aa, >ENST00000641081.179... at 96.14%
5       299aa, >ENST00000641310.395... at 100.00%
6       311aa, >ENST00000641581.842... at 96.14%
7       311aa, >ENST00000641668.842... at 96.14%
>Cluster 4
0       380aa, >ENST00000641310.583... *
1       314aa, >ENST00000332238.915... at 95.86%
2       310aa, >ENST00000641310.583... at 97.10%
>Cluster 5
0       370aa, >ref_ENST00000314644... *
1       316aa, >ref_ENST00000642128... at 100.00%
>Cluster 6
0       367aa, >ENST00000641310.213... *
1       326aa, >ENST00000531945.112... at 96.32%
2       319aa, >ENST00000641123.112... at 98.12%
3       313aa, >ENST00000641310.213... at 99.68%
>Cluster 7
0       367aa, >ENST00000641310.284... *

期望输出 #

ENST00000641310.394
ENST00000641310.284

尝试的脚本（无法正常运行） #

#!/bin/bash

# Assuming the dataset is stored in a file called "dataset.txt"
input_file="dataset.txt"

# Loop through each line in the input file
while IFS= read -r line; do
  # Check if the line starts with ">Cluster"
  if [[ $line == ">Cluster"* ]]; then
    cluster_number=${line#>Cluster }
    cluster_number=${cluster_number//[^0-9]/}
    cluster_found=false
  fi

  # Check if the line contains a singleton protein sequence
  if [[ $line == *"... *" ]]; then
    protein_sequence=$(echo "$line" | awk -F"[>, ]" '{print $4}')
    cluster_found=true
  fi

  # Print the singleton protein sequence if a cluster was found
  if [[ $cluster_found == true ]]; then
    echo "$protein_sequence"
  fi
done < "$input_file"

问题分析与解决方案 #

原脚本的核心逻辑错误：误将带*的聚类种子序列当成Singleton，而*仅标记该序列为聚类的代表，和序列是否唯一无关。根据示例输出，我们需要提取的是所在聚类仅有一条序列的ID（Cluster2、7均为单条目聚类），以下是可行的实现：

方案：提取单条目聚类的序列ID #

#!/bin/bash
input_file="dataset.txt"

# 存储当前聚类的序列数量和对应ID
current_count=0
current_id=""

while IFS= read -r line; do
    # 遇到新聚类时，检查上一个聚类是否为单条目
    if [[ $line == ">Cluster"* ]]; then
        if [[ $current_count -eq 1 ]]; then
            echo "$current_id"
        fi
        # 重置计数器和ID
        current_count=0
        current_id=""
    # 提取序列ID（匹配>到...之间的内容）
    elif [[ $line =~ >([^.]+.[^.]+....) ]]; then
        current_id=$(echo "$line" | sed -E 's/.*>([^.]+.[^.]+....)\.\.\..*/\1/')
        ((current_count++))
    fi
done < "$input_file"

# 处理最后一个聚类
if [[ $current_count -eq 1 ]]; then
    echo "$current_id"
fi

运行该脚本将得到与示例一致的输出。

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