2024.11.27【流程搭建L2】3句话搭建工业级生信流程（以RNA-seq为例）

在生物信息学领域，RNA-seq技术已成为研究基因表达变化的重要工具。为了高效地处理和分析RNA-seq数据，构建一个高效、可复用的分析流程至关重要。本文将详细介绍如何使用Cursor工具结合Snakemake语法框架来搭建一个RNA-seq生物信息分析流程，从理论到实践，为生物信息初学者提供一个全面的学习资源。

1.Cursor简介

Cursor是基于VSCode技术构建的AI代码编辑器，由Kite公司开发，旨在通过集成人工智能技术来提高开发者的编程效率和体验。

Cursor是基于VSCode技术构建的高级代码编辑器。它从VSCode的代码库派生而来，继承了VSCode的所有基础功能和用户界面，布局和操作基本一致，使得VSCode用户可以无缝迁移到Cursor，无需重新学习新的操作习惯。

最大的不同在于它内置了AI进行代码协作的功能，这是对VSCode的扩展和增强。Cursor的AI功能提供了代码补全、对话窗口和代码生成重写等辅助编程功能，这些在VSCode中通常需要通过插件如GitHub Copilot来实现。尽管Cursor基于VSCode，但它是一个独立的集成开发环境（IDE），提供了更深层次的AI集成和控制，这是作为VSCode插件无法实现的。

2.Cursor快速使用

要快速上手使用Cursor，可以按照以下步骤进行：

1. 下载和安装：

   • 访问Cursor官方网站cursor.com下载Cursor编辑器。
   • 根据你的操作系统（Windows、Mac或Linux）完成安装，我们以Windows为例。

2. 账号注册与登录：

   • 在第一次打开Cursor前，请先打开梯子。Cursor支持google与github账号信息登录。新用户可以享受2周的免费Pro计划。

3. 熟悉界面：

   • 如果你熟悉VS Code，使用Cursor会让你感到得心应手。你可以导入现有的VS Code配置到Cursor中。初学者也不用担心，使用默认设置也同样完成任务。

4. 核心功能：

   • 智能代码补全（Cursor Tab）：在编码时，Cursor会以灰色文本显示建议。按Tab接受建议，或按Esc取消。
   • 代码生成和编辑（Cmd K）：使用快捷键Ctrl/⌘ + K呼出对话框，输入你的要求，按Enter键开始编译代码。
   • Chat功能（重点！！！）（Cursor Chat）：使用Ctrl/⌘ + L开启Chat。如果鼠标选中了代码，使用Ctrl/⌘ + L会自动把代码插入到Chat中。

5. 多文件编辑：

   • 使用Composer进行多文件编辑，这适用于需要在代码库的多个部分进行修改的任务。

3.流程搭建（以RNA-seq为例）

RNA-seq分析流程概述

正常的RNA-seq分析流程包括数据前处理（qc）、比对到参考序列（align）、转录本定量（express）以及差异表达分析等步骤。使用AI编辑器后，读者将不必注意每个分析流程的步骤，而是需要了解提示工程的输入细节，我会介绍如何通过几句话就可以实现流程框架的搭建。

3.1 具体化（Specificity）

使用Ctrl/⌘ + L开启Chat后，我们可以向AI提出搭建需求。但是提示内容将会决定输出脚本的质量。

第一句话就是让Chat搭建一个流程，但是需求需要具体化。这里举个例子：

请帮我搭建一个rna-seq流程

和

请基于snakemake语言框架，帮我搭建一个rna-seq流程，包含质控，比对，定量，差异分析，注释分析等内容

区别非常大。经过测试，前者只输出到定量的分析，并且各个软件是单独的命令，无法连在一起。而后者可以通过snakemake框架将各环节进行连接。

# 导入必要的包
import os
from pathlib import Path

# 加载配置文件
configfile: "config.yaml"

# 获取所有样本名称
SAMPLES = config["samples"]

# 定义最终输出文件
rule all:
    input:
        "results/multiqc/multiqc_report.html",
        expand("results/counts/{sample}_counts.txt", sample=SAMPLES)

# FastQC质控
rule fastqc:
    input:
        "raw_data/{sample}.fastq.gz"
    output:
        html="results/fastqc/{sample}_fastqc.html",
        zip="results/fastqc/{sample}_fastqc.zip"
    params:
        outdir="results/fastqc"
    threads: 2
    shell:
        """
        fastqc {input} \
            --outdir {params.outdir} \
            --threads {threads}
        """

# MultiQC汇总报告
rule multiqc:
    input:
        expand("results/fastqc/{sample}_fastqc.zip", sample=SAMPLES)
    output:
        "results/multiqc/multiqc_report.html"
    params:
        fastqc_dir="results/fastqc",
        outdir="results/multiqc"
    shell:
        """
        multiqc {params.fastqc_dir} \
            -o {params.outdir}
        """

# STAR比对
rule star_align:
    input:
        reads="raw_data/{sample}.fastq.gz"
    output:
        bam="results/aligned/{sample}.Aligned.sortedByCoord.out.bam"
    params:
        index=config["star_index"],
        outprefix="results/aligned/{sample}."
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
            --genomeDir {params.index} \
            --readFilesIn {input.reads} \
            --readFilesCommand zcat \
            --outFileNamePrefix {params.outprefix} \
            --outSAMtype BAM SortedByCoordinate
        """

# featureCounts计数
rule featurecounts:
    input:
        bam="results/aligned/{sample}.Aligned.sortedByCoord.out.bam"
    output:
        counts="results/counts/{sample}_counts.txt"
    params:
        gtf=config["gtf"]
    threads: 4
    shell:
        """
        featureCounts -T {threads} \
            -a {params.gtf} \
            -o {output.counts} \
            {input.bam}
        """

Cursor多文件处理的特性，也帮开发人员生成了流程相应的configs.yaml，大大缩短了开发时间。

# 文件命名模式
naming_pattern:
    paired_end:
        r1: "{sample}_R1"  # 双端测序R1命名模式
        r2: "{sample}_R2"  # 双端测序R2命名模式
    single_end:
        r1: "{sample}"     # 单端测序命名模式

# 参考基因组和注释
ref:
    star_index: "/path/to/star/index"
    gtf: "/path/to/annotation.gtf"
    genome_fasta: "/path/to/genome.fa"

# 文件设置
file_params:
    fastq_ext: ".fastq.gz"
    results_dir: "results"
    
# 质控参数
qc_params:
    trim_quality: 20
    min_length: 36
    adapter: "AGATCGGAAGAG"
    # 双端测序特定参数
    pe_adapter1: "AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"
    pe_adapter2: "AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT"
    # 单端测序特定参数
    se_adapter: "AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"

# STAR参数
star_params:
    sjdbOverhang: 100
    outFilterMismatchNmax: 999
    outFilterMismatchNoverReadLmax: 0.04
    outFilterMultimapNmax: 20
    alignSJoverhangMin: 8
    alignSJDBoverhangMin: 1
    outFilterType: "BySJout"
    alignIntronMin: 20
    alignIntronMax: 1000000
    alignMatesGapMax: 1000000

# featureCounts参数
featurecounts_params:
    feature_type: "exon"
    attribute_type: "gene_id"
    strand_specific: 0  # 0=unstranded, 1=stranded, 2=reversely stranded
    min_overlap: 1

注：该文件后期经过一些改动，部分参数可能与上述流程代码无法对应。

确定CHAT写的内容没有问题，可以直接点击accept all进行引入。如果前面已经有大量代码的话，建议先选择save all，这样如果新生成的代码有问题，方便后撤一步。

3.2 展示范例（Few-shot prompting）

尽管通过snakemake连接起来了，并且也有config.yaml作为参数和数据输入。但作为工业级生信流程，我们是需要将每个大的分析步骤进行模块化设计，即单独写成一个文件，方便执行一些与当前分析结果相关的统计代码或者绘图脚本。在自身没有对流程建立整体构思框架前，参考其他同行的写作框架就显得非常重要。因此，第二句话就是引用github上写作框架比较完备的链接喂给AIchat，使其根据引用框架对现有代码进行多文件协作处理。

第二句话：

参考@https://github.com/zhxiaokang/RASflow/blob/master/workflow/quality_control.rules ，将各个分析步骤单独写一个文本，并在Snakefile中调用

根据我的引用，Chat明白了要将流程拆分成多个rules文件这个深层含义。
第一步，创建一个workflow汇总文件夹；
第二步，并将流程拆分成多个rules放入其中；
这里以quality_control.rules为例进行展示：

import pandas as pd

# FastQC质控
rule fastqc:
    input:
        unpack(get_fastq)
    output:
        html=f"{RESULTS_DIR}/fastqc/{{sample}}_fastqc.html",
        zip=f"{RESULTS_DIR}/fastqc/{{sample}}_fastqc.zip"
    params:
        outdir=f"{RESULTS_DIR}/fastqc"
    log:
        "logs/fastqc/{sample}.log"
    threads: 2
    shell:
        """
        mkdir -p {params.outdir}
        fastqc {input} \
            --outdir {params.outdir} \
            --threads {threads} \
            2> {log}
        """

# Trimmomatic质控
rule trimmomatic:
    input:
        unpack(get_fastq)
    output:
        r1=f"{RESULTS_DIR}/trimmed/{{sample}}_R1_trimmed.fastq.gz",
        r2=f"{RESULTS_DIR}/trimmed/{{sample}}_R2_trimmed.fastq.gz" if config["samples"]["{sample}"]["type"] == "paired" else []
    params:
        quality=config["qc_params"]["trim_quality"],
        minlen=config["qc_params"]["min_length"],
        adapter=config["qc_params"]["adapter"]
    log:
        "logs/trimmomatic/{sample}.log"
    threads: 4
    run:
        if config["samples"][wildcards.sample]["type"] == "paired":
            shell(
                """
                trimmomatic PE -threads {threads} \
                    {input.r1} {input.r2} \
                    {output.r1} {output.r1}.unpaired \
                    {output.r2} {output.r2}.unpaired \
                    ILLUMINACLIP:{params.adapter}:2:30:10 \
                    LEADING:{params.quality} \
                    TRAILING:{params.quality} \
                    SLIDINGWINDOW:4:{params.quality} \
                    MINLEN:{params.minlen} \
                    2> {log}
                """
            )
        else:
            shell(
                """
                trimmomatic SE -threads {threads} \
                    {input.r1} \
                    {output.r1} \
                    ILLUMINACLIP:{params.adapter}:2:30:10 \
                    LEADING:{params.quality} \
                    TRAILING:{params.quality} \
                    SLIDINGWINDOW:4:{params.quality} \
                    MINLEN:{params.minlen} \
                    2> {log}
                """
            )

这个比我自己写的参数要详细得多，看着也更舒服。

第三步；修改Snakefile源文件改为调用各个rules文件,不再展示完整流程；

import os
from pathlib import Path

# 加载配置文件
configfile: "config.yaml"

# 获取样本名称
SAMPLES = list(config["samples"].keys())
RESULTS_DIR = config["file_params"]["results_dir"]

def get_fastq(wildcards):
    """根据样本类型返回fastq文件（这里没调教好CHAT，对单双端数据做了一个不太友好的判定）"""
    if config["samples"][wildcards.sample]["type"] == "paired":
        return {
            "r1": f"raw_data/{config['samples'][wildcards.sample]['r1']}{config['file_params']['fastq_ext']}", 
            "r2": f"raw_data/{config['samples'][wildcards.sample]['r2']}{config['file_params']['fastq_ext']}"
        }
    else:
        return {
            "r1": f"raw_data/{config['samples'][wildcards.sample]['r1']}{config['file_params']['fastq_ext']}"
        }

# 定义最终输出文件
rule all:
    input:
        f"{RESULTS_DIR}/multiqc/multiqc_report.html",
        expand(f"{RESULTS_DIR}/counts/{{sample}}_counts.txt", sample=SAMPLES),
        f"{RESULTS_DIR}/qc/fastqc_report.html",
        f"{RESULTS_DIR}/star/alignment_stats.txt"

# 导入各个步骤的规则
include: "workflow/quality_control.rules"
include: "workflow/alignment.rules"
include: "workflow/quantification.rules"
include: "workflow/report.rules"

最后，为了方便用户了解代码修改后的情况，还专门输出代码库的结构目录：

├── config.yaml
├── Snakefile
└── workflow/
    ├── quality_control.rules
    ├── alignment.rules
    ├── quantification.rules
    └── report.rules

以上内容全是由Cursor自动生成，看着输出的内容很多，但用户确实只需要输入两句话就可以生成这么多高质量的代码框架。

3.3 上下文提示（context）

优化是最重要的一步，而由于Cursor给予了CHAT更高的权限，对于整个代码的上下文提示都比copliot做得更好。
前面已经介绍到，在编码时，Cursor会以灰色文本显示建议。按Tab接受建议，或按Esc取消，这个功能就是基于上下文代码进行补足。另外，在对话框中，新的对话会保留历史记录，其中也包括了生成的代码。

这里举个简单的例子，比如我们为了方便输入数据信息，要单独创建一个metadata.tsv文件，Cursor会基于上下文帮助我们创建这个文件，并给出格式参考。

第三句话

请参考@https://github.com/zhxiaokang/RASflow/blob/master/configs/metadata.tsv 单独写一个文档用于样品和分组信息输入，并在config.yaml中调用

注：考虑到第二句话AI没有自己设置数据配置文件，我重新引用了链接，并给出了提示。

Chat回答如图：

metadata.tsv格式如下：

sample_id	type	fq1	fq2	condition	batch
sample1	paired	sample1_R1	sample1_R2	control	batch1
sample2	single	sample2	NA	treatment	batch1
sample3	paired	sample3_R1	sample3_R2	treatment	batch2

数据文件，参数文件，流程模块，作为分析的主要框架基本就算搭起来了。后面还需要一些脚本库添加统计脚本，参考基因组文件汇总统计这些，都可以让Cursor帮忙整理和编写。

3.4 流程的使用方法

这个流程搭建完成后，可能你需要上传github文档，或者给同事使用，都应该写一个安装使用说明，这个cursor也帮你想到了。

3.4.1 创建环境

# 创建conda环境
conda create -n rna_seq python=3.8
conda activate rna_seq

# 安装依赖
conda install -c bioconda snakemake
conda install -c bioconda fastqc trimmomatic star subread multiqc

3.4.2 运行流程

1. 准备配置文件

# 编辑样本信息
vim metadata.tsv

# 修改配置参数
vim config.yaml

2. 运行程序

snakemake -j 4 --use-conda

3. 输出结果

results/
├── fastqc/          # 质控报告
├── trimmed/         # 过滤后数据
├── aligned/         # 比对结果
├── counts/          # 基因表达量
└── multiqc/         # 汇总报告

4.结论与展望

看到这里的读者已经知道，题目中说的三句话其实是指提示工程中的三个提示要素，理解并用好这些要素可以显著提升我们的生信分析效率。当然，这并不意味着CHAT能够完全替代人工，因为代码最后仍然需要由人来检查bug，和优化细节，如软件更换，参数调整，多因素判断（单\双端,真\原核）等要求对CHAT仍然具有一定挑战性，也对用户的知识储备有了更高的要求。因此，对于初学者而言，python和snakemake的语法基础仍然是需要掌握的。

我始终相信，无论技术如何进步，工具永远不会替代人类，只有会用工具的人替代不会工具的人。祖先如此，后代亦如此

参考文献

[1] Snakemake workflow for RNA-seq analysis. GitHub.

附录

Snakemake语法快速参考：snakemake官方文档

Cursor工具使用指南：Cursor官方文档

---

🌟 非常感谢您抽出宝贵的时间阅读我的文章。如果您觉得这篇文章对您有所帮助，或者激发了您对生物信息学的兴趣，我诚挚地邀请您：

👍 点赞这篇文章，让更多人看到我们共同的热爱和追求。

🔔 关注我的账号，不错过每一次知识的分享和探索的旅程。

📢 您的每一个点赞和关注都是对我最大的支持和鼓励，也是推动我继续创作优质内容的动力。

📚 我承诺，将持续为您带来深度与广度兼具的生物信息学内容，让我们一起在知识的海洋中遨游，发现更多未知的奇迹。

💌 如果您有任何问题或想要进一步交流，欢迎在评论区留言，我会尽快回复您。

🌐 点击下方的微信名片，获取本书资料，加入交流群，与志同道合的朋友们一起探讨、学习和成长。