Linux系统下，计算并输出宏组装数据N50/N90值

#宏基因组数据组装质量评估时，N50/N90是重要的参考值#

N50值和N90值： 这两个值表示组装后的contig或scaffold的长度分布情况。N50值是指在所有contig（或scaffold）中，长度大于等于N50值的contig（或scaffold）的总长度占总长度的50%。N90值类似，只是占总长度的90%。一般来说，N50和N90值越高，组装质量越好。

N50和N90是基于组装后的contig或scaffold长度排序而得出的统计指标。以下是计算N50和N90的一般步骤：

1. 将contig或scaffold按长度从大到小排序。

2. 计算总长度。 将所有contig或scaffold的长度相加，得到总长度。

3. 计算N50。 将所有contig或scaffold的长度从大到小排列，累加长度直到累加长度达到总长度的50%。此时，累加的最后一个contig或scaffold的长度即为N50值。

4. 计算N90。 同样的方法，将长度从大到小排列的contig或scaffold长度进行累加，直到累加长度达到总长度的90%。最后一个累加的contig或scaffold的长度即为N90值。

举例说明：

假设有以下 contig（或 scaffold）的长度（单位：bp）：

500, 400, 300, 200, 100

1. 总长度 = 500 + 400 + 300 + 200 + 100 = 1500 bp

2. 计算 N50： 依次累加长度：500, 900, 1200, 1400, 1500 在累加到1500时，最后一个累加的 contig 长度是 200 bp。 因此，N50 = 200 bp。

3. 计算 N90： 依次累加长度：500, 900, 1200, 1400, 1500 在累加到1350时，最后一个累加的 contig 长度是 300 bp。 因此，N90 = 300 bp。

N50 和 N90 值越高，表示组装的 contig（或 scaffold）长度分布越均匀，组装质量越好。

---

在linux上使用bioawk计算N50的方法：

bioawk -c fastx '{ print length($seq) }' your_file.fasta | sort -rn | awk '{ sum += $1 } END { print "Total length: " sum; half = sum / 2; curr = 0; while (curr <= half) { curr += $1; if (curr > half) { print "N50: " $1; exit } } }'

这个命令的步骤如下：

1. bioawk -c fastx '{ print length($seq) }' your_file.fasta: 使用 bioawk 提取每个序列的长度，并输出到标准输出。
2. sort -rn: 对序列长度进行降序排序，以便后续计算 N50。
3. awk '{ sum += $1 } END { print "Total length: " sum; half = sum / 2; curr = 0; while (curr <= half) { curr += $1; if (curr > half) { print "N50: " $1; exit } } }': 使用 AWK 计算总长度和 N50。在此 AWK 脚本中，首先计算总长度，然后找到总长度的一半，接着迭代排序后的序列长度，直到当前累加长度大于一半的总长度，此时即可确定 N50 的值并输出。

执行这个命令后，你将会得到序列的总长度以及N50值。

要计算N90值，你可以在 AWK 脚本中稍作修改，以迭代累加序列长度直到当前累加长度超过总长度的90%。以下是对应的命令：

bioawk -c fastx '{ print length($seq) }' your_file.fasta | sort -rn | awk '{ sum += $1 } END { print "Total length: " sum; ninety = sum * 0.9; curr = 0; for (i = 1; i <= NR; i++) { curr += $i; if (curr >= ninety) { print "N90: " $i; exit } } }'

---

批量处理多个.fa文件并输出对应的N50/N90值，并带有对应的文件名作为前缀以避免混淆：

#!/bin/bash

# 循环遍历所有.fa文件
for file in *.fa; do
    echo "Processing file: $file"
    
    # 使用 bioawk 计算 N50 和 N90 值
    n_values=$(bioawk -c fastx '{ print length($seq) }' "$file" | sort -rn)
    total_length=$(echo "$n_values" | awk '{ sum += $1 } END { print sum }')
    n50=$(echo "$n_values" | awk -v total="$total_length" '{ sum += $1; if (sum >= total * 0.5) { print $1; exit } }')
    n90=$(echo "$n_values" | awk -v total="$total_length" '{ sum += $1; if (sum >= total * 0.9) { print $1; exit } }')

    # 输出文件名及 N50 和 N90 值
    echo "File: $file"
    echo "Total length: $total_length"
    echo "N50: $n50"
    echo "N90: $n90"
    echo ""
done

将以上代码保存为一个名为 calculate_N50_N90.sh 的文件，并将其放置在你的 .fa 文件所在的目录中。然后在终端中进入该目录，并执行以下命令来运行脚本：

bash calculate_N50_N90.sh