Linux系统下，提取.fasta文件中序列长度＞n的序列(举例：sqlen＞1000)

#在宏基因组/宏转录组数据进行组装后，常需要去除短片段，筛选出较长的片段以供后续分析#

---

在Linux系统中，您可以使用一些文本处理工具来提取长度大于n的序列。其中，常用的工具之一是awk命令，它可以用于处理文本文件并提取符合条件的行。

以下是在Linux系统中使用awk命令来提取长度大于1000的序列的示例命令： 

awk '/^>/ {if (seqlen > 1000) {if (seqname != "") {print seqname; print seq} } seq=""; seqlen=0; seqname=$0; next} {seq = seq $0; seqlen += length($0)} END {if (seqlen > 1000) {if (seqname != "") {print seqname; print seq} } }' final_assembled.fasta > extracted_sequences.fasta

详细的：

awk '/^>/ {
    if (seqlen > 1000) {
        if (seqname != "") {
            print seqname;  # 打印上一个序列的名称
            print seq;      # 打印上一个序列的内容
        }
    }
    seq = "";            # 重置当前序列内容为空
    seqlen = 0;          # 重置当前序列长度为0
    seqname = $0;        # 保存当前序列名称行，包括 ">" 符号
    next                 # 跳过后续的处理并继续到下一行
}
{
    seq = seq $0;         # 将当前行的内容追加到当前序列内容中
    seqlen += length($0);  # 计算当前序列的长度
}
END {
    if (seqlen > 1000) {
        if (seqname != "") {
            print seqname;   # 打印最后一个序列的名称
            print seq;       # 打印最后一个序列的内容
        }
    }
}' final_contigs.fasta > extracted_sequences.fasta

代码的执行流程如下：

1. awk 是一个用于文本处理的命令行工具。我们用 /^>/ 来匹配以 ">" 开头的行，即序列名称行。

2. 当遇到以 ">" 开头的行时，表示遇到了序列名称行。接下来执行大括号内的代码块。

3. 在代码块中，首先执行 if (seqlen > 1000) 条件判断，检查上一个序列的长度是否大于1000。seqlen 是一个变量，用于保存当前序列的长度。

4. 如果上一个序列的长度大于1000，则进入内层的 if (seqname != "") 条件判断。seqname 是一个变量，用于保存上一个序列的名称行。

5. 如果上一个序列的名称行非空，则执行 print seqname; print seq;，打印上一个序列的名称和内容。seq 是一个变量，用于保存上一个序列的内容。

6. 然后，执行 seq = ""; seqlen = 0;，将当前序列的内容和长度重置为空，为下一个序列做准备。

7. 在序列名称行的处理完成后，使用 next 命令跳过后续的处理，并继续到下一行。

8. 对于不是以 ">" 开头的行，表示是序列的内容行。程序将当前行的内容追加到 seq 变量中，同时增加当前序列的长度 seqlen。

9. 当处理完所有行后，在 END 部分进行最后的处理。这部分代码在整个文件处理完毕后执行一次。

10. 在 END 部分，首先执行 if (seqlen > 1000) 条件判断，检查最后一个序列的长度是否大于1000。

11. 如果最后一个序列的长度大于1000，则进入内层的 if (seqname != "") 条件判断。

12. 如果最后一个序列的名称行非空，则执行 print seqname; print seq;，打印最后一个序列的名称和内容。

13. 最终，通过重定向操作符 >，将提取的具有序列名称的序列保存到 "extracted_sequences.fasta" 文件中