鸟类基因组RepeatMasker重复序列的长度统计、分类筛选最长序列及FASTA文件按长度排序方法咨询

嗨，我来帮你搞定这个问题！你已经用RepeatMasker处理了鸟类基因组的重复序列，现在需要统计每个序列的长度、按长度排序，还要找出每个重复类的最长序列，下面是一步步的具体方法，都是基于命令行的工具，简单好用：

一、生成「标题行+序列长度」的目标文件 #

首先我们可以用awk脚本快速处理你的FASTA格式重复序列文件，生成你想要的格式——每个>开头的标题行下面跟着对应序列的长度：

1. 新建一个文本文件，比如命名为fasta_length.awk，把下面的代码复制进去：

/^>/ {if (seqlen) print seqlen; print; seqlen=0; next}
{seqlen += length($0)}
END {print seqlen}

2. 在终端运行下面的命令，记得把input.fasta替换成你实际的重复序列文件名，output.txt是最终生成的结果文件名：

awk -f fasta_length.awk input.fasta > output.txt

运行完成后，打开output.txt就能看到你想要的格式啦。

二、按序列长度对文件排序 #

上面生成的文件是按原FASTA的顺序排列的，要按长度递增（或递减）排序的话，因为标题和长度是成对的，我们需要先把它们合并成单行，排序后再恢复成你想要的两行格式：

步骤1：把标题和长度合并为单行 #

新建一个名为merge_id_len.awk的脚本文件，内容如下：

/^>/ {if (len) print id "\t" len; id=$0; len=0; next}
{len += length($0)}
END {print id "\t" len}

运行命令生成单行格式的文件：

awk -f merge_id_len.awk input.fasta > id_length_pairs.txt

这个文件里每一行都是类似>rnd-4_family-127#LTR/ERV1 112的格式，用制表符分隔标题和长度。

步骤2：按长度排序 #

用sort命令对这个单行文件进行排序：

# 按长度从小到大排序
sort -k2,2n id_length_pairs.txt > sorted_id_length.txt

# 如果想要从大到小排序，用这个命令：
# sort -k2,2nr id_length_pairs.txt > sorted_id_length.txt

这里-k2,2n表示按第二列（也就是长度）的数值从小到大排序，nr则是从大到小排序。

步骤3：恢复为「标题+长度」的两行格式 #

新建split_back.awk脚本文件：

BEGIN {FS="\t"}
{print $1; print $2}

运行命令得到最终的排序后文件：

awk -f split_back.awk sorted_id_length.txt > sorted_output.txt

现在sorted_output.txt就是按你想要的顺序排列的「标题行+长度」文件了。

三、筛选每个重复类的最长序列 #

如果要找出每个重复类（比如LTR/ERV1、Unknown）的最长序列，我们可以用一个awk脚本直接处理原始FASTA文件，一次性提取每个类的最长序列及其长度：

新建longest_per_class.awk脚本文件，内容如下：

/^>/ {
    # 检查上一个序列是否是当前类的最长，更新记录
    if (current_seq != "" && current_len > max_len[current_class]) {
        max_seq[current_class] = current_id "\n" current_seq
        max_len[current_class] = current_len
    }
    # 从标题里提取重复类（#后面的部分）
    split($0, arr, "#")
    current_class = arr[2]
    current_id = $0
    current_seq = ""
    current_len = 0
    next
}
# 拼接序列并统计总长度
{
    current_seq = current_seq $0
    current_len += length($0)
}
# 处理最后一个序列
END {
    if (current_seq != "" && current_len > max_len[current_class]) {
        max_seq[current_class] = current_id "\n" current_seq
        max_len[current_class] = current_len
    }
    # 输出每个类的最长序列和长度
    for (class in max_seq) {
        print "=== 重复类：" class " ==="
        print max_seq[class]
        print "序列长度：" max_len[class]
        print ""
    }
}

运行命令：

awk -f longest_per_class.awk input.fasta > longest_repeats_per_class.txt

运行完成后，longest_repeats_per_class.txt里会按重复类分组，清晰显示每个类的最长序列及其长度，方便你查看和后续分析。

备注：内容来源于stack exchange，提问作者Jimin