文章目录

大家好，我是穆易青，专注于生物信息领域的博主。近期，国产AI模型kimi-K2在生物信息代码写作方面引起了我的关注，之前我也对其他模型做过测评（2025.02.02【横向评测L1】| Deepseek与ChatGPT的生信使用体验），今天就来跟大家分享一下我对Kimi-K2的初步测评，希望能为初学者们提供一些参考。

kimi-K2模型简介

Kimi-K2 是月之暗面（Moonshot AI）于 2025 年 7 月正式发布的第二代通用大语言模型，定位为“超长上下文、多轮对话、工具调用”三位一体的智能基座。与初代相比，K2 有了很多升级：

  1. 超长上下文
    原生支持 2M token 无损窗口（≈ 300 万汉字），可在单轮对话内完整阅读《三体》三部曲并回答跨卷细节问题；通过分段压缩技术，可弹性扩展至 10M token 级别，满足整库代码检索、海量日志分析等企业场景。

  2. 推理能力跃迁
    在 GSM8K、MATH、HumanEval 等基准上相对 K1 平均提升 28%，首次在公开评测中达到 o1-preview 同级别水平；引入“自洽解码 + 过程奖励模型”双引擎，显著降低逻辑跳步与数值计算错误。

  3. 工具调用与智能体
    内建 380+ 函数调用模板，支持 REST API、SQL、Bash、数据分析（pandas/R）等即插即用；官方 Agent 框架 K-Agent 允许开发者通过 YAML 脚本把多工具链组合为“一句话工作流”。

  4. 多模态统一
    原生整合视觉编码器，可直接理解高分辨率图片、PDF、PPT、手写公式；视频输入支持最长 3 小时关键帧抽取，实现“看完整部电影并回答细节”。

请添加图片描述

本次测评主要围绕生信分析的两个核心场景展开：R语言绘图和Shell脚本编写。基础评测选择的是R语言绘图，考察模型作为生信初学者日常助手的代码写作能力。中层评测Shell脚本的分析流程编写则更加复杂，对模型的逻辑推理、代码写作、技能准确性等方面都有较高的要求。

鉴于目前kimi-K2模型线上是免费使用，因此，我选择了同样免费使用的gemini-2.5-flash和Chat-gpt4o模型进行横向对比。至于更好的收费模型，此次测评并没有纳入考虑。

场景一：R语言火山图绘制

在生物信息分析中，火山图（Volcano plot）是展示差异表达基因的常用工具。我测试了kimi-K2生成使用 ggplot2 包绘制火山图的R脚本能力。

Prompt： “写一个R脚本，使用ggplot2包画一个火山图。”

代码展示

#!/usr/bin/env Rscript
## 火山图绘制脚本 -------------------------------------------------
## 依赖：ggplot2、dplyr、readr、ggrepel
## 使用方法：
##   Rscript volcano_plot.R data.csv  # CSV 文件
##   Rscript volcano_plot.R data.tsv  # TSV 文件

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  # library(tidyr) 这个一开始没有，模型认为replace_na函数在dplyr中（可能是版本问题），实际情况是我加载了tidyr包后才能加载该函数
  library(readr)
  library(ggrepel)
})

## 1. 读取命令行参数 -------------------------------------------------
args <- commandArgs(trailingOnly = TRUE)
if (length(args) == 0) {
  stop("请提供数据文件路径，例如：Rscript volcano_plot.R result.csv")
}

file_path <- args[1]
file_ext  <- tools::file_ext(file_path)

## 2. 读入数据 -------------------------------------------------------
if (tolower(file_ext) %in% c("tsv", "txt")) {
  df <- read_tsv(file_path)
} else {
  df <- read_csv(file_path)
}

## 3. 检查必要列 -----------------------------------------------------
need_cols <- c("log2FoldChange", "padj")
if (any(!need_cols %in% names(df))) {
  stop(sprintf("数据必须包含列：%s", paste(need_cols, collapse = ", ")))
}

## 4. 预处理 ---------------------------------------------------------
df <- df %>%
  mutate(padj = replace_na(padj, 1),
         log10Padj = -log10(padj),
         threshold = case_when(
           abs(log2FoldChange) >= 1 & padj < 0.05 ~ "Up/Down",
           TRUE ~ "NotSig"
         ),
         threshold = factor(threshold, levels = c("Up/Down", "NotSig")))

## 5. 绘图 ------------------------------------------------------------
p <- ggplot(df, aes(x = log2FoldChange, y = log10Padj)) +
  geom_point(aes(color = threshold), alpha = 0.8, size = 1.5) +
  scale_color_manual(values = c("Up/Down" = "#E64B35FF", "NotSig" = "grey70")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  labs(title = "Volcano plot",
       x = expression(log[2]~Fold~Change),
       y = expression(-log[10]~italic(p)-adj)) +
  theme_bw() +
  theme(legend.title = element_blank(),
        plot.title   = element_text(hjust = 0.5))

## 6. 标记前 20 个最显著基因 -----------------------------------------
top_genes <- df %>%
  filter(threshold == "Up/Down") %>%
  arrange(padj) %>%
  slice_head(n = 20)

if (nrow(top_genes) > 0) {
  p <- p + geom_text_repel(data = top_genes,
                           aes(label = rownames(top_genes)),
                           size = 3, max.overlaps = 20)
}

## 7. 保存 & 显示 ----------------------------------------------------
outfile <- sub(paste0("\\.", file_ext, "$"), "_volcano.png", basename(file_path))
ggsave(outfile, plot = p, width = 6, height = 5, dpi = 300)
message("火山图已保存到：", normalizePath(outfile))

## 在交互式环境（RStudio）直接显示
if (interactive()) print(p)

生成图片如下：
请添加图片描述

讨论

先说结论，这个脚本有四处问题提高了初学者的使用成本。

  1. 虽然生成了绘制火山图的R脚本，但是并没有生成示例数据，无法直接运行；
  2. 因缺少加载 tidyr 包而导致 replace_na 报错；
  3. Y轴的-log[10]~italic(p)-adj这个图例没有看懂，前面计算没有问题，可能是表达有误；
  4. 上下调没有分开颜色，个人来看，绘图完成度不够。

预设标记了top20的基因，倒是比gemini将显著基因全部标记要更好一些。

相应的解决方案也有，我附在下方，就不再解释了
设置prompt：”请为我生成一个包含gene_name, log2FoldChange, padj列的示例文件DE_results.csv“
安装、加载tidyr包：

install.packages("tidyr")
library（tidyr）

另外，我还用gemini和gpt4o都绘制了图片，先看下gemini的表现：
请添加图片描述

gemini在上下调的颜色选择上和测序公司日常交付的结果最接近。当然，因为我是中文提问，gemini也很“贴心”地将图例设置成中文，所以看着有些怪。在绘制图片前，gemini自行随机生成示例数据，包含基因名、log2(倍数变化) 和 p值，因此，这里Y轴显示P值是正确的。

gemini的脚本生成后运行是报错的，因为ggplot以图层形式添加内容，最后有一个“+”号，但是后面的突出显示显著性基因的代码是处于注释状态（可选），这让R语言不熟悉的初学者很难找到问题所在。生成的脚本也不会输出图片print(p)，需要用户手动保存（在Rstudio运行比较方便）。取消掉注释后，完整的图片展示如下：
请添加图片描述

GPT-4o是我体验比较好的工具，他和gemini一样也随机生成了示例数据，同样使用的是p值。不过，由于绘制比较简单，没有考虑标记显著基因，他反而因为加载R包比较少，成为三个模型中唯一一个生成的代码无报错，直接绘图的模型。并且在我中文提示词的前提下，使用了英文图例，也优于gemini。但他也没有做到上下调的颜色区分，尽管后期可以通过优化提示很容易改进，但和gemini相比还是稍逊一筹。
请添加图片描述

对于R语言绘图，kimi-K2能提供基础代码和数据，但仍需用户进行一定的调试和优化，尤其是在包的加载和图表美观度方面。

场景二：全基因组重测序分析Shell脚本编写（质控+比对）

全基因组重测序（WGS）是生物信息学中的一项重要技术，其分析流程通常涉及多个步骤。我测试了kimi-K2在生成WGS分析Shell脚本方面的能力。

Prompt： “请为我设计一个简单的全基因组重测序分析流程的Shell脚本。流程应包括：1. 用FastQC对原始测序数据（R1.fq.gz, R2.fq.gz）进行质控；2. 用BWA-MEM将读段比对到参考基因组上；3. 用Samtools对生成的SAM文件进行格式转换、排序和去重。”

代码展示

#!/usr/bin/env bash
# -------------------------------------------------------------
# 极简全基因组重测序分析流程
# 用法:
#   bash wgs_pipeline.sh R1.fq.gz R2.fq.gz reference.fa
# 结果目录: ./results/
# -------------------------------------------------------------

set -euo pipefail

# 检查参数
if [[ $# -ne 3 ]]; then
    echo "Usage: $0 <R1.fq.gz> <R2.fq.gz> <reference.fa>"
    exit 1
fi

R1=$1
R2=$2
REF=$3
PREFIX=$(basename ${R1%%_R1*})   # 提取样本前缀，例如 sampleX
OUTDIR="results/${PREFIX}"
THREADS=8                        # 可按需调整

mkdir -p "${OUTDIR}"

echo "[$(date)] 开始分析样本: ${PREFIX}"

# -------------------------------------------------
# 1) FastQC 质控
# -------------------------------------------------
echo "[$(date)] Step1: FastQC 质控..."
fastqc -t ${THREADS} -o "${OUTDIR}/fastqc" "${R1}" "${R2}"

# -------------------------------------------------
# 2) BWA-MEM 比对
# -------------------------------------------------
echo "[$(date)] Step2: BWA-MEM 比对..."
BAM_UNSORTED="${OUTDIR}/${PREFIX}.unsorted.bam"
bwa index "${REF}" 2>/dev/null || echo "参考基因组索引已存在"
bwa mem -M -t ${THREADS} "${REF}" "${R1}" "${R2}" | \
    samtools view -bS - > "${BAM_UNSORTED}"

# -------------------------------------------------
# 3) Samtools 排序、去重
# -------------------------------------------------
echo "[$(date)] Step3: 格式转换、排序、去重..."
BAM_SORTED="${OUTDIR}/${PREFIX}.sorted.bam"
BAM_DEDUP="${OUTDIR}/${PREFIX}.dedup.bam"
DEDUP_METRICS="${OUTDIR}/${PREFIX}.dedup.metrics"

samtools sort -@ ${THREADS} -o "${BAM_SORTED}" "${BAM_UNSORTED}"
samtools index "${BAM_SORTED}"
samtools markdup -r -s -@ ${THREADS} "${BAM_SORTED}" "${BAM_DEDUP}"
samtools index "${BAM_DEDUP}"

# 清理中间文件
rm "${BAM_UNSORTED}" "${BAM_SORTED}" "${BAM_SORTED}.bai"

echo "[$(date)] 流程完成！主要结果:"
echo "  质控报告: ${OUTDIR}/fastqc/"
echo "  去重后BAM: ${BAM_DEDUP}"

总体来说，kimi-K2与Gemini 2.5flash对于复杂代码的流程编写能力相近，命令行工具的整合能力比GPT-4o更好。

讨论

脚本结构与说明： 与GPT-4o相比，生成的脚本具有清晰的用法说明和目录展示，但没有Gemini详尽；和Gemini都采用了动态传参的方式，方便用户在执行脚本时传递参数。
质控分析 (FastQC)： 在分析前会通过 echo 输出每个步骤的描述。它使用双引号将变量括起来，有效防止了变量中可能出现的特殊符号导致语法错误。
比对环节 (BWA-MEM)： K2提供了构建参考基因组索引的命令，但没有如Gemini那样，使用if判断是否存在索引文件。值得一提的是，kimi-K2采用了管道符将比对（bwa mem）与samtools view进行格式转换的步骤连接起来，跳过sam文件，直接输出未排序的BAM文件。这种做法非常专业，能够有效减少数据量较大时服务器存储的占用，体现了其实战性。
排序与去重 (Samtools)： kimi-K2的脚本将排序和去重作为独立且连续的步骤处理。它在去重步骤中增加了 -s 参数，可以输出去重报告。此外，Kimi2K在去重结束后还包含清理中间文件的步骤，这也是一个非常实用的操作。

相比之下，Gemini没有对中间文件进行删除，而GPT-4o则使用了已经被弃用的samtools rmdup命令搭建流程（且模型自己也知道），没有使用samtools markdup，显然GPT的知识库和逻辑存在一些问题。

请添加图片描述

在Shell脚本编写方面，Kimi-K2展现出了一定专业水准和实战性。尤其是在管道符的运用和中间文件清理等方面，体现了其对生物信息分析流程的深入理解。

总结与展望

作为一个模型助手，对真实世界问题的解决能力到底有多少，对已有知识的理解能否构思出新的方法，又或是跨领域知识的结合能力又有几何，我们还可以继续探索。考虑到篇幅所限，更深度的评测将另起一篇文章进行介绍。就目前来看，Kimi-K2在生物信息代码写作方面也具备独特的产品竞争力，给我留下了深刻的印象。

关于模型幻觉这个问题，这次没有专门设计prompt去测试。但我们在使用其他模型时会发现一些无中生有，移花接木的情况。有时候AI模型提供的命令行中会出现一些本身没有的参数，或者把另一个命令行工具的参数放到想执行的工具命令中，这往往出现在一些特定情境下，让用户陷入困境。

优点：

  • 专业性强： 在Shell脚本编写中，Kimi-K2展现了专业的流程设计和命令运用，例如管道符的使用和中间文件的清理，对于提升分析效率和节省资源非常有帮助。
  • 用户友好： 脚本带有明确的用法说明和动态传参，方便用户根据需求调整。
  • 生成代码即用性高： 大部分生成的代码可以直接运行，减少了用户的调试时间。

待改进点：

  • R绘图细节： 在R语言绘图方面，Kimi-K2在包的加载和图表美观度上仍有提升空间，例如自动区分上下调基因颜色。
  • 错误提示： 对于R脚本，如果能更明确地指出可能缺失的包并进行提示，会进一步提升用户体验。

总的来说，Kimi2K作为一款国产AI模型，在生物信息代码写作领域已经具备了相当的实力，特别是在Shell脚本的生成方面，其专业性和实战性令人称赞。对于生物信息初学者而言，Kimi2K无疑是一个强大的辅助工具，可以帮助大家更快地入门和理解分析流程。未来，期待Kimi2K能在R语言绘图的细节和用户体验方面继续优化，成为生物信息领域更全面、更智能的AI助手。

希望这篇文章能对大家有所启发！如果您有任何问题或想分享您的使用体验，欢迎在评论区留言交流！

附录 1：

gemini 火山图代码：

# 1. 安装和加载必要的包
# 如果你还没有安装这些包，请取消注释下面的行并运行它们
# install.packages("ggplot2")
# install.packages("dplyr") # 用于数据处理，可选但推荐

library(ggplot2)
library(dplyr) # 加载dplyr包，用于数据处理

# 2. 创建示例数据
# 在实际应用中，你会从文件中加载你的数据 (例如，CSV文件)。
# 你的数据应该至少包含三列：基因名、log2(倍数变化) 和 p值。

set.seed(123) # 设置随机种子以保证结果可重现

num_genes <- 1000 # 假设有1000个基因

sample_data <- data.frame(
  Gene = paste0("Gene_", 1:num_genes),
  log2FoldChange = rnorm(num_genes, mean = 0, sd = 1.5), # 模拟log2倍数变化，正态分布
  p_value = runif(num_genes, min = 0, max = 1) # 模拟p值，均匀分布
)

# 为了更好地模拟真实数据，我们让一些基因有较低的p值和较大的log2倍数变化。
# 这些将是我们的“显著”基因。
significant_indices <- sample(1:num_genes, 50) # 随机选择50个基因作为显著基因
sample_data$log2FoldChange[significant_indices] <- rnorm(50, mean = c(-3, 3), sd = 0.5) # 给显著基因更大的log2FC
sample_data$p_value[significant_indices] <- runif(50, min = 0, max = 0.001) # 给显著基因更小的p值

# 3. 数据预处理：添加显著性标签和-log10(p值)
# 定义显著性阈值
padj_threshold <- 0.05 # P值阈值，通常为0.05
log2fc_threshold <- 1.0 # Log2倍数变化阈值，通常为1 (表示2倍变化)

sample_data <- sample_data %>%
  mutate(
    neg_log10_p_value = -log10(p_value), # 计算-log10(p值)，用于Y轴
    Significance = case_when(
      p_value < padj_threshold & log2FoldChange > log2fc_threshold ~ "上调", # 显著上调
      p_value < padj_threshold & log2FoldChange < -log2fc_threshold ~ "下调", # 显著下调
      TRUE ~ "无显著性" # 不显著
    )
  )

# 确保“Significance”列是因子，并设置顺序，以便在图中正确排序颜色和图例。
sample_data$Significance <- factor(sample_data$Significance,
                                   levels = c("上调", "下调", "无显著性"))

# 4. 绘制火山图
ggplot(sample_data, aes(x = log2FoldChange, y = neg_log10_p_value, color = Significance)) +
  geom_point(alpha = 0.8, size = 1.5) + # 绘制散点图，设置透明度和点的大小
  scale_color_manual(values = c("上调" = "red", "下调" = "blue", "无显著性" = "grey")) + # 自定义颜色
  geom_vline(xintercept = c(-log2fc_threshold, log2fc_threshold), linetype = "dashed", color = "black") + # 添加log2FC阈值垂直虚线
  geom_hline(yintercept = -log10(padj_threshold), linetype = "dashed", color = "black") + # 添加p值阈值水平虚线
  labs(
    title = "基因表达火山图", # 图表主标题
    x = "Log2(倍数变化)", # X轴标签
    y = "-Log10(P值)", # Y轴标签
    caption = paste0("P值阈值 = ", padj_threshold, ", Log2FC阈值 = ", log2fc_threshold) # 图注
  ) +
  theme_minimal() + # 使用简洁主题
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), # 标题居中加粗
    legend.position = "right", # 图例位置
    text = element_text(family = "sans") # 设置字体为无衬线字体
  ) + #作者注释：这里如果不取消后面的代码注释，直接运行就会报错
  # 突出显示最显著的基因 (可选)
  # 如果要添加基因标签，需要安装 'ggrepel' 包。
  # install.packages("ggrepel")
  # library(ggrepel)
  # geom_text_repel(
  #   data = filter(sample_data, Significance != "无显著性" & (abs(log2FoldChange) > 2 | neg_log10_p_value > 5)), # 过滤要标记的基因
  #   aes(label = Gene),
  #   size = 3, # 标签字体大小
  #   box.padding = unit(0.35, "lines"), # 标签与点之间的距离
  #   point.padding = unit(0.5, "lines"), # 标签与点之间的最小距离
  #   max.overlaps = Inf # 允许所有标签显示，避免重叠
  # )

GPT-4o火山图代码：

# 加载必要的包
library(ggplot2)
library(dplyr)

# 模拟数据：包含基因名称、log2FoldChange 和 pvalue
set.seed(123)
data <- data.frame(
  gene = paste0("Gene", 1:1000),
  log2FoldChange = rnorm(1000, mean = 0, sd = 2),
  pvalue = runif(1000, min = 0, max = 1)
)

# 计算 -log10(pvalue)
data$negLog10Pval <- -log10(data$pvalue)

# 添加分组列，根据阈值定义显著性
data <- data %>%
  mutate(
    significance = case_when(
      pvalue < 0.05 & abs(log2FoldChange) >= 1 ~ "Significant",
      TRUE ~ "Not Significant"
    )
  )

# 绘制火山图
volcano_plot <- ggplot(data, aes(x = log2FoldChange, y = negLog10Pval)) +
  geom_point(aes(color = significance), alpha = 0.7, size = 2) +
  scale_color_manual(values = c("Significant" = "red", "Not Significant" = "grey")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "blue") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "blue") +
  labs(
    title = "Volcano Plot",
    x = "log2(Fold Change)",
    y = "-log10(p-value)",
    color = "Significance"
  ) +
  theme_minimal()

# 显示图形
print(volcano_plot)

附录2：

gemini 2.5 flash 流程代码：

#!/bin/bash

# 全基因组重测序分析流程Shell脚本

# -----------------------------------------------------------------------------
# 脚本说明：
# 这是一个用于全基因组重测序数据分析的简单Shell脚本。
# 它包括以下核心步骤：
# 1. 使用 FastQC 对原始测序数据进行质量控制。
# 2. 使用 BWA-MEM 将测序读段比对到参考基因组。
# 3. 使用 Samtools 对比对结果进行排序、格式转换 (SAM 到 BAM) 和去重。
#
# 前提条件：
# - 确保以下工具已安装并添加到系统PATH中：
#   - FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
#   - BWA (http://bio-bwa.sourceforge.net/)
#   - Samtools (http://www.htslib.org/)
# - 拥有参考基因组文件 (例如：reference.fasta)。
# - 拥有原始测序数据文件 (例如：sample_R1.fq.gz, sample_R2.fq.gz)。
#
# 使用方法：
# 1. 将脚本保存为 .sh 文件 (例如：run_reseq_pipeline.sh)。
# 2. 赋予执行权限：chmod +x run_reseq_pipeline.sh
# 3. 运行脚本：./run_reseq_pipeline.sh <参考基因组路径> <R1测序文件路径> <R2测序文件路径> <输出目录> <样本名称>
#    例如：./run_reseq_pipeline.sh /path/to/genome/ref.fasta /path/to/reads/sample_R1.fq.gz /path/to/reads/sample_R2.fq.gz ./results sample_A
# -----------------------------------------------------------------------------

# 检查输入参数数量
if [ "$#" -ne 5 ]; then
    echo "用法: $0 <参考基因组路径> <R1测序文件路径> <R2测序文件路径> <输出目录> <样本名称>"
    echo "示例: $0 /path/to/genome/ref.fasta /path/to/reads/sample_R1.fq.gz /path/to/reads/sample_R2.fq.gz ./results sample_A"
    exit 1
    
fi

# 定义变量
REFERENCE_GENOME=$1
READ1=$2
READ2=$3
OUTPUT_DIR=$4
SAMPLE_NAME=$5

# 创建输出目录，如果不存在的话
mkdir -p "${OUTPUT_DIR}"
mkdir -p "${OUTPUT_DIR}/fastqc_reports"
mkdir -p "${OUTPUT_DIR}/bwa_output"
mkdir -p "${OUTPUT_DIR}/samtools_output"

echo "----------------------------------------------------"
echo "开始全基因组重测序分析流程..."
echo "参考基因组: ${REFERENCE_GENOME}"
echo "R1 测序文件: ${READ1}"
echo "R2 测序文件: ${READ2}"
echo "输出目录: ${OUTPUT_DIR}"
echo "样本名称: ${SAMPLE_NAME}"
echo "----------------------------------------------------"

# 1. FastQC 质量控制
echo "步骤 1/3: 运行 FastQC 进行质量控制..."
fastqc -o "${OUTPUT_DIR}/fastqc_reports" "${READ1}" "${READ2}"
if [ $? -ne 0 ]; then echo "FastQC 运行失败。"; exit 1; fi
echo "FastQC 报告已生成在 ${OUTPUT_DIR}/fastqc_reports"
echo "----------------------------------------------------"

# 2. BWA-MEM 比对
echo "步骤 2/3: 运行 BWA-MEM 将读段比对到参考基因组..."

# 检查参考基因组是否已建立索引，如果没有则建立
if [ ! -f "${REFERENCE_GENOME}.sa" ]; then
    echo "BWA 索引不存在，正在建立索引..."
    bwa index "${REFERENCE_GENOME}"
    if [ $? -ne 0 ]; then echo "BWA 索引建立失败。"; exit 1; fi
    echo "BWA 索引建立完成。"
else
    echo "BWA 索引已存在，跳过索引建立。"
fi

# 执行比对
# -M: 将短的拆分比对标记为次要比对，兼容 Picard 的 MarkDuplicates
# -t: 使用的线程数 (根据你的CPU核心数调整)
bwa mem -M -t 8 "${REFERENCE_GENOME}" "${READ1}" "${READ2}" > "${OUTPUT_DIR}/bwa_output/${SAMPLE_NAME}.sam"
if [ $? -ne 0 ]; then echo "BWA-MEM 比对失败。"; exit 1; fi
echo "BWA-MEM 比对完成，SAM 文件保存为 ${OUTPUT_DIR}/bwa_output/${SAMPLE_NAME}.sam"
echo "----------------------------------------------------"

# 3. Samtools 处理 (排序、格式转换、去重)
echo "步骤 3/3: 运行 Samtools 进行排序、格式转换和去重..."

# 3.1 排序并转换为 BAM 格式
# view -bS: 将 SAM 转换为 BAM
# sort: 对 BAM 文件进行排序
samtools view -bS "${OUTPUT_DIR}/bwa_output/${SAMPLE_NAME}.sam" | samtools sort -o "${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.sorted.bam"
if [ $? -ne 0 ]; then echo "Samtools 排序或格式转换失败。"; exit 1; fi
echo "Samtools 排序并转换为 BAM 完成，文件保存为 ${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.sorted.bam"

# 3.2 去重 (使用 picard MarkDuplicates 或 samtools markdup)
# 这里使用 samtools markdup，因为它通常更方便且功能强大
# 首先需要对排序后的BAM文件建立索引，以便markdup使用
samtools index "${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.sorted.bam"
if [ $? -ne 0 ]; then echo "Samtools 索引建立失败。"; exit 1; fi
echo "Samtools 索引建立完成。"

samtools markdup -r "${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.sorted.bam" "${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.dedup.bam"
if [ $? -ne 0 ]; then echo "Samtools 去重失败。"; exit 1; fi
echo "Samtools 去重完成，文件保存为 ${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.dedup.bam"

# 建立去重后 BAM 文件的索引
samtools index "${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.dedup.bam"
if [ $? -ne 0 ]; then echo "Samtools 去重后 BAM 索引建立失败。"; exit 1; fi
echo "去重后 BAM 文件的索引建立完成。"
echo "----------------------------------------------------"

echo "全基因组重测序分析流程完成！"
echo "最终去重后的比对文件: ${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.dedup.bam"
echo "及其索引文件: ${OUTPUT_DIR}/samtools_output/${SAMPLE_NAME}.dedup.bam.bai"

gpt-4o 流程代码：

#!/bin/bash

# 设置脚本参数
REF="ref.fa"
R1="R1.fq.gz"
R2="R2.fq.gz"
THREADS=8
PREFIX="sample"  # 输出前缀，可根据样品名调整

# 创建输出目录
mkdir -p qc bam

echo "=== Step 1: FastQC 质控检查 ==="
fastqc -t $THREADS -o qc $R1 $R2

echo "=== Step 2: BWA-MEM 比对到参考基因组 ==="
# 如果参考基因组没有建索引，请先运行 bwa index $REF
bwa mem -t $THREADS $REF $R1 $R2 > ${PREFIX}.sam

echo "=== Step 3: Samtools 处理 SAM 文件 ==="
# 3.1 转换 SAM -> BAM
samtools view -@ $THREADS -bS ${PREFIX}.sam > ${PREFIX}.bam

# 3.2 排序 BAM 文件
samtools sort -@ $THREADS -o ${PREFIX}.sorted.bam ${PREFIX}.bam

# 3.3 去重（Samtools rmdup 适用于已排序的 BAM 文件）
samtools rmdup ${PREFIX}.sorted.bam bam/${PREFIX}.sorted.rmdup.bam

# 可选：索引 BAM 文件
samtools index bam/${PREFIX}.sorted.rmdup.bam

# 清理中间文件（根据需要保留）
rm ${PREFIX}.sam ${PREFIX}.bam ${PREFIX}.sorted.bam

echo "=== 全流程完成！结果保存在 qc/ 和 bam/ 目录中 ==="
# 人工智能
Logo
智能体开发者社区

中国智能体开发者社区，聚焦智能体与大模型开发，提供前沿资讯、实用工具链、开源项目及行业案例。通过技术沙龙、开发者大赛等活动，促进经验交流与协作，助力开发者快速构建创新智能应用。

更多推荐

cover

【AI智能体】Codex 制作RAG知识库应用实战操作详解

avatar 智能体开发者社区
cover

22道AI Agent 工程师必会的面试题，别啥都不会就出去找工作

avatar 智能体开发者社区
cover

基于AI agent的智能图书馆

avatar 智能体开发者社区