计算机与编程基础
================

本章将系统性地介绍生物信息学研究中必需的计算机与编程基础知识。对于刚接触生物信息学的研究人员来说，掌握这些基础知识是开展后续数据分析工作的前提条件。生物信息学本质上是一门交叉学科，它需要研究者同时具备生物学知识和计算机技能。在实际工作中，无论是处理高通量测序数据、构建分析流程，还是开发新的算法工具，都离不开扎实的编程基础和对计算环境的深入理解。

现代生物信息学工作流程通常涉及多个环节，从原始数据的获取与预处理，到中间结果的存储与管理，再到最终的可视化与解读，每个环节都需要相应的技术支撑。Linux操作系统因其稳定性和强大的命令行工具生态，成为生物信息学分析的主流平台。Shell脚本编程能力使我们能够自动化处理重复性任务，提高工作效率。Python语言凭借其简洁的语法和丰富的科学计算库，在数据处理和分析中扮演着核心角色。数据库技术则帮助我们高效地存储和查询大规模生物学数据。版本控制系统确保我们的代码和分析流程可追溯、可重复。容器化技术解决了软件环境配置的难题，使得分析流程在不同计算平台间具有更好的可移植性。

本章内容按照从基础到进阶的顺序组织，首先介绍Linux操作系统的基本概念和常用命令，然后逐步深入到Shell脚本编程、Python编程、数据库操作等主题。每个知识点都配有实际可运行的代码示例，读者可以在自己的计算机上实践这些操作，通过动手实践来加深理解。需要特别说明的是，本章专注于计算机技术层面，不涉及具体的生物学背景知识，这些生物学背景将在后续章节中结合具体应用场景进行讲解。

.. _linux_os:

Linux 操作系统
--------------

Linux发行版选择
~~~~~~~~~~~~~~~

Linux操作系统并非一个单一的产品，而是由Linux内核配合各种系统工具、应用程序和桌面环境组成的完整系统。不同的组织和个人将这些组件以不同方式组合，形成了众多的Linux发行版。对于生物信息学工作者来说，选择合适的Linux发行版是开始工作的第一步。

Ubuntu是目前使用最广泛的Linux发行版之一，它以易用性和丰富的软件仓库著称。Ubuntu每六个月发布一个新版本，每两年发布一个长期支持版本，长期支持版本提供五年的安全更新和维护。对于生物信息学初学者，Ubuntu是一个很好的起点，因为它拥有庞大的用户社区，遇到问题时很容易找到解决方案。Ubuntu的软件包管理系统apt使得安装软件变得简单，许多生物信息学软件都提供了Ubuntu版本的安装包或安装说明。

CentOS和RockyLinux属于企业级Linux发行版，它们源自Red Hat Enterprise Linux的源代码。这类发行版以稳定性为首要目标，软件版本相对保守但经过充分测试，适合在生产环境中使用。许多高性能计算集群和服务器都运行CentOS或类似的发行版。需要注意的是，CentOS 8之后，CentOS项目改变了发布策略，因此RockyLinux作为社区驱动的替代方案应运而生。如果你的工作环境使用的是这类企业级发行版，建议熟悉其软件包管理工具yum或dnf。

选择发行版时需要考虑实际工作环境。如果你所在的研究机构或实验室已经有现成的计算平台，那么最好使用与平台一致的发行版，这样可以避免环境差异带来的问题。如果你在自己的个人计算机上搭建分析环境，Ubuntu或其衍生版本如Linux Mint会是不错的选择。对于需要在服务器上部署的分析流程，CentOS或RockyLinux更为合适。无论选择哪个发行版，掌握基本的命令行操作和系统管理知识都是必要的。

文件系统结构
~~~~~~~~~~~~

Linux文件系统采用树状层次结构，所有文件和目录都从根目录开始。理解这个结构对于正确组织项目和查找文件至关重要。根目录用正斜杠表示，是整个文件系统的起点。在根目录下，有几个重要的系统目录，每个目录都有特定的用途。

home目录是普通用户的主目录，每个用户在home目录下都有一个以用户名命名的子目录。用户对自己的主目录拥有完全的读写执行权限，可以在这里创建文件、安装软件和组织项目。在生物信息学工作中，你的分析项目、数据文件和脚本通常都存放在主目录下的各个子目录中。主目录的路径可以用波浪号表示，例如当前用户的主目录可以简写为波浪号加斜杠，或者直接用波浪号。

usr目录存放系统范围的程序和数据。usr这个名字来源于Unix系统的早期，最初表示用户相关的程序，现在则表示Unix系统资源。usr目录下的bin子目录包含大多数用户级程序，lib目录包含程序运行所需的库文件，local目录用于存放本地安装的软件。当你从源代码编译安装生物信息学软件时，通常会将程序安装到usr/local目录下。

var目录存放系统运行过程中会频繁变化的文件，如日志文件、缓存文件和临时文件。对于生物信息学工作来说，var目录下的log子目录特别重要，系统服务和应用程序的日志文件存放在这里。当分析流程出现问题时，查看日志文件是排查问题的重要手段。tmp目录用于存放临时文件，系统重启后这个目录的内容会被清空。在处理大型数据文件时，一些程序会使用tmp目录作为临时存储空间，因此需要确保tmp目录有足够的可用空间。

etc目录存放系统配置文件。大多数系统服务和应用程序的配置文件都存放在etc目录或其子目录中。例如，环境变量的系统级配置、网络配置、用户信息等都存放在这里。了解etc目录的结构有助于你理解和修改系统配置。不过修改系统配置文件通常需要管理员权限，操作不当可能导致系统问题，因此需要谨慎行事。

proc和sys目录是虚拟文件系统，它们不对应实际的存储设备，而是内核提供的接口。通过读写这些目录下的文件，可以查看和修改内核参数、进程信息等系统状态。例如，proc目录下的meminfo文件包含内存使用信息，cpuinfo文件包含处理器信息。这些信息在性能监控和系统调优时非常有用。

用户与权限管理
~~~~~~~~~~~~~~

Linux是一个多用户操作系统，每个用户都有自己的账户和权限。这种设计保证了系统的安全性，不同用户之间的文件和进程相互隔离。理解用户和权限管理对于在共享计算环境中工作至关重要。

useradd命令用于创建新用户。创建用户时可以指定用户的主目录、默认shell、所属组等信息。例如，要创建一个名为bioinfo的用户，可以使用useradd命令配合适当的选项。创建用户后，通常需要使用passwd命令为用户设置密码。在实际工作中，你可能没有创建用户的权限，这需要系统管理员来执行。但是了解用户创建过程有助于理解系统的工作方式。

.. code-block:: bash

    # 创建新用户bioinfo
    sudo useradd -m -s /bin/bash bioinfo
    
    # 为用户设置密码
    sudo passwd bioinfo
    
    # 查看用户信息
    id bioinfo

每个文件和目录都有三组权限，分别对应文件所有者、文件所属组和其他用户。每组权限包含读、写、执行三种权限。读权限允许查看文件内容或列出目录内容，写权限允许修改文件或在目录中创建删除文件，执行权限允许运行程序或进入目录。使用ls命令的长格式选项可以查看文件权限，权限显示为一串十个字符，第一个字符表示文件类型，后九个字符分为三组，分别对应所有者、组和其他用户的权限。

.. code-block:: bash

    # 查看文件权限
    ls -l sample.fastq
    
    # 输出示例：
    # -rw-r--r-- 1 user group 1024 Jan 1 10:00 sample.fastq
    # 第一个字符-表示普通文件，d表示目录
    # 后九个字符分为三组：所有者权限、组权限、其他用户权限

chmod命令用于修改文件权限。权限可以用数字表示，读权限对应数字4，写权限对应数字2，执行权限对应数字1，三种权限相加得到一个0到7之间的数字。例如，chmod 755 file将文件权限设置为所有者可读写执行，组和其他用户可读执行。权限也可以用字母表示，使用u、g、o分别代表所有者、组和其他用户，使用加号和减号添加或移除权限。例如，chmod u+x file为所有者添加执行权限。

.. code-block:: bash

    # 使用数字设置权限
    chmod 755 script.sh
    
    # 使用字母添加权限
    chmod u+x script.sh
    
    # 为所有用户添加执行权限
    chmod a+x script.sh
    
    # 移除组的写权限
    chmod g-w data.txt

chown命令用于修改文件所有者和所属组。只有文件所有者或root用户才能修改文件的所有者。这个命令在协作项目中很有用，当你需要将文件的所有权转移给其他用户时使用。chgrp命令专门用于修改文件的所属组。在实际工作中，合理设置文件权限和所有者，既能保证数据安全，又能方便团队协作。

.. code-block:: bash

    # 修改文件所有者
    sudo chown bioinfo data.txt
    
    # 修改文件所有者和所属组
    sudo chown bioinfo:bioinfo data.txt
    
    # 递归修改目录及其内容的所有者
    sudo chown -R bioinfo:bioinfo /home/bioinfo/project

umask命令设置创建文件时的默认权限掩码。当创建新文件或目录时，系统会根据umask值自动设置权限。默认情况下，文件的初始权限是666，目录的初始权限是777，然后减去umask值得到最终权限。例如，umask值为022时，创建的文件权限为644，目录权限为755。可以在shell配置文件中设置umask值，使其在每次登录时自动生效。

.. code-block:: bash

    # 查看当前umask值
    umask
    
    # 设置umask值
    umask 022
    
    # 创建新文件，权限会是644
    touch newfile.txt
    ls -l newfile.txt

进程管理
~~~~~~~~

在Linux系统中，每个运行的程序都是一个进程。进程管理是系统管理的重要组成部分，了解如何查看和控制进程对于运行生物信息学分析流程至关重要。长时间运行的数据分析任务需要正确管理，避免意外中断或占用过多系统资源。

ps命令用于查看当前系统中的进程状态。不带任何选项的ps命令只显示当前终端启动的进程。使用aux选项可以显示系统中所有进程的详细信息，包括进程ID、CPU使用率、内存使用率、启动时间、运行时间等。进程ID是系统分配给每个进程的唯一标识符，在控制进程时需要使用这个ID。在生物信息学工作中，经常需要使用ps命令查看分析任务是否正在运行，以及资源使用情况。

.. code-block:: bash

    # 查看当前终端的进程
    ps
    
    # 查看所有进程的详细信息
    ps aux
    
    # 查找特定进程
    ps aux | grep python
    
    # 以用户格式显示进程
    ps -u username

top命令提供了一个动态的实时视图，显示系统中进程的资源使用情况。top命令默认每几秒刷新一次，按照CPU使用率排序显示进程。在top界面中，可以看到系统的总体资源使用情况，包括CPU使用率、内存使用率、交换空间使用率等。按下特定的键可以对进程进行排序、过滤或终止。top命令是监控系统状态和排查性能问题的有力工具。htop是top的增强版本，提供了更友好的界面和更多功能，如鼠标操作、颜色显示、进程树视图等。htop通常需要单独安装，但在许多系统上已经成为标准工具。

.. code-block:: bash

    # 启动top监控
    top
    
    # 在top界面中常用按键：
    # P - 按CPU使用率排序
    # M - 按内存使用率排序
    # k - 终止进程（需要输入PID）
    # q - 退出top
    
    # 安装htop
    sudo apt install htop
    
    # 启动htop
    htop

kill命令用于向进程发送信号，最常用的信号是终止信号。使用kill命令时需要指定进程ID。默认情况下，kill发送TERM信号，请求进程正常退出。如果进程没有响应，可以发送KILL信号强制终止。需要注意的是，强制终止进程可能导致数据丢失，应该作为最后的手段使用。除了终止进程，kill还可以发送其他信号，如让进程重新加载配置、暂停进程、继续运行暂停的进程等。

.. code-block:: bash

    # 正常终止进程
    kill 12345
    
    # 强制终止进程
    kill -9 12345
    
    # 查看可用的信号
    kill -l
    
    # 暂停进程
    kill -STOP 12345
    
    # 继续运行暂停的进程
    kill -CONT 12345

jobs命令显示当前shell会话中的后台任务。当你在命令后加上与符号时，命令会在后台运行，shell会立即返回提示符，你可以继续输入其他命令。jobs命令列出这些后台任务及其状态。使用fg命令可以将后台任务调到前台，使用bg命令可以在后台继续运行暂停的任务。这些命令在管理长时间运行的分析任务时非常有用。

.. code-block:: bash

    # 在后台运行命令
    python analysis.py &
    
    # 查看后台任务
    jobs
    
    # 将后台任务调到前台
    fg %1
    
    # 在后台继续运行暂停的任务
    bg %1

nohup命令让进程在用户退出登录后继续运行。默认情况下，当用户退出登录时，该用户启动的所有进程都会收到挂起信号并终止。使用nohup命令启动的进程会忽略挂起信号，即使你关闭终端或断开SSH连接，进程仍会继续运行。nohup通常与后台运行符号一起使用，并将输出重定向到文件。这对于运行需要数小时甚至数天的分析任务非常重要。

.. code-block:: bash

    # 使用nohup在后台运行任务
    nohup python analysis.py > output.log 2>&1 &
    
    # 查看任务是否在运行
    ps aux | grep analysis.py
    
    # 查看输出日志
    tail -f output.log

磁盘与文件管理
~~~~~~~~~~~~~~

磁盘空间管理和文件操作是日常工作中最频繁的任务。生物信息学数据文件通常很大，合理管理磁盘空间和熟练使用文件操作命令可以大大提高工作效率。

df命令显示文件系统的磁盘空间使用情况。使用h选项可以以人类可读的格式显示，将字节数转换为KB、MB、GB等单位。df命令显示每个文件系统的总容量、已用空间、可用空间和挂载点。在开始大型分析任务之前，应该检查目标文件系统的可用空间，确保有足够的空间存储中间文件和结果文件。如果磁盘空间不足，分析任务可能中途失败，浪费已经投入的计算时间。

.. code-block:: bash

    # 查看磁盘空间使用情况
    df -h
    
    # 查看特定文件系统的空间
    df -h /home
    
    # 查看inode使用情况
    df -i

du命令显示目录或文件的磁盘使用量。使用sh选项可以显示目录的总大小，而不是递归显示每个子目录的大小。du命令对于找出占用大量空间的文件和目录非常有用。当磁盘空间紧张时，可以使用du命令逐级查找，定位占用空间最大的文件或目录，然后决定是否删除或移动这些文件。结合sort命令，可以按大小排序显示目录。

.. code-block:: bash

    # 查看当前目录大小
    du -sh .
    
    # 查看各子目录大小
    du -h --max-depth=1
    
    # 查看当前目录下最大的10个文件或目录
    du -ah | sort -rh | head -n 10

ls命令是最基本的文件列表命令。使用l选项以长格式显示，可以看到文件的权限、所有者、大小、修改时间等详细信息。使用h选项以人类可读的格式显示文件大小。使用a选项显示包括隐藏文件在内的所有文件。ls命令支持多种排序方式，如按修改时间排序、按文件大小排序等。熟练使用ls命令的各种选项可以快速获取所需的文件信息。

.. code-block:: bash

    # 长格式显示文件信息
    ls -lh
    
    # 显示所有文件包括隐藏文件
    ls -la
    
    # 按修改时间排序
    ls -lt
    
    # 按文件大小排序
    ls -lhS
    
    # 递归列出子目录
    ls -R

cp命令用于复制文件和目录。复制文件时需要指定源文件和目标位置。使用r选项可以递归复制整个目录。使用p选项可以保留文件的属性，如修改时间、权限等。在复制大型文件时，可以显示进度信息。mv命令用于移动或重命名文件。移动文件时，文件从一个位置转移到另一个位置，原位置的文件消失。重命名文件实际上是移动文件到同一目录下的不同名称。rm命令用于删除文件和目录。使用r选项递归删除目录，使用f选项强制删除而不提示确认。删除操作不可恢复，使用时需要特别小心，尤其是在使用通配符时。

.. code-block:: bash

    # 复制文件
    cp source.txt destination.txt
    
    # 复制目录
    cp -r source_dir/ destination_dir/
    
    # 保留文件属性复制
    cp -p source.txt destination.txt
    
    # 移动文件
    mv oldname.txt newname.txt
    
    # 移动目录
    mv old_dir/ new_dir/
    
    # 删除文件
    rm file.txt
    
    # 递归删除目录
    rm -r directory/
    
    # 强制删除（谨慎使用）
    rm -rf directory/

ln命令用于创建链接。Linux支持两种类型的链接，硬链接和符号链接。硬链接是文件的另一个名称，指向相同的数据块。删除原文件后，硬链接仍然可以访问文件内容。符号链接类似于Windows的快捷方式，是一个指向另一个文件的特殊文件。如果原文件被删除，符号链接将失效。在生物信息学工作中，经常使用符号链接来组织文件，例如将分散在不同目录下的数据文件链接到统一的工作目录，而不需要复制文件本身。

.. code-block:: bash

    # 创建硬链接
    ln source.txt hardlink.txt
    
    # 创建符号链接
    ln -s source.txt symlink.txt
    
    # 创建目录的符号链接
    ln -s /path/to/data data_link
    
    # 查看链接信息
    ls -l symlink.txt

文本编辑器
~~~~~~~~~~

在Linux环境中工作，熟练使用命令行文本编辑器是必不可少的技能。虽然图形界面的编辑器更加友好，但在远程服务器上工作时，命令行编辑器往往是唯一的选择。vim和nano是两个最常用的命令行编辑器，各有特点。

vim是vi编辑器的增强版本，是Linux系统中最强大的文本编辑器之一。vim采用模式编辑的设计，有不同的工作模式，如普通模式、插入模式、命令模式等。在普通模式下，按键被解释为命令而不是文本输入，这使得编辑操作非常高效。按i键进入插入模式，可以正常输入文本。按Esc键返回普通模式。在普通模式下，可以使用h、j、k、l键移动光标，使用dd删除整行，使用yy复制整行，使用p粘贴。在命令模式下，可以执行保存、退出、搜索替换等操作。按冒号进入命令模式，输入w保存文件，输入q退出，输入wq保存并退出，输入q感叹号强制退出不保存。

vim的学习曲线较陡峭，但掌握后编辑效率极高。vim支持语法高亮、代码补全、多文件编辑等高级功能。通过配置文件可以定制vim的行为，安装插件可以扩展vim的功能。对于需要频繁编辑文本文件的用户，投入时间学习vim是值得的。vim还提供了丰富的在线帮助文档，在vim中输入冒号help可以查看帮助。

nano是一个简单易用的文本编辑器，适合初学者使用。nano没有复杂的模式概念，打开文件后就可以直接编辑。屏幕底部显示常用的快捷键提示，使用Control键配合字母键执行操作。例如，Control加O保存文件，Control加X退出编辑器，Control加W搜索文本。nano的功能相对简单，但对于快速编辑配置文件或脚本已经足够。许多Linux发行版将nano设置为默认编辑器。

选择哪个编辑器取决于个人偏好和工作需求。如果只是偶尔编辑简单的配置文件，nano就足够了。如果需要长时间编写代码或处理大量文本，vim的高效编辑能力会带来明显的优势。建议初学者先掌握nano的基本使用，然后逐步学习vim。在实际工作中，两个编辑器都了解会更有灵活性。

环境变量
~~~~~~~~

环境变量是操作系统中用于存储配置信息的动态值。程序在运行时可以读取环境变量的值来获取配置信息。理解环境变量的概念和使用方法对于配置分析环境和解决软件依赖问题非常重要。

PATH是最重要的环境变量之一，它定义了系统查找可执行文件的目录列表。当你在命令行输入一个命令时，系统会在PATH变量指定的目录中依次查找对应的可执行文件。PATH变量的值是一系列目录路径，用冒号分隔。如果某个程序的可执行文件不在PATH指定的目录中，你需要输入完整路径才能运行该程序，或者将该程序所在目录添加到PATH中。在生物信息学工作中，经常需要将自行安装的软件目录添加到PATH，以便在任何位置都能直接运行这些软件。

LD_LIBRARY_PATH变量指定动态链接库的搜索路径。许多程序在运行时需要加载共享库文件，系统会在标准库目录和LD_LIBRARY_PATH指定的目录中查找这些库文件。如果程序依赖的库文件不在标准位置，需要通过设置LD_LIBRARY_PATH来告诉系统库文件的位置。不过，过度依赖LD_LIBRARY_PATH可能导致库版本冲突，更好的做法是在编译时指定库路径或使用ldconfig配置系统库缓存。

PYTHONPATH变量指定Python模块的搜索路径。当Python程序导入模块时，解释器会在PYTHONPATH指定的目录中查找模块文件。如果你开发了自定义的Python模块或安装了不在标准位置的第三方库，可以通过设置PYTHONPATH让Python能够找到这些模块。

查看环境变量的值可以使用echo命令配合美元符号。例如，echo $PATH显示PATH变量的值。设置环境变量使用export命令。例如，export PATH=/new/path:$PATH将新路径添加到PATH变量的开头。这种设置只在当前shell会话中有效，退出登录后会失效。要使环境变量设置永久生效，需要将export命令添加到shell配置文件中。对于bash shell，通常编辑主目录下的bashrc文件。

.. code-block:: bash

    # 查看PATH环境变量
    echo $PATH
    
    # 查看所有环境变量
    env
    
    # 设置环境变量（临时）
    export PATH=/home/user/bin:$PATH
    
    # 设置PYTHONPATH
    export PYTHONPATH=/home/user/python_modules:$PYTHONPATH
    
    # 设置LD_LIBRARY_PATH
    export LD_LIBRARY_PATH=/home/user/lib:$LD_LIBRARY_PATH
    
    # 在bashrc中添加永久设置
    echo 'export PATH=/home/user/bin:$PATH' >> ~/.bashrc
    source ~/.bashrc

压缩与归档
~~~~~~~~~~

生物信息学数据文件通常很大，使用压缩可以节省存储空间并加快文件传输速度。Linux提供了多种压缩和归档工具，了解它们的特点和使用方法可以更有效地管理数据文件。

tar命令是最常用的归档工具，它可以将多个文件和目录打包成一个归档文件。tar本身只进行打包，不进行压缩，但可以与其他压缩工具配合使用。使用c选项创建归档，x选项解包，v选项显示处理过程，f选项指定归档文件名。例如，tar cvf archive.tar directory将directory目录打包成archive.tar文件。使用z选项可以在打包的同时使用gzip压缩，生成tar.gz或tgz文件。使用j选项使用bzip2压缩，生成tar.bz2文件。tar命令保留文件的权限和时间戳信息，适合用于备份和传输完整的目录结构。

.. code-block:: bash

    # 创建tar归档文件
    tar cvf archive.tar directory/
    
    # 创建gzip压缩的tar文件
    tar czvf archive.tar.gz directory/
    
    # 创建bzip2压缩的tar文件
    tar cjvf archive.tar.bz2 directory/
    
    # 解压tar文件
    tar xvf archive.tar
    
    # 解压tar.gz文件
    tar xzvf archive.tar.gz
    
    # 解压tar.bz2文件
    tar xjvf archive.tar.bz2
    
    # 解压到指定目录
    tar xzvf archive.tar.gz -C /target/directory/

gzip是最常用的压缩工具，使用DEFLATE算法进行压缩。gzip只能压缩单个文件，不能压缩目录。使用gzip filename压缩文件，生成filename.gz文件，原文件被删除。使用gunzip或gzip -d解压文件。gzip提供了不同压缩级别选项，级别1压缩速度最快但压缩率较低，级别9压缩率最高但速度最慢，默认使用级别6。对于大型数据文件，可以选择较低的压缩级别以加快处理速度，或者选择较高的压缩级别以节省存储空间。

.. code-block:: bash

    # 压缩文件
    gzip file.txt
    
    # 解压文件
    gunzip file.txt.gz
    
    # 使用最快压缩速度
    gzip -1 file.txt
    
    # 使用最高压缩率
    gzip -9 file.txt
    
    # 保留原文件
    gzip -c file.txt > file.txt.gz

bzip2使用Burrows-Wheeler算法进行压缩，通常比gzip有更高的压缩率，但压缩和解压速度较慢。bzip2适合用于需要长期存储且不频繁访问的文件。zip和unzip命令处理与Windows兼容的zip格式文件。zip格式在跨平台文件交换中很常用。与gzip不同，zip可以将多个文件压缩到一个zip文件中。

.. code-block:: bash

    # 使用bzip2压缩
    bzip2 file.txt
    
    # 解压bzip2文件
    bunzip2 file.txt.bz2
    
    # 创建zip文件
    zip archive.zip file1.txt file2.txt
    
    # 压缩目录
    zip -r archive.zip directory/
    
    # 解压zip文件
    unzip archive.zip
    
    # 解压到指定目录
    unzip archive.zip -d /target/directory/

在实际工作中，选择压缩工具需要考虑压缩率、速度和兼容性。对于日常使用，tar配合gzip是最常见的选择。如果需要更高的压缩率且不介意较慢的速度，可以使用bzip2或xz。如果需要与Windows用户交换文件，使用zip格式。

远程访问
~~~~~~~~

生物信息学分析通常在远程服务器或计算集群上进行，掌握远程访问和文件传输的方法是必备技能。SSH协议提供了安全的远程登录和文件传输功能。

ssh命令用于建立安全的远程连接。使用ssh username@hostname连接到远程主机，其中username是远程主机上的用户名，hostname是远程主机的地址。首次连接时，ssh会显示远程主机的指纹信息并询问是否继续连接，确认后指纹信息会被保存，后续连接不再提示。连接建立后，你可以在远程主机上执行命令，就像在本地终端一样。使用exit命令断开连接。ssh支持密钥认证，使用ssh-keygen生成密钥对，将公钥复制到远程主机的authorized_keys文件中，之后连接时不需要输入密码。密钥认证不仅方便，而且比密码认证更安全。

.. code-block:: bash

    # 连接到远程主机
    ssh username@hostname
    
    # 使用特定端口连接
    ssh -p 2222 username@hostname
    
    # 在远程主机上执行单个命令
    ssh username@hostname ls -l
    
    # 生成SSH密钥对
    ssh-keygen -t rsa -b 4096
    
    # 将公钥复制到远程主机
    ssh-copy-id username@hostname
    
    # 使用密钥文件连接
    ssh -i ~/.ssh/mykey username@hostname

scp命令用于在本地和远程主机之间复制文件。scp使用与ssh相同的认证方式，提供安全的文件传输。复制文件到远程主机使用scp localfile username@hostname:/remote/path。从远程主机复制文件使用scp username@hostname:/remote/file /local/path。使用r选项可以递归复制整个目录。scp命令简单易用，适合传输少量文件。

.. code-block:: bash

    # 复制文件到远程主机
    scp localfile.txt username@hostname:/remote/path/
    
    # 从远程主机复制文件
    scp username@hostname:/remote/file.txt /local/path/
    
    # 递归复制目录
    scp -r localdir/ username@hostname:/remote/path/
    
    # 使用特定端口
    scp -P 2222 file.txt username@hostname:/remote/path/

rsync是一个更强大的文件同步工具，它只传输文件的变化部分，支持断点续传，适合传输大量文件或大文件。rsync的选项非常丰富，常用的组合是avz选项，a选项表示归档模式，保留文件属性，v选项显示详细信息，z选项启用压缩传输。rsync source destination格式与cp命令类似，但source和destination可以是本地路径或远程路径。远程路径使用username@hostname:path格式。使用rsync进行备份或同步时，可以删除目标位置多余的文件，使目标与源完全一致，这需要使用delete选项。

.. code-block:: bash

    # 同步本地目录到远程主机
    rsync -avz localdir/ username@hostname:/remote/path/
    
    # 从远程主机同步到本地
    rsync -avz username@hostname:/remote/path/ localdir/
    
    # 显示传输进度
    rsync -avz --progress localdir/ username@hostname:/remote/path/
    
    # 删除目标位置多余的文件
    rsync -avz --delete localdir/ username@hostname:/remote/path/
    
    # 排除特定文件
    rsync -avz --exclude '*.log' localdir/ username@hostname:/remote/path/

sftp是一个交互式的安全文件传输程序，提供类似FTP的命令界面。使用sftp username@hostname连接后，可以使用get命令下载文件，put命令上传文件，ls命令列出远程文件，lls命令列出本地文件。sftp适合需要浏览远程文件系统并选择性传输文件的场景。

.. code-block:: bash

    # 连接到sftp
    sftp username@hostname
    
    # 在sftp交互界面中常用命令：
    # get remote_file - 下载文件
    # put local_file - 上传文件
    # ls - 列出远程文件
    # lls - 列出本地文件
    # cd - 切换远程目录
    # lcd - 切换本地目录
    # mkdir - 创建远程目录
    # exit - 退出sftp

管道与重定向
~~~~~~~~~~~~

管道和重定向是Linux命令行的强大特性，它们允许你将一个命令的输出作为另一个命令的输入，或者将输出保存到文件。掌握这些技术可以将多个简单命令组合成复杂的数据处理流程。

重定向操作符将命令的输入或输出指向文件。大于号将命令的标准输出重定向到文件，如果文件已存在则覆盖。两个大于号将输出追加到文件末尾而不覆盖原有内容。小于号将文件内容作为命令的标准输入。例如，ls > filelist.txt将当前目录的文件列表保存到filelist.txt文件。echo "New log entry" >> logfile.txt向日志文件追加一条记录。在生物信息学工作中，经常需要将分析结果保存到文件，或者从文件读取输入数据。

.. code-block:: bash

    # 将输出重定向到文件
    ls > filelist.txt
    
    # 追加输出到文件
    echo "New log entry" >> logfile.txt
    
    # 从文件读取输入
    wc -l < file.txt
    
    # 同时重定向标准输出和标准错误
    command > output.log 2>&1
    
    # 或者使用更简洁的写法
    command &> output.log
    
    # 丢弃输出
    command > /dev/null 2>&1

标准错误输出是另一个输出流，用于输出错误信息和诊断信息。使用2>将标准错误重定向到文件。例如，command 2> error.log将命令的错误信息保存到error.log文件。使用2>&1将标准错误重定向到标准输出，这样可以将所有输出保存到同一个文件。command > output.log 2>&1将标准输出和标准错误都保存到output.log文件。另一种写法是command &> output.log，效果相同。

管道符将一个命令的标准输出连接到另一个命令的标准输入，使数据在命令之间流动。通过管道，可以将多个命令串联起来，每个命令处理前一个命令的输出。例如，cat file.txt | grep "pattern" | wc -l计算文件中包含特定模式的行数。cat命令输出文件内容，grep命令过滤包含模式的行，wc命令统计行数。管道使得复杂的数据处理任务可以分解为一系列简单的步骤，每个步骤使用专门的工具完成。

.. code-block:: bash

    # 基本管道使用
    cat file.txt | grep "pattern"
    
    # 多个管道串联
    cat file.txt | grep "pattern" | wc -l
    
    # 提取特定列并排序
    cut -f1 data.txt | sort | uniq
    
    # 查找并处理文件
    find . -name "*.fastq" | xargs wc -l
    
    # 组合重定向和管道
    cat file.txt | grep "pattern" > results.txt

管道和重定向可以灵活组合。例如，find . -name "*.txt" | xargs grep "keyword" > results.txt查找当前目录下所有txt文件中包含关键词的行，并将结果保存到文件。这个命令组合了文件查找、内容搜索和结果保存三个步骤。在构建生物信息学分析流程时，管道和重定向是连接各个处理步骤的基本手段。

查找与过滤
~~~~~~~~~~

在处理大量文件和数据时，快速查找文件和过滤内容是常见需求。Linux提供了多种查找和过滤工具，各有特点和适用场景。

find命令是功能最强大的文件查找工具，它可以根据文件名、类型、大小、修改时间等多种条件查找文件。find命令的基本格式是find path expression，path指定查找的起始目录，expression指定查找条件。例如，find . -name "*.fastq"查找当前目录及其子目录下所有fastq文件。find /data -type f -size +100M查找data目录下大于100MB的文件。find . -mtime -7查找最近7天内修改过的文件。find命令还可以对找到的文件执行操作，例如find . -name "*.log" -delete删除所有log文件，find . -name "*.txt" -exec wc -l {} \;统计每个txt文件的行数。

.. code-block:: bash

    # 按文件名查找
    find . -name "*.fastq"
    
    # 按文件类型查找
    find . -type f  # 查找文件
    find . -type d  # 查找目录
    
    # 按文件大小查找
    find . -size +100M  # 大于100MB
    find . -size -1M    # 小于1MB
    
    # 按修改时间查找
    find . -mtime -7   # 最近7天内修改
    find . -mtime +30  # 超过30天前修改
    
    # 对找到的文件执行操作
    find . -name "*.log" -delete
    find . -name "*.txt" -exec wc -l {} \;

locate命令通过搜索预先建立的文件数据库快速查找文件。与find不同，locate不需要实时遍历文件系统，因此速度非常快。但数据库不是实时更新的，新创建的文件可能找不到。updatedb命令更新数据库，通常由系统定期自动执行。locate命令的使用很简单，locate filename查找文件名包含指定字符串的所有文件。在需要快速查找文件且对实时性要求不高时，locate是很好的选择。

.. code-block:: bash

    # 查找文件
    locate sample.fastq
    
    # 更新数据库
    sudo updatedb
    
    # 统计匹配文件数量
    locate -c pattern
    
    # 限制输出数量
    locate -n 10 pattern

grep命令是最常用的文本搜索工具，它在文件中搜索匹配指定模式的行。grep支持基本正则表达式和扩展正则表达式。例如，grep "gene" file.txt查找文件中包含gene的行。使用i选项忽略大小写，v选项显示不匹配的行，c选项只显示匹配行数，n选项显示行号。使用r选项递归搜索目录下的所有文件。grep "error" /var/log递归搜索日志目录中的错误信息。在生物信息学中，grep常用于从序列文件或注释文件中提取特定信息。

.. code-block:: bash

    # 基本搜索
    grep "gene" file.txt
    
    # 忽略大小写
    grep -i "GENE" file.txt
    
    # 显示行号
    grep -n "gene" file.txt
    
    # 显示不匹配的行
    grep -v "gene" file.txt
    
    # 只显示匹配行数
    grep -c "gene" file.txt
    
    # 递归搜索目录
    grep -r "error" /var/log/
    
    # 使用扩展正则表达式
    grep -E "gene|protein" file.txt

ag和rg是现代化的文本搜索工具，它们比grep更快，默认支持递归搜索，自动忽略版本控制目录和二进制文件。ag是The Silver Searcher，rg是ripgrep。这两个工具在大型代码库或数据目录中搜索时优势明显。它们的使用方式与grep类似，但提供了更好的默认行为和更快的搜索速度。如果你的系统上安装了这些工具，建议在日常工作中使用它们替代grep。

.. code-block:: bash

    # 使用ag搜索
    ag "pattern" directory/
    
    # 使用rg搜索
    rg "pattern" directory/
    
    # 只显示文件名
    ag -l "pattern"
    
    # 显示上下文行
    ag -C 3 "pattern"  # 显示匹配行前后各3行

计划任务
~~~~~~~~

计划任务功能允许你在指定的时间自动执行命令或脚本，这对于定期运行分析任务、执行系统维护非常有用。cron是Linux系统中最常用的计划任务工具。

cron守护进程在后台运行，每分钟检查一次是否有任务需要执行。用户的cron任务配置称为crontab，使用crontab命令管理。crontab -e编辑当前用户的cron配置，crontab -l列出当前用户的cron任务，crontab -r删除所有cron任务。cron配置文件的每一行定义一个任务，格式为五个时间字段后跟要执行的命令。五个时间字段分别表示分钟、小时、日、月、星期，使用星号表示任意值，使用逗号分隔多个值，使用短横线表示范围，使用斜杠表示间隔。

.. code-block:: bash

    # 编辑cron任务
    crontab -e
    
    # 列出当前用户的cron任务
    crontab -l
    
    # 删除所有cron任务
    crontab -r
    
    # cron任务示例（在crontab文件中添加）
    # 每天凌晨2点执行备份脚本
    0 2 * * * /home/user/backup.sh
    
    # 每30分钟执行一次检查脚本
    */30 * * * * /home/user/check.sh
    
    # 周一到周五每天上午9点执行工作脚本
    0 9 * * 1-5 /home/user/work.sh
    
    # 将输出重定向到日志文件
    0 2 * * * /home/user/backup.sh > /home/user/backup.log 2>&1

例如，0 2 * * * /home/user/backup.sh表示每天凌晨2点执行backup.sh脚本。*/30 * * * * /home/user/check.sh表示每30分钟执行一次check.sh脚本。0 9 * * 1-5 /home/user/work.sh表示周一到周五每天上午9点执行work.sh脚本。cron任务的输出会通过邮件发送给用户，可以重定向到文件或丢弃。例如，0 2 * * * /home/user/backup.sh > /home/user/backup.log 2>&1将输出和错误信息保存到日志文件。

at命令用于一次性计划任务，在指定的时间执行一次命令。使用at time指定执行时间，然后输入要执行的命令，按Ctrl+D结束输入。例如，at 10pm tomorrow在明天晚上10点执行。atq命令列出待执行的任务，atrm命令删除任务。at适合需要在特定时间执行一次的任务，而cron适合周期性执行的任务。

.. code-block:: bash

    # 在指定时间执行命令
    at 10pm tomorrow
    
    # 在at提示符下输入命令
    at> python analysis.py
    at> <EOT>  # 按Ctrl+D结束
    
    # 从文件读取任务
    at -f script.sh 10pm tomorrow
    
    # 列出待执行的at任务
    atq
    
    # 删除at任务
    atrm job_number
    
    # 在特定日期时间执行
    at 14:30 December 31

在使用计划任务时需要注意环境变量问题。cron执行任务时的环境变量与登录shell不同，PATH变量可能不包含你需要的目录。因此在cron任务中最好使用完整路径，或者在脚本开头设置必要的环境变量。另外，cron任务的标准输出和错误输出默认会通过邮件系统发送，如果邮件系统未配置，输出可能丢失。建议在cron任务中明确重定向输出到文件。

Shell 脚本编程
--------------

通配符与glob模式
~~~~~~~~~~~~~~~~

Shell脚本中最基本的功能之一是使用通配符来匹配文件名，这种模式称为glob模式。通配符可以让你用简单的模式匹配多个文件，而不需要逐个指定文件名。这在处理批量文件时非常有用，例如处理一组测序数据文件或日志文件。

星号是最常用的通配符，它匹配任意长度的任意字符。例如，*.fastq匹配当前目录下所有以fastq结尾的文件，sample_*匹配所有以sample_开头的文件。星号可以出现在模式的任何位置，也可以使用多个星号。例如，*/*.txt匹配所有子目录下的txt文件，sample*_R1*.fastq匹配文件名中包含sample和R1的fastq文件。需要注意的是，星号不匹配以点开头的隐藏文件，除非明确指定点字符。

问号匹配单个任意字符。例如，sample?.fastq匹配sample1.fastq、sampleA.fastq等，但不匹配sample10.fastq，因为问号只匹配一个字符。问号在需要匹配特定长度的文件名时很有用。

方括号用于匹配方括号内的任意一个字符。例如，sample[123].fastq匹配sample1.fastq、sample2.fastq和sample3.fastq，但不匹配sample4.fastq。可以在方括号中使用范围表示法，例如[a-z]匹配任意小写字母，[0-9]匹配任意数字。在方括号开头使用感叹号或脱字符表示不匹配方括号内的字符，例如sample[!0-9].fastq匹配sample后跟非数字字符的文件。

大括号扩展是另一种强大的模式匹配机制，它可以生成多个字符串。例如，{sample1,sample2,sample3}.fastq展开为sample1.fastq sample2.fastq sample3.fastq。大括号支持范围表示法，例如{1..10}展开为1到10的数字，{a..z}展开为所有小写字母。大括号可以嵌套使用，例如sample{A,B}{1,2}.fastq展开为sampleA1.fastq sampleA2.fastq sampleB1.fastq sampleB2.fastq。

.. code-block:: bash

    # 星号通配符示例
    ls *.fastq          # 列出所有fastq文件
    ls sample_*         # 列出所有以sample_开头的文件
    ls */*.txt          # 列出所有子目录下的txt文件
    
    # 问号通配符示例
    ls sample?.fastq    # 匹配sample1.fastq, sampleA.fastq等
    ls file?.txt        # 匹配file1.txt, fileA.txt等
    
    # 方括号通配符示例
    ls sample[123].fastq    # 匹配sample1, sample2, sample3
    ls file[a-z].txt        # 匹配filea.txt到filez.txt
    ls sample[!0-9].fastq   # 匹配sample后非数字的文件
    
    # 大括号扩展示例
    ls {sample1,sample2,sample3}.fastq
    ls sample{1..10}.fastq
    mkdir project_{A,B,C}_{1,2,3}
    
    # 组合使用
    ls {sample,control}_*_{1,2}.fastq

变量与数组
~~~~~~~~~~

Shell脚本中的变量用于存储数据，变量名以字母或下划线开头，可以包含字母、数字和下划线。变量赋值时等号两边不能有空格，使用变量时需要在变量名前加美元符号。Shell变量默认是字符串类型，但可以进行整数运算。

定义变量很简单，variable=value的格式即可。例如，data_dir=/data/project定义了data_dir变量。使用变量时，$variable或${variable}都可以，后者在变量名后紧跟其他字符时更清晰。例如，echo "Data directory is $data_dir"输出变量的值。变量可以在脚本中随时修改，新的赋值会覆盖旧的值。

.. code-block:: bash

    # 定义变量
    data_dir=/data/project
    sample_name="sample_001"
    
    # 使用变量
    echo $data_dir
    echo "Sample name is $sample_name"
    echo "${sample_name}_R1.fastq"  # 使用大括号避免歧义
    
    # 变量运算
    count=10
    count=$((count + 1))
    echo $count  # 输出11
    
    # 命令替换
    current_date=$(date +%Y%m%d)
    file_count=`ls | wc -l`  # 旧式写法

数组可以存储多个值，bash支持一维数组。数组定义使用圆括号，元素之间用空格分隔。例如，samples=(sample1 sample2 sample3)定义了一个包含三个元素的数组。访问数组元素使用${array[index]}的格式，索引从0开始。${samples[0]}访问第一个元素，${samples[@]}或${samples[*]}访问所有元素。${#samples[@]}获取数组长度。数组可以追加元素，samples+=(sample4)向数组末尾添加新元素。

.. code-block:: bash

    # 定义数组
    samples=(sample1 sample2 sample3)
    
    # 访问数组元素
    echo ${samples[0]}     # 输出sample1
    echo ${samples[@]}     # 输出所有元素
    
    # 获取数组长度
    echo ${#samples[@]}    # 输出3
    
    # 遍历数组
    for sample in ${samples[@]}; do
        echo "Processing $sample"
    done
    
    # 追加元素
    samples+=(sample4 sample5)
    
    # 通过索引赋值
    samples[10]=sample10

环境变量是特殊的变量，它们在所有子进程中都可用。使用export命令将变量导出为环境变量。例如，export PATH=/new/bin:$PATH将新路径添加到PATH环境变量。环境变量通常在shell配置文件中设置，以便每次登录时自动生效。常见的环境变量包括PATH、HOME、USER、PWD等。

.. code-block:: bash

    # 导出环境变量
    export MY_VAR="hello"
    export PATH=/home/user/bin:$PATH
    
    # 查看环境变量
    echo $PATH
    echo $HOME
    echo $USER
    
    # 在脚本中使用环境变量
    python_script_dir=/usr/local/bin
    export PATH=$python_script_dir:$PATH

特殊变量由shell自动设置，它们提供了有用的信息。$?保存上一个命令的退出状态，0表示成功，非0表示失败。$$是当前进程的进程ID。$!是最近一个后台进程的进程ID。$#是脚本参数的个数。$@和$*是所有的位置参数。

.. code-block:: bash

    # 特殊变量示例
    echo "Script name: $0"
    echo "First argument: $1"
    echo "Number of arguments: $#"
    echo "All arguments: $@"
    
    # 检查上一个命令的退出状态
    ls /nonexistent
    if [ $? -ne 0 ]; then
        echo "Command failed"
    fi
    
    # 使用进程ID创建临时文件
    temp_file="/tmp/script_$$"
    echo "Process ID: $$"

条件判断
~~~~~~~~

条件判断是控制脚本执行流程的重要手段。Shell提供了if语句和case语句两种条件判断结构，以及test命令和方括号用于条件测试。

if语句的基本结构是if condition; then commands; fi。condition是一个命令或测试表达式，如果命令的退出状态为0，则执行then后面的命令。可以使用else子句指定条件不满足时执行的命令，使用elif子句添加更多条件分支。例如，检查文件是否存在，如果存在则处理，否则输出错误信息。

.. code-block:: bash

    # 基本if语句
    if [ -f "data.txt" ]; then
        echo "File exists"
    fi
    
    # if-else语句
    if [ -f "data.txt" ]; then
        echo "File exists"
    else
        echo "File not found"
    fi
    
    # if-elif-else语句
    if [ $count -gt 100 ]; then
        echo "Large dataset"
    elif [ $count -gt 10 ]; then
        echo "Medium dataset"
    else
        echo "Small dataset"
    fi

test命令用于测试文件属性、字符串和数值。test -f file检查文件是否存在且是普通文件，test -d directory检查目录是否存在，test -r file检查文件是否可读，test -w file检查文件是否可写，test -x file检查文件是否可执行。test命令可以用方括号替代，[ condition ]与test condition等价。注意方括号内两侧必须有空格。

.. code-block:: bash

    # 文件测试
    if [ -f "file.txt" ]; then
        echo "文件存在且是普通文件"
    fi
    
    if [ -d "directory" ]; then
        echo "目录存在"
    fi
    
    if [ -r "file.txt" ]; then
        echo "文件可读"
    fi
    
    # 字符串测试
    if [ -z "$var" ]; then
        echo "字符串为空"
    fi
    
    if [ "$str1" = "$str2" ]; then
        echo "字符串相等"
    fi
    
    # 数值比较
    if [ $num -eq 10 ]; then
        echo "数值等于10"
    fi
    
    if [ $num -gt 100 ]; then
        echo "数值大于100"
    fi

字符串测试包括检查字符串是否为空、是否相等。test -z string检查字符串长度是否为0，test -n string检查字符串长度是否非0，test string1 = string2检查两个字符串是否相等，test string1 != string2检查两个字符串是否不相等。数值比较使用-eq、-ne、-lt、-le、-gt、-ge等操作符，分别表示等于、不等于、小于、小于等于、大于、大于等于。

case语句适合多分支条件判断，特别是当一个变量可能有多个离散值时。case语句的结构是case variable in pattern1) commands ;; pattern2) commands ;; esac。每个模式以右括号结束，命令以双分号结束。模式可以使用通配符，例如*.txt匹配所有txt文件。可以使用竖线分隔多个模式，例如yes|y|Y)匹配yes、y或Y。

.. code-block:: bash

    # case语句示例
    case $file_type in
        *.fastq)
            echo "Processing FASTQ file"
            ;;
        *.fasta)
            echo "Processing FASTA file"
            ;;
        *.bam)
            echo "Processing BAM file"
            ;;
        *)
            echo "Unknown file type"
            ;;
    esac
    
    # 多模式匹配
    case $answer in
        yes|y|Y)
            echo "Confirmed"
            ;;
        no|n|N)
            echo "Cancelled"
            ;;
        *)
            echo "Invalid input"
            ;;
    esac

循环
~~~~

循环结构允许重复执行一组命令，Shell提供了for、while和until三种循环。循环在处理批量文件、迭代计算等场景中非常有用。

for循环遍历一个列表中的每个元素，对每个元素执行相同的命令。基本格式是for variable in list; do commands; done。list可以是空格分隔的值列表、命令的输出、或文件通配符匹配的结果。例如，for file in *.fastq; do echo "Processing $file"; done遍历当前目录下所有fastq文件。for循环也可以使用C风格的语法，for ((i=0; i<10; i++)); do commands; done，这在需要数值迭代时更方便。

.. code-block:: bash

    # 遍历文件列表
    for file in *.fastq; do
        echo "Processing $file"
        wc -l "$file"
    done
    
    # 遍历自定义列表
    for sample in sample1 sample2 sample3; do
        echo "Analyzing $sample"
    done
    
    # 遍历命令输出
    for dir in $(ls -d */); do
        echo "Directory: $dir"
    done
    
    # C风格for循环
    for ((i=1; i<=10; i++)); do
        echo "Number: $i"
    done
    
    # 遍历数组
    samples=(sample1 sample2 sample3)
    for sample in ${samples[@]}; do
        echo "Processing $sample"
    done

while循环在条件为真时重复执行命令。基本格式是while condition; do commands; done。condition是一个命令或测试表达式，如果退出状态为0，则继续执行循环体。例如，逐行读取文件内容，while read line; do echo "$line"; done < file.txt。while循环常用于需要持续执行直到某个条件满足的场景，例如等待某个文件出现或某个进程结束。

.. code-block:: bash

    # 逐行读取文件
    while read line; do
        echo "Line: $line"
    done < data.txt
    
    # 计数器循环
    count=0
    while [ $count -lt 10 ]; do
        echo "Count: $count"
        count=$((count + 1))
    done
    
    # 等待文件出现
    while [ ! -f "output.txt" ]; do
        echo "Waiting for output..."
        sleep 5
    done
    echo "File is ready"
    
    # 无限循环（需要手动退出）
    while true; do
        echo "Running..."
        sleep 1
    done

until循环与while循环相反，它在条件为假时重复执行命令。until condition; do commands; done等价于while ! condition; do commands; done。until循环在需要执行命令直到某个条件成立时使用。

.. code-block:: bash

    # until循环示例
    count=0
    until [ $count -ge 10 ]; do
        echo "Count: $count"
        count=$((count + 1))
    done
    
    # 等待进程结束
    until ! ps -p $pid > /dev/null; do
        echo "Process still running"
        sleep 5
    done
    echo "Process finished"

break命令用于跳出循环，continue命令用于跳过当前迭代进入下一次迭代。break n可以跳出n层循环。在嵌套循环中，这些命令提供了灵活的流程控制。

.. code-block:: bash

    # 使用break跳出循环
    for i in {1..10}; do
        if [ $i -eq 5 ]; then
            break
        fi
        echo "Number: $i"
    done
    
    # 使用continue跳过迭代
    for i in {1..10}; do
        if [ $i -eq 5 ]; then
            continue
        fi
        echo "Number: $i"
    done
    
    # 嵌套循环中使用break
    for i in {1..5}; do
        for j in {1..5}; do
            if [ $j -eq 3 ]; then
                break 2  # 跳出两层循环
            fi
            echo "i=$i, j=$j"
        done
    done

函数定义与调用
~~~~~~~~~~~~~~

函数是将一组命令组织在一起以便重复使用的机制。使用函数可以使脚本更加模块化和易于维护。函数在生物信息学脚本中特别有用，可以将常用的数据处理步骤封装成函数。

函数定义有两种格式。一种是使用function关键字，function name { commands; }。另一种是使用圆括号，name() { commands; }。两种格式功能相同，后者更通用。函数定义后，可以像使用命令一样调用函数，直接使用函数名即可。

函数可以接受参数，在函数体内使用$1、$2等位置参数访问参数值。$#表示参数个数，$@表示所有参数。这些位置参数是函数的局部参数，不会影响脚本的全局位置参数。例如，定义一个处理文件的函数，process_file() { local file=$1; echo "Processing $file"; }，然后调用process_file sample.fastq。

函数可以使用local关键字定义局部变量，局部变量只在函数内部可见，不会影响全局变量。这是避免变量名冲突的重要手段。函数可以使用return命令返回一个退出状态，return value返回指定的值，如果不指定则返回最后一个命令的退出状态。

函数可以返回字符串值，但需要使用命令替换来捕获输出。例如，get_file_count() { echo $(ls -1 | wc -l); }，然后count=$(get_file_count)获取返回值。另一种方法是使用全局变量来传递返回值。

.. code-block:: bash

    # 函数定义的两种格式
    # 格式一：使用function关键字
    function greet {
        echo "Hello, $1!"
    }
    
    # 格式二：使用圆括号（更通用）
    greet() {
        echo "Hello, $1!"
    }
    
    # 调用函数
    greet "World"
    greet "Bioinformatics"
    
    # 带参数的函数
    process_file() {
        local file=$1
        local output_dir=$2
        
        echo "Processing file: $file"
        echo "Output directory: $output_dir"
        
        if [ -f "$file" ]; then
            wc -l "$file"
            return 0
        else
            echo "Error: File not found"
            return 1
        fi
    }
    
    # 调用带参数的函数
    process_file "data.fastq" "results/"
    
    # 函数返回字符串值
    get_file_count() {
        local dir=$1
        echo $(ls -1 "$dir" 2>/dev/null | wc -l)
    }
    
    # 捕获函数返回值
    count=$(get_file_count "data/")
    echo "Found $count files"
    
    # 生物信息学常用函数示例
    calculate_gc_content() {
        local fasta_file=$1
        local total=$(grep -v "^>" "$fasta_file" | tr -d '\n' | wc -c)
        local gc=$(grep -v "^>" "$fasta_file" | tr -d '\n' | grep -o "[GCgc]" | wc -l)
        
        if [ $total -gt 0 ]; then
            echo "scale=2; $gc * 100 / $total" | bc
        else
            echo "0"
        fi
    }
    
    gc_percent=$(calculate_gc_content "genome.fasta")
    echo "GC content: $gc_percent%"

输入输出
~~~~~~~~

Shell脚本需要与用户交互或处理输入输出。read命令用于从标准输入读取数据，echo和printf用于输出数据。

read命令从标准输入读取一行，并将输入赋值给变量。基本格式是read variable。如果输入多个词，可以使用多个变量，read var1 var2，多余的词会赋值给最后一个变量。read -p "prompt" variable显示提示信息后读取输入。read -s variable在输入时不显示字符，适合读取密码。read -t timeout variable设置超时时间，如果用户在指定时间内没有输入，命令返回非零状态。read -a array将输入读入数组。

echo命令是最简单的输出命令，它将参数打印到标准输出。echo "Hello World"输出字符串。echo支持一些转义字符，如\n表示换行，\t表示制表符，需要使用-e选项启用转义。echo -n选项不输出末尾的换行符。echo命令简单易用，但对输出格式的控制有限。

printf命令提供更强大的格式化输出功能，语法类似于C语言的printf函数。printf format arguments，格式字符串中的格式说明符会被参数替换。常见的格式说明符包括%s表示字符串，%d表示整数，%f表示浮点数。可以指定宽度和精度，例如%10s表示字符串占10个字符宽度，右对齐，%-10s表示左对齐，%.2f表示保留两位小数。printf不会自动添加换行符，需要在格式字符串中明确指定。

.. code-block:: bash

    # read命令基本用法
    echo "Enter your name:"
    read name
    echo "Hello, $name!"
    
    # 带提示信息的read
    read -p "Enter file path: " filepath
    echo "You entered: $filepath"
    
    # 读取多个变量
    read -p "Enter first name and last name: " first last
    echo "First: $first, Last: $last"
    
    # 读取密码（不显示输入）
    read -s -p "Enter password: " password
    echo ""
    echo "Password received"
    
    # 设置超时
    if read -t 5 -p "Enter value (5s timeout): " value; then
        echo "You entered: $value"
    else
        echo "Timeout!"
    fi
    
    # 读取数组
    echo "Enter multiple values:"
    read -a values
    echo "First value: ${values[0]}"
    echo "All values: ${values[@]}"
    
    # echo命令示例
    echo "Simple output"
    echo -n "No newline at end"
    echo " - continued"
    
    # echo使用转义字符
    echo -e "Line1\nLine2\nLine3"
    echo -e "Column1\tColumn2\tColumn3"
    
    # printf命令示例
    printf "Name: %s\n" "John"
    printf "Age: %d\n" 25
    printf "Score: %.2f\n" 95.567
    
    # 格式化表格输出
    printf "%-15s %10s %10s\n" "Name" "Score" "Grade"
    printf "%-15s %10d %10s\n" "Alice" 95 "A"
    printf "%-15s %10d %10s\n" "Bob" 87 "B"
    printf "%-15s %10d %10s\n" "Charlie" 92 "A"
    
    # 生物信息学示例：格式化输出序列统计
    printf "%-20s %10s %10s %10s\n" "Sample" "Reads" "Bases" "GC%"
    printf "%-20s %10d %10d %10.2f\n" "sample1.fastq" 10000 1500000 45.32
    printf "%-20s %10d %10d %10.2f\n" "sample2.fastq" 8500 1275000 48.76
    printf "%-20s %10d %10d %10.2f\n" "sample3.fastq" 12000 1800000 46.89
    
    # 交互式脚本示例
    read -p "Enter input FASTQ file: " input_file
    read -p "Enter output directory: " output_dir
    
    if [ ! -f "$input_file" ]; then
        echo "Error: Input file not found!"
        exit 1
    fi
    
    mkdir -p "$output_dir"
    echo "Processing $input_file..."
    echo "Results will be saved to $output_dir"

特殊变量
~~~~~~~~

Shell脚本中有一些特殊变量，它们由shell自动设置，提供了脚本执行环境的重要信息。理解这些特殊变量对于编写健壮的脚本至关重要。

$0是脚本本身的名称，包含脚本的路径。这在脚本需要知道自己被如何调用时很有用，例如显示用法信息时包含脚本名。$1到$9是脚本的位置参数，表示传递给脚本的第一个到第九个参数。如果参数超过9个，需要使用${10}、${11}等格式。$#是位置参数的个数，可以用来检查脚本是否收到了足够的参数。$@和$*都表示所有的位置参数，但在双引号中有区别。"$@"将每个参数作为独立的字符串，而"$*"将所有参数作为一个字符串。

$?是上一个命令的退出状态。在脚本中经常需要检查上一个命令是否成功执行，if [ $? -eq 0 ]; then表示上一个命令执行成功。$$是当前shell进程的进程ID，常用于生成临时文件名，例如temp_file="/tmp/script_$$"。$!是最近一个后台进程的进程ID，在启动后台任务后可以用这个变量监控或控制后台进程。

.. code-block:: bash

    # 创建示例脚本 demo.sh
    cat > demo.sh << 'EOF'
    #!/bin/bash
    
    # $0 - 脚本名称
    echo "Script name: $0"
    
    # $# - 参数个数
    echo "Number of arguments: $#"
    
    # 检查参数个数
    if [ $# -lt 2 ]; then
        echo "Usage: $0 <input_file> <output_file>"
        exit 1
    fi
    
    # $1, $2, ... - 位置参数
    echo "First argument: $1"
    echo "Second argument: $2"
    
    # ${10}, ${11}, ... - 第10个及以后的参数
    echo "Tenth argument: ${10}"
    
    # $@ - 所有参数（每个参数独立）
    echo "All arguments (using \$@):"
    for arg in "$@"; do
        echo "  - $arg"
    done
    
    # $* - 所有参数（作为一个字符串）
    echo "All arguments (using \$*): $*"
    
    # $$ - 当前进程ID
    echo "Process ID: $$"
    temp_file="/tmp/script_$$_temp.txt"
    echo "Temp file: $temp_file"
    
    # $? - 上一个命令的退出状态
    ls /nonexistent 2>/dev/null
    if [ $? -ne 0 ]; then
        echo "Command failed"
    fi
    
    # 成功的命令
    echo "Hello"
    if [ $? -eq 0 ]; then
        echo "Previous command succeeded"
    fi
    EOF
    
    chmod +x demo.sh
    
    # 运行脚本
    ./demo.sh file1.txt file2.txt file3.txt
    
    # 生物信息学脚本示例
    cat > process_samples.sh << 'EOF'
    #!/bin/bash
    
    # 检查参数
    if [ $# -eq 0 ]; then
        echo "Usage: $0 <sample1> <sample2> ..."
        exit 1
    fi
    
    echo "Processing $# samples..."
    
    # 使用$@遍历所有样本
    for sample in "$@"; do
        echo "Processing sample: $sample"
        
        # 创建临时文件
        temp_log="/tmp/process_$$_${sample}.log"
        
        # 执行处理
        fastqc "${sample}.fastq" > "$temp_log" 2>&1
        
        # 检查执行状态
        if [ $? -eq 0 ]; then
            echo "  Success: $sample"
        else
            echo "  Failed: $sample (check $temp_log)"
        fi
    done
    
    echo "All samples processed"
    EOF
    
    # 后台进程示例
    sleep 100 &
    background_pid=$!
    echo "Background process started with PID: $background_pid"
    
    # 检查后台进程是否还在运行
    if ps -p $background_pid > /dev/null; then
        echo "Process $background_pid is still running"
        kill $background_pid
    fi

作业控制与后台运行
~~~~~~~~~~~~~~~~~~

Shell支持作业控制，允许在后台运行任务，并在前台和后台之间切换。这对于运行长时间的数据分析任务非常有用。

在命令后加上&符号，命令会在后台运行，shell立即返回提示符，你可以继续输入其他命令。后台任务会输出一个作业号和进程ID。例如，python analysis.py &会在后台运行分析脚本。后台任务的输出仍然会显示在终端，可以重定向到文件避免干扰。

jobs命令列出当前shell会话中的所有作业。作业号用方括号括起来，后面是作业状态和命令。作业状态可以是Running、Stopped或Done。使用fg命令将后台作业调到前台运行，fg %job_number指定作业号。使用bg命令在后台继续运行暂停的作业。

Ctrl+Z可以暂停当前前台任务，任务会进入Stopped状态。然后可以使用bg命令让它在后台继续运行，或使用fg命令调回前台。这对于临时暂停一个长时间运行的任务，执行其他操作后再继续很有用。

wait命令等待后台任务完成。wait等待所有后台任务，wait PID等待指定进程。这在脚本中启动多个后台任务，然后等待它们全部完成时很有用。例如，并行处理多个文件，然后等待所有处理完成。

nohup命令让进程在用户退出登录后继续运行。nohup command &启动的进程会忽略挂起信号。默认情况下，nohup将输出重定向到nohup.out文件。建议明确指定输出文件，nohup command > output.log 2>&1 &。

.. code-block:: bash

    # 后台运行命令
    python analysis.py &
    # 输出: [1] 12345
    
    # 查看当前作业
    jobs
    # 输出: [1]+  Running                 python analysis.py &
    
    # 将后台作业调到前台
    fg %1
    
    # 暂停当前前台任务（Ctrl+Z）
    # 然后在后台继续运行
    bg %1
    
    # 启动多个后台任务
    fastqc sample1.fastq &
    fastqc sample2.fastq &
    fastqc sample3.fastq &
    
    # 等待所有后台任务完成
    wait
    
    # 等待特定进程
    sleep 100 &
    pid=$!
    wait $pid
    echo "Process $pid completed"
    
    # nohup示例
    nohup python long_analysis.py > analysis.log 2>&1 &
    
    # 查看nohup进程
    ps aux | grep long_analysis.py
    
    # 生物信息学并行处理示例
    cat > parallel_qc.sh << 'EOF'
    #!/bin/bash
    
    # 并行运行FastQC
    for fastq in *.fastq; do
        echo "Starting FastQC for $fastq"
        fastqc "$fastq" > "${fastq}.log" 2>&1 &
    done
    
    # 等待所有FastQC完成
    echo "Waiting for all FastQC jobs to complete..."
    wait
    
    echo "All FastQC jobs completed"
    
    # 检查结果
    ls -lh *fastqc.zip
    EOF
    
    # 后台运行并记录PID
    nohup bash parallel_qc.sh > parallel_qc.log 2>&1 &
    echo $! > parallel_qc.pid
    echo "Job started with PID: $(cat parallel_qc.pid)"
    
    # 监控后台任务
    tail -f parallel_qc.log
    
    # 检查进程是否还在运行
    if ps -p $(cat parallel_qc.pid) > /dev/null 2>&1; then
        echo "Job is still running"
    else
        echo "Job completed"
    fi

调试技巧
~~~~~~~~

调试是脚本开发中不可避免的环节。Shell提供了一些选项和技巧来帮助调试脚本。

set -x选项开启执行跟踪，shell会在执行每个命令前打印该命令，命令前会有加号标记。这对于查看脚本的实际执行流程非常有用。可以在脚本开头使用set -x，也可以在脚本的特定部分前后开启和关闭。set +x关闭执行跟踪。例如，在脚本中插入set -x和set +x可以只跟踪特定部分的执行。

set -e选项让脚本在命令返回非零状态时立即退出。默认情况下，脚本会继续执行后续命令，即使前面的命令失败了。set -e可以防止错误累积，让问题尽早暴露。但有时命令的失败是预期的，可以使用command || true来避免触发退出，或者临时关闭set -e。

set -u选项将未定义变量的使用视为错误。默认情况下，使用未定义的变量会返回空字符串，这可能导致难以发现的错误。set -u可以帮助发现变量名拼写错误等问题。如果某个变量可能未定义，可以使用${variable:-default}提供默认值。

set -o pipefail选项让管道的退出状态为最后一个返回非零状态的命令的退出状态。默认情况下，管道的退出状态是最后一个命令的退出状态，即使前面的命令失败了。这对于检测管道中的错误很有用。

可以在脚本开头组合使用这些选项，set -euxo pipefail开启所有常用的调试选项。也可以在执行脚本时通过命令行选项指定，bash -x script.sh开启执行跟踪。

.. code-block:: bash

    # 创建带调试选项的脚本
    cat > debug_demo.sh << 'EOF'
    #!/bin/bash
    
    # 开启所有调试选项
    set -euxo pipefail
    
    # 脚本内容
    file="data.txt"
    
    # 这个命令会失败，因为文件不存在
    cat "$file"
    
    echo "This line won't execute due to set -e"
    EOF
    
    # 运行脚本（会失败并显示详细信息）
    bash debug_demo.sh
    
    # 只开启执行跟踪
    cat > trace_demo.sh << 'EOF'
    #!/bin/bash
    
    echo "Starting script"
    
    set -x  # 开启跟踪
    for i in {1..3}; do
        echo "Number: $i"
    done
    set +x  # 关闭跟踪
    
    echo "Script finished"
    EOF
    
    bash trace_demo.sh
    
    # 使用set -e防止错误累积
    cat > safe_script.sh << 'EOF'
    #!/bin/bash
    set -e
    
    # 如果任何命令失败，脚本立即退出
    mkdir -p results
    cd results
    
    # 这个命令可能失败，但我们想继续
    test -f config.txt || echo "Config not found, using defaults"
    
    # 或者临时关闭set -e
    set +e
    some_risky_command
    set -e
    
    echo "Processing completed"
    EOF
    
    # 使用set -u检测未定义变量
    cat > check_undefined.sh << 'EOF'
    #!/bin/bash
    set -u
    
    name="Alice"
    echo "Hello, $name"
    
    # 这行会报错，因为变量未定义
    echo "Hello, $undefined_var"
    EOF
    
    # 使用默认值避免未定义变量错误
    cat > safe_variables.sh << 'EOF'
    #!/bin/bash
    
    # 使用默认值
    input_file="${input_file:-default.txt}"
    output_dir="${output_dir:-output}"
    threads="${threads:-4}"
    
    echo "Input: $input_file"
    echo "Output: $output_dir"
    echo "Threads: $threads"
    EOF
    
    # 使用pipefail检测管道错误
    cat > pipe_check.sh << 'EOF'
    #!/bin/bash
    set -o pipefail
    
    # 如果grep失败，整个管道会返回非零状态
    cat data.txt | grep "pattern" | wc -l
    
    if [ $? -eq 0 ]; then
        echo "Pattern found"
    else
        echo "Pattern not found or file error"
    fi
    EOF
    
    # 命令行调试选项
    # 执行跟踪
    bash -x script.sh
    
    # 检查语法错误
    bash -n script.sh
    
    # 详细模式
    bash -v script.sh
    
    # 组合调试脚本
    cat > bioinfo_debug.sh << 'EOF'
    #!/bin/bash
    # 调试选项：错误即退出、未定义变量报错、执行跟踪、管道错误检测
    set -euxo pipefail
    
    input_fastq="${1:?Error: Input FASTQ file required}"
    output_dir="${2:-results}"
    
    echo "Processing $input_fastq"
    mkdir -p "$output_dir"
    
    # FastQC质量控制
    fastqc "$input_fastq" -o "$output_dir"
    
    # 检查结果
    if [ -f "${output_dir}/${input_fastq}_fastqc.html" ]; then
        echo "FastQC completed successfully"
    else
        echo "FastQC failed"
        exit 1
    fi
    EOF

常用单行命令
~~~~~~~~~~~~

Shell命令行工具提供了丰富的文本处理能力，熟练使用这些工具可以快速完成许多数据处理任务。这些工具在生物信息学数据处理中经常使用。

awk是一个强大的文本处理工具，它可以按列处理文本文件。awk的基本格式是awk 'pattern {action}' file。pattern是匹配条件，action是要执行的操作。例如，awk '{print $1}' file打印每行的第一列。awk -F, '{print $1,$3}' file使用逗号作为分隔符，打印第一列和第三列。awk '$3 > 100 {print $0}' file打印第三列大于100的行。awk可以计算统计值，awk '{sum+=$1} END {print sum}' file计算第一列的总和。awk支持变量、条件判断、循环等编程结构，可以完成复杂的文本处理任务。

sed是流编辑器，用于文本替换和转换。sed的基本格式是sed 'command' file。最常用的命令是替换，sed 's/old/new/g' file将所有old替换为new。g标志表示全局替换，不加g只替换每行的第一个匹配。sed 's/old/new/2' file只替换每行的第二个匹配。sed可以删除行，sed '/pattern/d' file删除匹配pattern的行。sed -n '10,20p' file打印第10到20行。sed -i选项直接修改文件，sed -i 's/old/new/g' file原地修改文件。

cut命令从每行提取指定的列或字符。cut -d, -f1,3 file使用逗号作为分隔符，提取第一列和第三列。cut -c1-10 file提取每行的第1到第10个字符。cut命令简单快速，适合提取固定格式的数据。

sort命令对文本文件排序。sort file按字典序排序。sort -n file按数值排序。sort -k2,2 file按第二列排序。sort -r file逆序排序。sort -u file排序并去重。sort常与其他命令配合使用，例如提取唯一值。

uniq命令去除连续的重复行。uniq file去除连续的重复行，只保留一个。uniq -c file显示每行出现的次数。uniq -d file只显示重复的行。uniq通常与sort配合使用，sort file | uniq先排序再去重。

wc命令统计文件的行数、词数和字节数。wc -l file统计行数。wc -w file统计词数。wc -c file统计字节数。wc常用于快速查看文件大小或记录数。

head和tail命令查看文件的开头或结尾部分。head -n 20 file显示文件的前20行。tail -n 20 file显示文件的后20行。tail -f file持续显示文件末尾的新增内容，常用于查看日志文件。

paste命令将多个文件按列合并。paste file1 file2将两个文件并列显示，用制表符分隔。paste -d, file1 file2使用逗号作为分隔符。paste适合合并相关的数据文件。

join命令根据共同字段连接两个文件。join -1 1 -2 2 file1 file2根据file1的第一列和file2的第二列连接。join要求输入文件已按连接字段排序。join类似于数据库的连接操作。

.. code-block:: bash

    # awk命令示例
    
    # 打印第一列
    awk '{print $1}' data.txt
    
    # 指定分隔符，打印多列
    awk -F, '{print $1, $3}' data.csv
    
    # 条件过滤
    awk '$3 > 100 {print $0}' data.txt
    
    # 计算统计值
    awk '{sum+=$1} END {print "Sum:", sum, "Average:", sum/NR}' numbers.txt
    
    # 生物信息学示例：计算FASTQ文件的平均长度
    awk 'NR%4==2 {sum+=length($0); count++} END {print "Average read length:", sum/count}' sample.fastq
    
    # 提取FASTA序列名
    awk '/^>/ {print substr($1, 2)}' genome.fasta
    
    # 统计每条染色体上的基因数量（GFF文件）
    awk '$3=="gene" {count[$1]++} END {for (chr in count) print chr, count[chr]}' genes.gff
    
    # sed命令示例
    
    # 替换所有匹配
    sed 's/old/new/g' file.txt
    
    # 只替换每行第一个匹配
    sed 's/old/new/' file.txt
    
    # 替换每行第二个匹配
    sed 's/old/new/2' file.txt
    
    # 删除匹配行
    sed '/pattern/d' file.txt
    
    # 打印指定行范围
    sed -n '10,20p' file.txt
    
    # 原地修改文件
    sed -i 's/old/new/g' file.txt
    
    # 生物信息学示例：修改染色体名称
    sed 's/^>chr/>/g' genome.fasta
    
    # cut命令示例
    
    # 提取指定列（制表符分隔）
    cut -f1,3 data.txt
    
    # 指定分隔符
    cut -d, -f1,3 data.csv
    
    # 提取指定字符范围
    cut -c1-10 data.txt
    
    # 提取FASTQ质量值（第4行）
    sed -n '4~4p' sample.fastq | cut -c1-50
    
    # sort命令示例
    
    # 字典序排序
    sort data.txt
    
    # 数值排序
    sort -n numbers.txt
    
    # 按指定列排序
    sort -k2,2 data.txt
    
    # 逆序排序
    sort -r data.txt
    
    # 排序并去重
    sort -u data.txt
    
    # 按数值逆序排序
    sort -k2,2nr data.txt
    
    # uniq命令示例
    
    # 去除连续重复行
    uniq data.txt
    
    # 显示重复次数
    uniq -c data.txt
    
    # 只显示重复行
    uniq -d data.txt
    
    # 统计唯一值（需要先排序）
    sort data.txt | uniq -c | sort -nr
    
    # wc命令示例
    
    # 统计行数
    wc -l data.txt
    
    # 统计词数
    wc -w data.txt
    
    # 统计字节数
    wc -c data.txt
    
    # 统计FASTQ文件的reads数量
    wc -l sample.fastq | awk '{print $1/4}'
    
    # head和tail命令示例
    
    # 查看文件前20行
    head -n 20 data.txt
    
    # 查看文件后20行
    tail -n 20 data.txt
    
    # 实时查看日志
    tail -f analysis.log
    
    # 查看FASTA文件前10条序列
    head -n 40 genome.fasta
    
    # paste命令示例
    
    # 合并两个文件
    paste file1.txt file2.txt
    
    # 指定分隔符
    paste -d',' file1.txt file2.txt
    
    # 合并多个文件
    paste sample1.txt sample2.txt sample3.txt > combined.txt
    
    # join命令示例
    
    # 根据第一列连接
    join file1.txt file2.txt
    
    # 指定不同的连接列
    join -1 1 -2 2 file1.txt file2.txt
    
    # 指定输出字段
    join -o 1.1,2.2 file1.txt file2.txt
    
    # 生物信息学综合示例
    
    # 统计FASTQ文件中reads的平均长度和质量
    awk 'NR%4==2 {len+=length($0); count++} NR%4==0 {qual+=length($0)} END {
        print "Total reads:", count
        print "Average length:", len/count
        print "Average quality length:", qual/count
    }' sample.fastq
    
    # 从GFF文件中提取基因ID和位置信息
    awk '$3=="gene" {
        match($9, /ID=([^;]+)/, arr);
        print arr[1]"\t"$1"\t"$4"\t"$5
    }' genes.gff | sort -k1,1
    
    # 统计VCF文件中每个样本的变异数量
    grep -v "^#" variants.vcf | cut -f10- | awk '{
        for(i=1; i<=NF; i++) {
            if($i !~ /^\./) count[i]++
        }
    } END {
        for(i in count) print "Sample", i, "Variants:", count[i]
    }'
    
    # 提取特定染色体的序列
    samtools faidx genome.fasta chr1 > chr1.fasta
    
    # 批量重命名文件
    for f in *.fastq; do
        mv "$f" "$(echo $f | sed 's/.fastq/_R1.fastq/')"
    done

环境与包管理
------------

Conda / Mamba 安装与配置
~~~~~~~~~~~~~~~~~~~~~~~~~

在生物信息学工作中，我们需要使用大量的软件工具和编程库，这些软件往往有不同的依赖关系和版本要求。Conda是一个开源的包管理系统和环境管理系统，它可以方便地安装、运行和更新软件包及其依赖。Mamba是Conda的快速替代品，使用C++实现，在解析依赖关系时速度更快。

Conda的安装相对简单，推荐安装Miniconda，它是Conda的最小安装版本，只包含Conda和Python。从官网下载对应系统的安装脚本，执行安装程序并按照提示操作即可。安装完成后，需要重新打开终端或执行source命令使配置生效。Conda安装后会自动配置shell环境，添加必要的初始化代码到shell配置文件中。

Mamba可以作为Conda的替代品安装，也可以在现有的Conda环境中安装。Mamba完全兼容Conda的命令和配置，可以作为drop-in替代品使用。安装Mamba后，可以将conda命令替换为mamba命令，享受更快的包解析和安装速度。在处理复杂的依赖关系时，Mamba的优势尤为明显。

Conda的配置文件位于用户主目录下的.condarc文件中，可以配置镜像源、代理设置、默认频道等选项。国内用户可以配置清华、中科大等镜像源，加快下载速度。配置镜像源可以通过命令行或直接编辑配置文件完成。合理的配置可以大大提高软件安装的成功率和速度。

环境创建、激活、导出
~~~~~~~~~~~~~~~~~~~~

Conda的核心功能之一是环境管理，它允许创建相互隔离的软件环境，每个环境可以有不同的Python版本和软件包。这种隔离机制避免了不同项目之间的依赖冲突，是保证分析可重复性的重要手段。

创建新环境使用conda create命令，需要指定环境名称和要安装的包。例如，conda create -n bioinfo python=3.9创建一个名为bioinfo的环境，安装Python 3.9。创建环境时可以同时安装多个包，conda create -n analysis python=3.9 numpy pandas matplotlib。环境创建后，会在Conda的envs目录下创建对应的目录，所有该环境的软件都安装在这个目录中。

激活环境使用conda activate命令，conda activate bioinfo激活名为bioinfo的环境。激活后，命令行提示符通常会显示当前环境名称，PATH环境变量也会相应调整，使得该环境中的命令优先被找到。在激活的环境中安装的软件只在该环境中可用，不会影响其他环境。退出当前环境使用conda deactivate命令。

查看所有环境使用conda env list或conda info --envs命令，会列出所有环境及其路径。当前激活的环境会有星号标记。删除环境使用conda env remove -n envname或conda remove -n envname --all命令。删除环境会删除该环境目录及其中的所有软件，释放磁盘空间。

导出环境配置对于分享和复现分析环境非常重要。conda env export > environment.yml将当前环境的所有包及其版本导出到YAML文件。这个文件包含了环境名称、频道信息和所有依赖包的精确版本。其他人可以使用conda env create -f environment.yml从配置文件创建完全相同的环境。也可以导出为文本格式，conda list --export > requirements.txt，然后使用conda create --name newenv --file requirements.txt创建环境。

频道管理
~~~~~~~~

Conda的软件包存放在不同的频道中，频道是软件包的仓库。默认频道是defaults，包含由Anaconda公司维护的常用软件包。Conda-forge是一个社区驱动的频道，包含大量科学计算和数据分析相关的软件包，更新频率高，版本较新。Bioconda是专门为生物信息学软件设立的频道，包含了数千个生物信息学工具。

频道的配置影响软件包的搜索和安装顺序。使用conda config --add channels channelname添加频道，频道会按照添加顺序的逆序搜索，最后添加的频道优先级最高。建议的频道顺序是conda-forge优先，然后是bioconda，最后是defaults。这个顺序确保了软件包的最新版本优先被安装。

频道优先级可以在配置文件中查看和修改。conda config --show channels显示当前配置的频道。conda config --remove channels channelname移除指定频道。conda config --set channel_priority strict设置严格优先级模式，确保高优先级频道的包不会被低优先级频道的同名包覆盖。严格优先级模式可以避免版本冲突，但可能导致某些包无法安装。

Bioconda频道是生物信息学工作者必须配置的频道。它包含了大量常用的生物信息学软件，如samtools、bcftools、bwa、hisat2等，以及R语言的Bioconductor包。使用Bioconda需要先配置conda-forge频道，因为Bioconda的许多依赖包来自conda-forge。正确配置频道后，安装生物信息学软件变得非常简单，conda install package-name即可。

虚拟环境与隔离
~~~~~~~~~~~~~~

除了Conda，Python还有其他虚拟环境管理工具，如venv和virtualenv。这些工具创建的虚拟环境更加轻量，只管理Python包，不管理系统级的软件包。对于纯Python项目，使用这些工具就足够了。

venv是Python标准库的一部分，从Python 3.3开始提供。创建虚拟环境使用python -m venv envname命令，会在当前目录下创建envname目录，包含Python解释器和pip包管理器。激活环境使用source envname/bin/activate命令，退出使用deactivate命令。在虚拟环境中安装的Python包存放在虚拟环境目录中，与系统Python隔离。

virtualenv是一个第三方工具，功能比venv更强大，支持Python 2和Python 3，可以为已有的Python解释器创建虚拟环境。安装使用pip install virtualenv命令。创建环境使用virtualenv envname，激活和退出方式与venv相同。virtualenv还支持指定Python版本，virtualenv -p python3.8 envname使用指定的Python版本创建环境。

选择哪种虚拟环境工具取决于项目需求。如果项目只需要Python包，venv或virtualenv就足够了。如果项目需要系统级的软件包，如C库或生物信息学工具，Conda是更好的选择。在实际工作中，Conda因其强大的跨语言包管理能力，在生物信息学领域使用最广泛。

容器化基础
~~~~~~~~~~

容器化技术提供了比虚拟环境更彻底的隔离机制。容器不仅隔离软件包，还隔离整个运行环境，包括操作系统、库文件、环境变量等。容器技术将在后续章节详细介绍，这里简要说明其与环境管理的关系。

Docker是最流行的容器化平台，它将应用程序及其所有依赖打包成镜像，镜像可以分发到任何支持Docker的系统上运行。容器化确保了分析环境的完全可重复性，无论在什么系统上运行，结果都是一致的。许多生物信息学软件都提供了Docker镜像，可以直接使用。

Singularity是另一种容器技术，特别适合在高性能计算环境中使用。与Docker需要root权限不同，Singularity允许普通用户运行容器，这使得它在共享的计算集群中更受欢迎。Singularity可以使用Docker镜像，也可以创建自己的镜像格式。

容器化技术与Conda可以结合使用。在Docker或Singularity镜像中使用Conda管理软件包，既保证了环境的可移植性，又保留了Conda方便的包管理能力。许多生物信息学的容器镜像都采用这种方式构建。

.. code-block:: bash

    # Conda安装和配置示例
    
    # 下载并安装Miniconda（Linux/macOS）
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    bash Miniconda3-latest-Linux-x86_64.sh
    
    # 初始化conda
    conda init bash
    source ~/.bashrc
    
    # 配置镜像源（国内用户）
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
    conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
    
    # 设置搜索时显示频道地址
    conda config --set show_channel_urls yes
    
    # 查看配置
    conda config --show channels
    
    # 环境管理示例
    
    # 创建新环境
    conda create -n bioinfo python=3.9
    conda create -n rna_seq python=3.9 numpy pandas matplotlib
    
    # 激活环境
    conda activate bioinfo
    
    # 在环境中安装软件
    conda install -c bioconda samtools
    conda install -c bioconda bwa
    conda install -c bioconda fastqc
    
    # 查看已安装的包
    conda list
    
    # 查看所有环境
    conda env list
    conda info --envs
    
    # 导出环境配置
    conda env export > environment.yml
    
    # 从配置文件创建环境
    conda env create -f environment.yml
    
    # 导出为文本格式
    conda list --export > requirements.txt
    conda create --name newenv --file requirements.txt
    
    # 删除环境
    conda deactivate
    conda env remove -n bioinfo
    
    # Mamba使用示例
    
    # 安装Mamba
    conda install -c conda-forge mamba
    
    # 使用Mamba（命令与conda相同，但速度更快）
    mamba create -n analysis python=3.9
    mamba install -c bioconda star hisat2
    
    # 生物信息学环境配置示例
    
    # 创建RNA-seq分析环境
    conda create -n rna_seq_analysis python=3.9
    conda activate rna_seq_analysis
    conda install -c bioconda fastqc trimmomatic star hisat2 samtools subread
    conda install -c bioconda bioconductor-deseq2 bioconductor-edger
    
    # 创建基因组分析环境
    conda create -n genome_analysis python=3.9
    conda activate genome_analysis
    conda install -c bioconda bwa samtools bcftools gatk4 snpeff
    
    # 创建Python数据分析环境
    conda create -n pydata python=3.9 numpy pandas matplotlib seaborn scikit-learn jupyter
    
    # venv虚拟环境示例
    
    # 创建虚拟环境
    python -m venv myproject_env
    
    # 激活环境（Linux/macOS）
    source myproject_env/bin/activate
    
    # 激活环境（Windows）
    myproject_env\Scripts\activate
    
    # 安装包
    pip install numpy pandas matplotlib
    
    # 导出依赖
    pip freeze > requirements.txt
    
    # 从requirements.txt安装
    pip install -r requirements.txt
    
    # 退出环境
    deactivate

Python编程基础
--------------

Python是一种高级编程语言，以其简洁清晰的语法和强大的功能而闻名。在生物信息学领域，Python是最常用的编程语言之一，广泛用于数据处理、分析和可视化。Python的设计哲学强调代码的可读性和简洁性，这使得它成为初学者学习编程的理想选择，同时也足够强大，能够处理复杂的科学计算任务。

Python的基本语法
~~~~~~~~~~~~~~~~

Python的语法设计追求简洁和可读性。与许多其他编程语言使用大括号来表示代码块不同，Python使用缩进来表示代码块的层次结构。这种设计强制程序员编写格式整齐的代码，提高了代码的可读性。通常建议使用4个空格作为一级缩进，避免使用制表符以防止在不同编辑器中显示不一致。

变量是存储数据的容器。Python是动态类型语言，变量不需要预先声明类型，解释器会根据赋值自动推断类型。变量名应该具有描述性，使用小写字母，单词之间用下划线分隔。Python支持多种数据类型，包括整数、浮点数、字符串、布尔值等基本类型，以及列表、元组、字典、集合等复合类型。

注释是代码中不被执行的部分，用于解释代码的功能和逻辑。Python使用井号表示单行注释，三个引号表示多行注释或文档字符串。良好的注释习惯是编写可维护代码的重要组成部分，注释应该解释为什么这样做，而不是做了什么。

Python的运算符包括算术运算符、比较运算符、逻辑运算符和赋值运算符等。算术运算符包括加、减、乘、除、整除、取余和幂运算。比较运算符用于比较两个值，返回布尔值。逻辑运算符包括and、or和not，用于组合布尔表达式。赋值运算符用于给变量赋值，还支持复合赋值运算符如+=、-=等。

数据类型与结构
~~~~~~~~~~~~~~

Python提供了丰富的内置数据类型。数值类型包括整数和浮点数，整数可以是任意大小，浮点数是双精度浮点数。字符串是不可变的字符序列，可以使用单引号或双引号定义，支持转义字符和原始字符串。布尔类型只有True和False两个值，常用于条件判断。

列表是Python中最常用的数据结构，它是一个有序的可变序列，可以包含不同类型的元素。列表使用方括号定义，元素之间用逗号分隔。列表支持索引和切片操作，索引从0开始，负数索引从末尾开始计数。列表提供了丰富的方法，如append添加元素、extend扩展列表、insert插入元素、remove删除元素、pop弹出元素等。

元组与列表类似，但是不可变的。元组使用圆括号定义，一旦创建就不能修改。元组通常用于存储不应该被改变的数据，如坐标点、配置信息等。由于元组不可变，它可以作为字典的键，而列表不能。元组在内存占用和访问速度上比列表更高效。

字典是键值对的集合，使用花括号定义。字典是无序的，通过键来访问值。键必须是不可变类型，如字符串、数字或元组，值可以是任意类型。字典提供了丰富的方法，如keys获取所有键、values获取所有值、items获取所有键值对、get安全获取值等。字典在生物信息学中广泛用于存储映射关系，如基因ID到基因名的映射。

集合是无序的唯一元素集合，使用花括号或set函数创建。集合支持数学集合运算，如并集、交集、差集等。集合常用于去重和成员测试，成员测试的时间复杂度是O(1)，比列表快得多。

控制流程
~~~~~~~~

控制流程语句决定了程序的执行顺序。条件语句使用if、elif和else关键字，根据条件执行不同的代码块。条件表达式可以是任何能够转换为布尔值的表达式，非零数值、非空字符串、非空列表等都被视为True。条件语句可以嵌套，但应该避免过深的嵌套层次。

循环语句用于重复执行代码块。for循环遍历可迭代对象，如列表、字符串、字典等。for循环可以与range函数配合使用，生成指定范围的数字序列。while循环在条件为真时重复执行代码块，需要注意避免无限循环，确保循环条件最终会变为假。

循环控制语句break和continue用于改变循环的执行流程。break语句立即退出循环，不再执行循环体中剩余的语句。continue语句跳过当前迭代的剩余语句，直接进入下一次迭代。这些语句应该谨慎使用，过度使用会降低代码的可读性。

Python提供了列表推导式、字典推导式和集合推导式，可以用简洁的语法创建新的数据结构。推导式是Python的特色语法，它将循环和条件判断组合在一行中，使代码更加简洁。但复杂的推导式会降低可读性，应该在简洁和可读性之间找到平衡。

函数定义与调用
~~~~~~~~~~~~~~

函数是组织好的、可重复使用的代码块，用于执行特定的任务。函数使用def关键字定义，后跟函数名和参数列表。函数名应该具有描述性，使用小写字母，单词之间用下划线分隔。参数列表可以为空，也可以包含多个参数，参数可以有默认值。

函数的参数有多种类型。位置参数是最基本的参数类型，按照位置顺序传递。关键字参数在调用时指定参数名，可以不按照定义顺序传递。默认参数在定义时指定默认值，调用时可以省略。可变参数使用星号接收任意数量的位置参数，双星号接收任意数量的关键字参数。

函数可以有返回值，使用return语句返回。如果没有return语句或return后没有值，函数返回None。函数可以返回多个值，实际上是返回一个元组。返回值可以被调用者接收和使用，也可以忽略。

变量作用域是指变量的可见范围。函数内部定义的变量是局部变量，只在函数内部可见。函数外部定义的变量是全局变量，可以在整个程序中访问。在函数内部可以使用global关键字声明使用全局变量，但这通常不是好的做法，应该尽量避免。

Python支持匿名函数，使用lambda关键字定义。匿名函数是一种简洁的函数定义方式，通常用于简单的、只用一次的函数。匿名函数的函数体只能是单个表达式，不能包含复杂语句。匿名函数常用于高阶函数的参数，如map、filter、sorted等函数的key参数。

模块与包
~~~~~~~~

模块是包含Python代码的文件，文件名就是模块名加上.py扩展名。模块可以定义函数、类和变量，也可以包含可执行的代码。使用import语句导入模块，可以使用模块中定义的函数和变量。导入时可以使用as关键字给模块起别名，使用from关键字导入模块中的特定对象。

包是包含多个模块的目录，目录中必须包含__init__.py文件。包提供了一种组织模块的方式，避免模块名冲突。包可以嵌套，形成层级结构。导入包中的模块可以使用点号表示层级关系，如import package.module。

Python标准库提供了大量有用的模块，如os模块用于操作系统接口、sys模块用于Python解释器相关功能、re模块用于正则表达式、json模块用于JSON数据处理等。标准库是Python安装的一部分，可以直接导入使用，无需额外安装。

第三方包是Python生态系统的重要组成部分。PyPI是Python包索引，包含了数十万个第三方包。使用pip工具可以方便地安装、升级和卸载第三方包。pip install package-name安装包，pip install package-name==version安装指定版本，pip uninstall package-name卸载包，pip list列出已安装的包。

文件操作
~~~~~~~~

文件操作是程序与外部数据交互的重要方式。Python使用open函数打开文件，返回文件对象。open函数接受文件路径和模式参数，模式包括只读r、写入w、追加a、读写r+等。对于文本文件，还需要指定编码方式，推荐使用UTF-8编码。

文件对象提供read、readline、readlines方法读取文件内容。read一次性读取整个文件，readline每次读取一行，readlines读取所有行并返回列表。对于大文件，应该使用迭代方式逐行读取，避免一次性加载到内存。

写入文件使用write方法，该方法返回写入的字符数。写入后需要调用close方法关闭文件，或者使用with语句自动管理文件资源。with语句是推荐的文件操作方式，它会在代码块结束时自动关闭文件，即使发生异常也能正确关闭。

Python还提供了os和shutil模块用于文件系统操作。os模块提供了文件和目录操作的函数，如创建目录、删除文件、重命名等。shutil模块提供了高级的文件操作，如复制文件和目录、移动文件、删除目录树等。这些模块在处理大量文件时非常有用。

异常处理
~~~~~~~~

异常是程序运行时发生的错误，如果不处理，程序会终止。Python使用try-except语句捕获和处理异常。try块中放置可能引发异常的代码，except块中放置异常处理代码。可以捕获特定类型的异常，也可以捕获所有异常。

except语句可以指定异常类型，多个except块可以处理不同类型的异常。可以使用as关键字将异常对象赋值给变量，获取异常的详细信息。except块后面可以有else块，当try块没有发生异常时执行。finally块无论是否发生异常都会执行，通常用于清理资源。

Python内置了许多异常类型，如ValueError、TypeError、FileNotFoundError等。程序也可以使用raise语句主动抛出异常，这在参数验证和错误处理时很有用。自定义异常类应该继承Exception类或其子类，按照惯例异常类名以Error结尾。

良好的异常处理应该遵循一些原则。不要使用空的except块捕获所有异常，这会隐藏错误。尽量捕获具体的异常类型，而不是笼统的Exception。在except块中应该记录错误信息或采取恢复措施，而不是简单地忽略。finally块用于释放资源，如关闭文件、数据库连接等。

面向对象编程
~~~~~~~~~~~~

面向对象编程是一种编程范式，它将数据和操作数据的方法封装在对象中。Python是一种多范式语言，支持面向对象编程。类是对象的模板，定义了对象的属性和方法。对象是类的实例，具有类定义的属性和方法。

类使用class关键字定义，后跟类名和冒号。类名通常使用驼峰命名法，首字母大写。类的属性在类体中定义，方法使用def关键字定义。方法的第一个参数通常是self，表示对象本身。通过self可以访问对象的属性和其他方法。

构造方法__init__在创建对象时自动调用，用于初始化对象的属性。析构方法__del__在对象被销毁时调用，但Python的垃圾回收机制使得析构方法的调用时机不确定，不应该依赖析构方法进行资源清理。

继承是面向对象编程的重要特性，子类可以继承父类的属性和方法。子类可以添加新的属性和方法，也可以重写父类的方法。Python支持多重继承，一个类可以继承多个父类。多重继承增加了灵活性，但也带来了复杂性，应该谨慎使用。

封装是隐藏对象的内部实现细节，只暴露必要的接口。Python通过命名约定实现封装，以单下划线开头的属性和方法表示受保护的，以双下划线开头的属性和方法表示私有的。这只是一种约定，Python并没有真正的访问控制机制。

多态是指不同类型的对象可以对相同的方法调用做出不同的响应。Python的鸭子类型使得多态变得自然，只要对象具有相应的方法，就可以调用，而不需要继承关系。这种灵活性是Python的优点，但也需要程序员注意类型安全。

生物信息学常用库
~~~~~~~~~~~~~~~~

Python在生物信息学领域有许多专门开发的库。Biopython是最著名的生物信息学Python库，提供了处理序列、结构、进化树等生物数据的工具。Biopython可以读写多种格式的生物数据文件，如FASTA、GenBank、PDB等，还可以访问NCBI、UniProt等在线数据库。

pandas是数据分析的核心库，提供了DataFrame数据结构，类似于R语言的data.frame。DataFrame是一个二维表格，每列可以是不同的数据类型。pandas提供了强大的数据操作功能，如数据清洗、转换、聚合、合并等。在生物信息学中，pandas常用于处理表达矩阵、样本信息表等表格数据。

NumPy是科学计算的基础库，提供了高效的多维数组对象和数学函数。NumPy数组在内存中是连续存储的，比Python列表更高效。NumPy提供了大量的数学函数，可以进行向量化运算，避免显式循环。在生物信息学中，NumPy常用于数值计算，如矩阵运算、统计分析等。

Matplotlib是Python的绘图库，可以生成各种静态、动态和交互式的图表。Matplotlib的设计参考了MATLAB，提供了类似MATLAB的绘图接口。Seaborn是基于Matplotlib的高级绘图库，提供了更美观的默认样式和更简单的接口。在生物信息学中，这些库常用于可视化基因表达、聚类结果、基因组特征等。

Scikit-learn是机器学习库，提供了各种机器学习算法的实现，包括分类、回归、聚类、降维等。Scikit-learn的API设计一致，易于使用。在生物信息学中，Scikit-learn常用于样本分类、特征选择、降维分析等任务。

.. code-block:: python

    # Python基本语法示例
    
    # 变量和数据类型
    gene_name = "TP53"
    expression_level = 45.6
    is_significant = True
    sample_ids = ["S1", "S2", "S3"]
    
    # 注释
    # 这是单行注释
    """
    这是多行注释
    可以跨越多行
    """
    
    # 运算符
    a, b = 10, 3
    print(a + b)   # 加法: 13
    print(a - b)   # 减法: 7
    print(a * b)   # 乘法: 30
    print(a / b)   # 除法: 3.333...
    print(a // b)  # 整除: 3
    print(a % b)   # 取余: 1
    print(a ** b)  # 幂运算: 1000
    
    # 比较运算符
    print(a > b)   # True
    print(a == b)  # False
    
    # 逻辑运算符
    print(True and False)  # False
    print(True or False)   # True
    print(not True)        # False

.. code-block:: python

    # 数据类型与结构示例
    
    # 列表操作
    genes = ["BRCA1", "TP53", "EGFR"]
    genes.append("MYC")
    genes.insert(0, "KRAS")
    print(genes)
    
    # 列表索引和切片
    print(genes[0])
    print(genes[-1])
    print(genes[1:3])
    
    # 元组
    coordinates = (100, 200, 300)
    print(coordinates[0])
    
    # 字典
    gene_info = {
        "name": "TP53",
        "chromosome": "chr17",
        "start": 7668421,
        "end": 7687490
    }
    print(gene_info["name"])
    print(gene_info.get("strand", "+"))
    
    # 集合
    unique_genes = set(["TP53", "BRCA1", "TP53", "EGFR"])
    print(unique_genes)
    
    # 生物信息学示例：基因表达数据
    expression_data = {
        "TP53": [10.5, 12.3, 8.7, 15.2],
        "BRCA1": [5.2, 6.1, 4.8, 7.3],
        "EGFR": [20.1, 18.5, 22.3, 19.8]
    }
    
    for gene, values in expression_data.items():
        avg = sum(values) / len(values)
        print(f"{gene}: average expression = {avg:.2f}")

.. code-block:: python

    # 控制流程示例
    
    # 条件语句
    expression = 15.5
    if expression > 20:
        status = "high"
    elif expression > 10:
        status = "medium"
    else:
        status = "low"
    print(f"Expression status: {status}")
    
    # for循环
    genes = ["TP53", "BRCA1", "EGFR"]
    for gene in genes:
        print(f"Processing gene: {gene}")
    
    # 使用range
    for i in range(5):
        print(f"Sample {i+1}")
    
    # while循环
    count = 0
    while count < 3:
        print(f"Count: {count}")
        count += 1
    
    # break和continue
    for gene in genes:
        if gene == "BRCA1":
            continue
        print(f"Analyzing {gene}")
        if gene == "EGFR":
            break
    
    # 列表推导式
    expression_values = [10.5, 12.3, 8.7, 15.2, 6.1]
    high_expression = [x for x in expression_values if x > 10]
    print(high_expression)
    
    # 字典推导式
    gene_lengths = {"TP53": 393, "BRCA1": 1863, "EGFR": 1210}
    long_genes = {k: v for k, v in gene_lengths.items() if v > 500}
    print(long_genes)

.. code-block:: python

    # 函数定义与调用示例
    
    # 基本函数
    def calculate_gc_content(sequence):
        gc_count = sequence.count('G') + sequence.count('C')
        return (gc_count / len(sequence)) * 100
    
    gc = calculate_gc_content("ATGCGATCGA")
    print(f"GC content: {gc:.2f}%")
    
    # 带默认参数的函数
    def analyze_sequence(sequence, min_length=10):
        if len(sequence) < min_length:
            return None
        return {
            "length": len(sequence),
            "gc_content": calculate_gc_content(sequence)
        }
    
    result = analyze_sequence("ATGCGATCGATCG")
    print(result)
    
    # 可变参数
    def calculate_mean(*values):
        return sum(values) / len(values)
    
    mean = calculate_mean(1, 2, 3, 4, 5)
    print(f"Mean: {mean}")
    
    # 关键字参数
    def process_fastq(file_path, quality_threshold=20, min_length=50):
        print(f"Processing {file_path}")
        print(f"Quality threshold: {quality_threshold}")
        print(f"Min length: {min_length}")
    
    process_fastq("sample.fastq", quality_threshold=30)
    
    # 匿名函数
    sequences = ["ATGC", "GCTA", "TTAA"]
    sorted_seqs = sorted(sequences, key=lambda x: x.count('G'))
    print(sorted_seqs)
    
    # 生物信息学函数示例
    def reverse_complement(sequence):
        complement = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
        return ''.join(complement[base] for base in reversed(sequence))
    
    rc = reverse_complement("ATGCGATCGA")
    print(f"Reverse complement: {rc}")

.. code-block:: python

    # 模块与包示例
    
    # 导入标准库模块
    import os
    import sys
    import json
    import re
    
    # 使用os模块
    current_dir = os.getcwd()
    print(f"Current directory: {current_dir}")
    
    # 列出文件
    files = os.listdir(".")
    for f in files[:5]:
        print(f)
    
    # 使用json模块
    gene_data = {
        "gene_id": "ENSG00000141510",
        "gene_name": "TP53",
        "biotype": "protein_coding"
    }
    
    # 写入JSON文件
    with open("gene_info.json", "w") as f:
        json.dump(gene_data, f, indent=2)
    
    # 读取JSON文件
    with open("gene_info.json", "r") as f:
        loaded_data = json.load(f)
    print(loaded_data)
    
    # 导入特定对象
    from collections import Counter
    
    bases = "ATGCGATCGATCGATCG"
    base_counts = Counter(bases)
    print(base_counts)

.. code-block:: python

    # 文件操作示例
    
    # 读取FASTA文件
    def read_fasta(file_path):
        sequences = {}
        current_id = None
        current_seq = []
        
        with open(file_path, 'r') as f:
            for line in f:
                line = line.strip()
                if line.startswith('>'):
                    if current_id:
                        sequences[current_id] = ''.join(current_seq)
                    current_id = line[1:].split()[0]
                    current_seq = []
                else:
                    current_seq.append(line)
            
            if current_id:
                sequences[current_id] = ''.join(current_seq)
        
        return sequences
    
    # 写入FASTA文件
    def write_fasta(sequences, output_file, line_width=60):
        with open(output_file, 'w') as f:
            for seq_id, sequence in sequences.items():
                f.write(f">{seq_id}\n")
                for i in range(0, len(sequence), line_width):
                    f.write(sequence[i:i+line_width] + "\n")
    
    # 示例使用
    sequences = {
        "gene1": "ATGCGATCGATCGATCGATCG",
        "gene2": "GCTAGCTAGCTAGCTAGCTA"
    }
    write_fasta(sequences, "output.fasta")
    
    # 处理大文件（逐行读取）
    def process_large_file(file_path):
        with open(file_path, 'r') as f:
            for line_num, line in enumerate(f, 1):
                if line_num % 1000000 == 0:
                    print(f"Processed {line_num} lines")
                # 处理每一行
                process_line(line.strip())
    
    def process_line(line):
        pass

.. code-block:: python

    # 异常处理示例
    
    # 基本异常处理
    try:
        with open("nonexistent.txt", 'r') as f:
            content = f.read()
    except FileNotFoundError:
        print("File not found!")
    except PermissionError:
        print("Permission denied!")
    except Exception as e:
        print(f"An error occurred: {e}")
    
    # 带else和finally的异常处理
    try:
        result = 10 / 2
    except ZeroDivisionError:
        print("Cannot divide by zero!")
    else:
        print(f"Result: {result}")
    finally:
        print("This always executes")
    
    # 主动抛出异常
    def validate_sequence(sequence):
        valid_bases = set('ATGC')
        if not all(base in valid_bases for base in sequence.upper()):
            raise ValueError("Invalid sequence: contains non-ATGC characters")
        return True
    
    try:
        validate_sequence("ATGXYZ")
    except ValueError as e:
        print(e)
    
    # 自定义异常
    class SequenceError(Exception):
        pass
    
    def check_sequence_length(sequence, min_length=10):
        if len(sequence) < min_length:
            raise SequenceError(f"Sequence too short: {len(sequence)} < {min_length}")
        return True

.. code-block:: python

    # 面向对象编程示例
    
    # 定义基因类
    class Gene:
        def __init__(self, gene_id, name, chromosome, start, end, strand='+'):
            self.gene_id = gene_id
            self.name = name
            self.chromosome = chromosome
            self.start = start
            self.end = end
            self.strand = strand
            self.transcripts = []
        
        def length(self):
            return self.end - self.start + 1
        
        def add_transcript(self, transcript_id):
            self.transcripts.append(transcript_id)
        
        def __str__(self):
            return f"{self.name} ({self.chromosome}:{self.start}-{self.end})"
        
        def __repr__(self):
            return f"Gene('{self.gene_id}', '{self.name}')"
    
    # 创建基因对象
    tp53 = Gene("ENSG00000141510", "TP53", "chr17", 7668421, 7687490)
    print(tp53)
    print(f"Length: {tp53.length()} bp")
    
    # 继承示例
    class ProteinCodingGene(Gene):
        def __init__(self, gene_id, name, chromosome, start, end, strand='+'):
            super().__init__(gene_id, name, chromosome, start, end, strand)
            self.protein_length = None
        
        def set_protein_length(self, length):
            self.protein_length = length
        
        def get_exon_count(self):
            return len([t for t in self.transcripts if 'exon' in t])
    
    # 使用子类
    brca1 = ProteinCodingGene("ENSG00000012048", "BRCA1", "chr17", 43044295, 43125483)
    brca1.set_protein_length(1863)
    print(f"Protein length: {brca1.protein_length}")

.. code-block:: python

    # 生物信息学常用库示例
    
    # Biopython示例
    from Bio import SeqIO
    from Bio.Seq import Seq
    from Bio.SeqRecord import SeqRecord
    
    # 创建序列对象
    seq = Seq("ATGCGATCGATCGATCG")
    print(f"Sequence: {seq}")
    print(f"Complement: {seq.complement()}")
    print(f"Reverse complement: {seq.reverse_complement()}")
    print(f"Translation: {seq.translate()}")
    
    # 读取FASTA文件
    for record in SeqIO.parse("example.fasta", "fasta"):
        print(f"ID: {record.id}")
        print(f"Length: {len(record.seq)}")
        print(f"GC content: {100 * (record.seq.count('G') + record.seq.count('C')) / len(record.seq):.2f}%")
    
    # pandas示例
    import pandas as pd
    
    # 创建DataFrame
    expression_df = pd.DataFrame({
        'gene': ['TP53', 'BRCA1', 'EGFR', 'MYC'],
        'sample1': [10.5, 5.2, 20.1, 15.3],
        'sample2': [12.3, 6.1, 18.5, 14.8],
        'sample3': [8.7, 4.8, 22.3, 16.1]
    })
    
    print(expression_df)
    
    # 计算每行的平均值
    expression_df['mean'] = expression_df[['sample1', 'sample2', 'sample3']].mean(axis=1)
    
    # 筛选高表达基因
    high_expr = expression_df[expression_df['mean'] > 10]
    print(high_expr)
    
    # NumPy示例
    import numpy as np
    
    # 创建数组
    expression_array = np.array([[10.5, 12.3, 8.7],
                                  [5.2, 6.1, 4.8],
                                  [20.1, 18.5, 22.3]])
    
    # 计算统计量
    print(f"Mean: {np.mean(expression_array)}")
    print(f"Std: {np.std(expression_array)}")
    print(f"Max: {np.max(expression_array)}")
    
    # 按行计算
    row_means = np.mean(expression_array, axis=1)
    print(f"Row means: {row_means}")
    
    # Matplotlib示例
    import matplotlib.pyplot as plt
    
    # 简单的折线图
    genes = ['TP53', 'BRCA1', 'EGFR', 'MYC']
    values = [10.5, 5.2, 20.1, 15.3]
    
    plt.figure(figsize=(10, 6))
    plt.bar(genes, values)
    plt.xlabel('Gene')
    plt.ylabel('Expression Level')
    plt.title('Gene Expression')
    plt.savefig('expression_plot.png', dpi=300)
    plt.close()

正则表达式
----------

正则表达式是一种强大的文本模式匹配工具，它使用特殊的语法来描述字符串的模式。在生物信息学中，正则表达式广泛用于序列模式搜索、数据提取、格式验证等任务。掌握正则表达式可以大大提高文本处理的效率和准确性。

正则表达式基础
~~~~~~~~~~~~~~

正则表达式由普通字符和特殊字符（元字符）组成。普通字符匹配它们自身，如字母a匹配字母a。元字符具有特殊含义，如点号匹配任意单个字符，星号匹配前面的字符零次或多次。理解这些基本元素是使用正则表达式的基础。

常用的元字符包括点号、星号、加号、问号等。点号匹配除换行符外的任意单个字符。星号匹配前面的字符零次或多次。加号匹配前面的字符一次或多次。问号匹配前面的字符零次或一次。这些量词可以组合使用，构建复杂的匹配模式。

字符类用于匹配一组字符中的任意一个。方括号定义字符类，如[abc]匹配a、b或c中的任意一个。可以使用连字符表示范围，如[a-z]匹配任意小写字母，[0-9]匹配任意数字。脱字符在字符类开头表示否定，如[^0-9]匹配任意非数字字符。

预定义字符类提供了常用字符类的简写形式。\\d匹配任意数字，等价于[0-9]。\\D匹配任意非数字字符。\\w匹配任意单词字符，包括字母、数字和下划线。\\W匹配任意非单词字符。\\s匹配任意空白字符，包括空格、制表符、换行符等。\\S匹配任意非空白字符。

锚点用于指定匹配的位置。脱字符匹配字符串的开头，美元符号匹配字符串的结尾。\\b匹配单词边界，\\B匹配非单词边界。锚点不消耗字符，只是指定匹配必须发生的位置。

分组与捕获
~~~~~~~~~~

圆括号用于分组，可以将多个字符作为一个整体处理。分组可以应用量词，如(abc)+匹配一个或多个abc序列。分组还可以用于捕获匹配的内容，捕获的内容可以在后续引用或提取。

捕获组按照左括号的顺序编号，第一个捕获组是组1，第二个是组2，依此类推。组0表示整个匹配。在正则表达式内部，可以使用反斜杠加数字引用之前的捕获组，如\\1引用第一个捕获组。这在匹配重复内容时很有用，如匹配引号内的字符串，(["'])\\1匹配成对的引号。

非捕获组使用(?:...)语法，它只分组不捕获。非捕获组在不需要提取匹配内容时使用，可以提高性能。命名捕获组使用(?P<name>...)语法，可以给捕获组命名，便于后续引用。

Python的re模块
~~~~~~~~~~~~~~

Python的re模块提供了正则表达式操作的功能。re模块的主要函数包括match、search、findall、finditer、sub、split等。match函数从字符串开头开始匹配，search函数在整个字符串中搜索第一个匹配，findall函数返回所有匹配的列表，finditer函数返回所有匹配的迭代器。

正则表达式可以使用原始字符串表示，在字符串前加r前缀。原始字符串中的反斜杠不需要转义，使正则表达式更加清晰。例如，r'\\d+'比'\\\\d+'更易读。建议在编写正则表达式时始终使用原始字符串。

编译正则表达式使用re.compile函数，返回一个正则表达式对象。编译后的正则表达式可以重复使用，提高性能。对于需要多次使用的正则表达式，建议先编译。正则表达式对象的方法与re模块的函数类似，但不需要传入正则表达式参数。

匹配对象是match和search函数返回的结果，包含了匹配的详细信息。匹配对象提供了group方法获取捕获组的内容，start和end方法获取匹配的位置，span方法返回匹配位置的元组。使用这些方法可以精确地处理匹配结果。

生物信息学应用
~~~~~~~~~~~~~~

正则表达式在生物信息学中有广泛的应用。在序列分析中，可以用来搜索特定的序列模式，如限制性酶切位点、转录因子结合位点、蛋白质结构域等。在数据提取中，可以从文本文件中提取基因ID、序列、注释等信息。在数据清洗中，可以验证和标准化数据格式。

处理FASTA和FASTQ文件时，正则表达式可以提取序列标识符、序列内容、质量值等信息。例如，从FASTA头部提取基因ID和描述信息，从FASTQ文件中提取序列和质量值。正则表达式还可以验证序列格式，检查是否包含非法字符。

处理基因组注释文件时，正则表达式可以解析GFF、GTF、BED等格式的文件。可以提取基因的位置信息、属性字段中的特定值。例如，从GFF文件的属性字段中提取gene_id、transcript_id等信息。

处理变异数据时，正则表达式可以解析VCF文件，提取变异位置、参考碱基、变异碱基等信息。可以验证变异格式的正确性，提取复杂的变异信息。正则表达式在处理复杂的变异表示时特别有用。

.. code-block:: python

    # 正则表达式基础示例
    import re
    
    # 基本匹配
    text = "Gene TP53 is located on chromosome 17"
    
    # search: 搜索第一个匹配
    match = re.search(r'TP53', text)
    if match:
        print(f"Found: {match.group()}")
        print(f"Position: {match.start()}-{match.end()}")
    
    # findall: 找到所有匹配
    numbers = re.findall(r'\d+', text)
    print(f"Numbers: {numbers}")
    
    # 元字符示例
    sequence = "ATGCGATCGATCGATCG"
    
    # 点号匹配任意字符
    print(re.findall(r'A.G', sequence))
    
    # 星号匹配零次或多次
    print(re.findall(r'AT*', sequence))
    
    # 加号匹配一次或多次
    print(re.findall(r'AT+', sequence))
    
    # 问号匹配零次或一次
    print(re.findall(r'AT?', sequence))
    
    # 字符类示例
    # 匹配任意碱基
    print(re.findall(r'[ATGC]', sequence))
    
    # 匹配非数字
    text2 = "chr17:7668421-7687490"
    print(re.findall(r'[^0-9]', text2))
    
    # 预定义字符类
    print(re.findall(r'\d+', text2))
    print(re.findall(r'\D+', text2))

.. code-block:: python

    # 分组与捕获示例
    
    # 捕获组
    gff_line = "chr17\tHAVANA\tgene\t7668421\t7687490\t.\t+\t.\tID=gene:TP53;Name=TP53"
    
    # 提取基因位置信息
    pattern = r'(\w+)\t(\w+)\t(\w+)\t(\d+)\t(\d+)'
    match = re.search(pattern, gff_line)
    if match:
        print(f"Chromosome: {match.group(1)}")
        print(f"Feature: {match.group(3)}")
        print(f"Start: {match.group(4)}")
        print(f"End: {match.group(5)}")
    
    # 提取属性中的ID和Name
    attr_pattern = r'ID=gene:(\w+);Name=(\w+)'
    match = re.search(attr_pattern, gff_line)
    if match:
        print(f"Gene ID: {match.group(1)}")
        print(f"Gene Name: {match.group(2)}")
    
    # 命名捕获组
    pattern = r'(?P<chr>\w+):(?P<start>\d+)-(?P<end>\d+)'
    text = "chr17:7668421-7687490"
    match = re.search(pattern, text)
    if match:
        print(f"Chromosome: {match.group('chr')}")
        print(f"Start: {match.group('start')}")
        print(f"End: {match.group('end')}")
    
    # 非捕获组
    text = "sample1_rep1 sample2_rep1 sample1_rep2"
    pattern = r'(?:sample\d+)_rep\d+'
    matches = re.findall(pattern, text)
    print(matches)

.. code-block:: python

    # Python re模块详解
    
    # compile: 编译正则表达式
    pattern = re.compile(r'\d+')
    text = "chr17:7668421-7687490"
    print(pattern.findall(text))
    
    # match: 从开头匹配
    text = "TP53 BRCA1 EGFR"
    match = re.match(r'\w+', text)
    if match:
        print(f"Matched: {match.group()}")
    
    # search: 搜索第一个匹配
    match = re.search(r'BRCA\d', text)
    if match:
        print(f"Found: {match.group()}")
    
    # findall: 返回所有匹配的列表
    genes = re.findall(r'[A-Z]+\d*', text)
    print(genes)
    
    # finditer: 返回匹配对象的迭代器
    for match in re.finditer(r'[A-Z]+\d*', text):
        print(f"{match.group()} at position {match.start()}")
    
    # sub: 替换
    text = "chr1 chr2 chrX chrY"
    result = re.sub(r'chr', 'chromosome ', text)
    print(result)
    
    # split: 分割
    text = "TP53,BRCA1;EGFR:MYC"
    genes = re.split(r'[,;:]', text)
    print(genes)

.. code-block:: python

    # 生物信息学应用示例
    
    # 处理FASTA文件
    fasta_header = ">ENSG00000141510 TP53 tumor protein p53 [Homo sapiens]"
    
    # 提取基因ID和名称
    pattern = r'>(\w+)\s+(\w+)\s+(.+)'
    match = re.match(pattern, fasta_header)
    if match:
        gene_id = match.group(1)
        gene_name = match.group(2)
        description = match.group(3)
        print(f"ID: {gene_id}, Name: {gene_name}, Desc: {description}")
    
    # 验证DNA序列
    def validate_dna(sequence):
        if re.match(r'^[ATGC]+$', sequence):
            return True
        return False
    
    print(validate_dna("ATGCGATCGA"))
    print(validate_dna("ATGCGATXGA"))
    
    # 搜索限制性酶切位点
    sequence = "ATGCGAATTCGATCGAATTCCGATCG"
    ecoRI_pattern = r'GAATTC'
    for match in re.finditer(ecoRI_pattern, sequence):
        print(f"EcoRI site at position {match.start()}")
    
    # 提取VCF文件信息
    vcf_line = "17\t7668421\t.\tG\tA\t100\tPASS\tAC=1;AF=0.5;AN=2"
    pattern = r'(\d+)\t(\d+)\t(\S+)\t(\w)\t(\w)\t(\d+)\t(\w+)\t(.+)'
    match = re.match(pattern, vcf_line)
    if match:
        chrom = match.group(1)
        pos = match.group(2)
        ref = match.group(4)
        alt = match.group(5)
        print(f"chr{chrom}:{pos} {ref}>{alt}")
    
    # 解析GFF属性字段
    def parse_gff_attributes(attr_string):
        attrs = {}
        for item in attr_string.split(';'):
            if '=' in item:
                key, value = item.split('=', 1)
                attrs[key] = value
        return attrs
    
    attrs = parse_gff_attributes("ID=gene:TP53;Name=TP53;biotype=protein_coding")
    print(attrs)
    
    # 提取基因ID格式
    text = "ENSG00000141510 ENSG00000012048 ENSMUSG00000059552"
    human_genes = re.findall(r'ENSG\d+', text)
    mouse_genes = re.findall(r'ENSMUSG\d+', text)
    print(f"Human: {human_genes}")
    print(f"Mouse: {mouse_genes}")
    
    # 处理FASTQ质量值
    fastq_record = """@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65"""
    
    lines = fastq_record.strip().split('\n')
    seq_id = re.match(r'@(\S+)', lines[0]).group(1)
    sequence = lines[1]
    quality = lines[3]
    print(f"ID: {seq_id}")
    print(f"Sequence length: {len(sequence)}")
    print(f"Quality string length: {len(quality)}")

版本控制
--------

版本控制是管理文件变更历史的系统，它记录文件的每一次修改，允许回溯到之前的版本，协调多人协作开发。在生物信息学项目中，版本控制是保证分析可重复性和代码可维护性的重要工具。Git是目前最流行的分布式版本控制系统，广泛应用于软件开发和科学研究领域。

Git基础概念
~~~~~~~~~~~

Git是一个分布式版本控制系统，每个开发者都有完整的代码库副本，包括完整的历史记录。这种设计使得Git在没有网络连接时也能工作，并且比集中式系统更安全。Git的核心概念包括仓库、提交、分支、远程仓库等。

仓库是Git管理的基本单位，包含项目的所有文件和完整的历史记录。仓库可以是本地的，也可以托管在远程服务器上。初始化仓库使用git init命令，会在当前目录创建.git目录，存储所有版本控制信息。克隆远程仓库使用git clone命令，会复制整个仓库到本地。

提交是Git的基本操作单元，每次提交记录了文件的一个快照。提交包含唯一的SHA-1哈希值、作者信息、时间戳和提交消息。提交消息应该简洁明了地描述本次修改的内容。良好的提交习惯是版本控制的基础，应该频繁提交，每次提交只包含逻辑相关的修改。

工作区、暂存区和仓库是Git的三个重要区域。工作区是实际的工作目录，包含项目的文件。暂存区是准备提交的文件列表，使用git add命令将工作区的修改添加到暂存区。仓库包含所有提交的历史记录，使用git commit命令将暂存区的内容提交到仓库。

分支是指向某个提交的可移动指针，允许在不同的开发线上并行工作。默认分支通常叫master或main。创建分支使用git branch命令，切换分支使用git checkout或git switch命令。分支使得开发新功能、修复bug等任务可以独立进行，不影响主分支。

基本操作流程
~~~~~~~~~~~~

Git的基本工作流程包括修改文件、暂存修改、提交修改。首先在工作区修改文件，然后使用git add命令将修改添加到暂存区，最后使用git commit命令将暂存区的内容提交到仓库。这个流程可以重复进行，每次提交都会创建一个新的版本。

查看状态使用git status命令，它会显示工作区和暂存区的状态，包括哪些文件被修改、哪些文件已暂存、哪些文件未被跟踪。查看提交历史使用git log命令，它会显示提交的哈希值、作者、日期和提交消息。

比较差异使用git diff命令。不带参数的git diff显示工作区和暂存区的差异。git diff --staged显示暂存区和最新提交的差异。git diff commit1 commit2显示两个提交之间的差异。这些命令帮助理解修改的内容。

撤销修改有几种方式。git checkout -- file丢弃工作区的修改。git reset HEAD file取消暂存。git reset --hard commit回退到指定的提交。注意，这些操作可能丢失修改，应该谨慎使用。在执行撤销操作前，确保理解其影响。

远程仓库操作
~~~~~~~~~~~~

远程仓库是托管在网络服务器上的仓库，用于多人协作和备份。常见的远程仓库托管服务包括GitHub、GitLab、Bitbucket等。添加远程仓库使用git remote add命令，查看远程仓库使用git remote -v命令。

推送使用git push命令，将本地的提交推送到远程仓库。首次推送需要使用-u参数设置上游分支，git push -u origin master。之后可以直接使用git push。推送前应该先拉取远程的更新，避免冲突。

拉取使用git pull命令，从远程仓库获取最新的提交并合并到当前分支。git fetch命令只获取远程的更新，不自动合并。建议先fetch查看更新内容，再决定是否合并。拉取时可能遇到冲突，需要手动解决。

克隆使用git clone命令，复制远程仓库到本地。克隆会自动设置远程仓库为origin，并创建默认分支的跟踪分支。克隆后可以正常进行修改、提交、推送等操作。克隆是参与已有项目的第一步。

分支管理
~~~~~~~~

分支是Git的核心功能，它允许在不同的开发线上并行工作。创建分支使用git branch branchname命令。切换分支使用git checkout branchname或git switch branchname命令。创建并切换分支可以使用git checkout -b branchname命令。

合并分支使用git merge命令。首先切换到目标分支，然后执行git merge source-branch。合并可能产生冲突，需要手动解决。冲突的文件会包含冲突标记，编辑文件解决冲突后，使用git add标记为已解决，然后提交。

删除分支使用git branch -d branchname命令。已合并的分支可以安全删除。如果分支未合并，需要使用-D参数强制删除。删除分支不会删除提交，只是删除分支指针。

分支策略是团队协作的重要规范。常见的策略包括Git Flow、GitHub Flow等。Git Flow使用master和develop两个长期分支，以及feature、release、hotfix等短期分支。GitHub Flow更简单，只使用master分支和feature分支。选择合适的分支策略可以提高协作效率。

协作工作流程
~~~~~~~~~~~~

多人协作时，通常使用fork和pull request模式。首先fork项目到自己的账户，然后克隆fork的仓库到本地。创建分支进行开发，完成后推送到自己的仓库。然后在GitHub上创建pull request，请求将修改合并到原项目。

代码审查是协作的重要环节。pull request提供了讨论和审查的平台。审查者可以评论代码，提出修改建议。作者可以根据反馈修改代码，新的提交会自动添加到pull request。审查通过后，维护者合并pull request。

保持同步是协作的常见需求。当原项目有新的提交时，需要将这些提交同步到自己的fork。可以添加原项目为远程仓库，然后fetch并merge。或者使用GitHub的sync fork功能。定期同步可以减少合并冲突。

解决冲突是协作中不可避免的。当多人修改同一文件的同一部分时，会产生冲突。Git会标记冲突的位置，需要手动编辑文件解决。解决冲突后，标记为已解决并提交。良好的沟通和协调可以减少冲突的发生。

Git在生物信息学中的应用
~~~~~~~~~~~~~~~~~~~~~~~~

在生物信息学项目中，Git用于管理分析脚本、配置文件、文档等。分析脚本应该纳入版本控制，记录每次修改的原因和内容。配置文件记录分析参数的变更。文档记录项目的进展和结果。

数据文件通常不纳入Git管理，因为它们通常很大。可以使用.gitignore文件排除数据文件。大型数据文件可以使用Git LFS或专门的数据管理系统。但小型的示例数据可以纳入版本控制，用于测试和演示。

分析流程的可重复性是科学研究的基本要求。使用Git可以记录分析流程的每一个版本，包括使用的脚本、参数和配置。当需要重复分析时，可以检出特定版本的代码。这保证了分析的可重复性。

协作开发生物信息学工具时，Git提供了完整的协作平台。可以创建组织来管理项目，使用issue跟踪bug和功能请求，使用wiki编写文档，使用GitHub Pages发布网站。这些功能使得Git成为生物信息学开发的中心平台。

.. code-block:: bash

    # Git基础操作示例
    
    # 初始化仓库
    mkdir bioinfo_project
    cd bioinfo_project
    git init
    
    # 配置用户信息
    git config --global user.name "Your Name"
    git config --global user.email "your.email@example.com"
    
    # 添加文件到暂存区
    git add analysis.py
    git add README.md
    
    # 添加所有修改的文件
    git add .
    
    # 提交修改
    git commit -m "Initial commit: add analysis script and README"
    
    # 查看状态
    git status
    
    # 查看提交历史
    git log
    git log --oneline
    git log --graph --oneline --all
    
    # 查看差异
    git diff
    git diff --staged
    git diff HEAD~1 HEAD

.. code-block:: bash

    # 分支操作示例
    
    # 创建新分支
    git branch feature-add-qc
    
    # 切换分支
    git checkout feature-add-qc
    # 或者使用新命令
    git switch feature-add-qc
    
    # 创建并切换分支
    git checkout -b feature-add-qc
    
    # 查看所有分支
    git branch -a
    
    # 在新分支上工作
    echo "def quality_control():" >> analysis.py
    git add analysis.py
    git commit -m "Add quality control function"
    
    # 切换回主分支
    git checkout master
    
    # 合并分支
    git merge feature-add-qc
    
    # 删除已合并的分支
    git branch -d feature-add-qc
    
    # 强制删除未合并的分支
    git branch -D feature-abandoned

.. code-block:: bash

    # 远程仓库操作示例
    
    # 克隆远程仓库
    git clone https://github.com/user/bioinfo_project.git
    cd bioinfo_project
    
    # 查看远程仓库
    git remote -v
    
    # 添加远程仓库
    git remote add origin https://github.com/user/bioinfo_project.git
    
    # 推送到远程仓库
    git push -u origin master
    git push
    
    # 拉取远程更新
    git pull origin master
    
    # 获取远程更新但不合并
    git fetch origin
    
    # 查看远程分支
    git branch -r
    
    # 创建本地分支跟踪远程分支
    git checkout -b feature-branch origin/feature-branch

.. code-block:: bash

    # 撤销操作示例
    
    # 丢弃工作区的修改
    git checkout -- analysis.py
    
    # 取消暂存
    git reset HEAD analysis.py
    
    # 修改最后一次提交
    git commit --amend -m "New commit message"
    
    # 回退到上一个版本
    git reset --hard HEAD^
    
    # 回退到指定版本
    git reset --hard a1b2c3d
    
    # 查看操作历史（用于恢复误删的提交）
    git reflog
    
    # 创建标签
    git tag v1.0.0
    git tag -a v1.0.0 -m "Release version 1.0.0"
    
    # 推送标签
    git push origin v1.0.0
    git push origin --tags

.. code-block:: bash

    # 生物信息学项目示例
    
    # 创建项目目录结构
    mkdir -p bioinfo_analysis/{scripts,data,results,docs}
    cd bioinfo_analysis
    
    # 初始化Git
    git init
    
    # 创建.gitignore文件
    cat > .gitignore << 'EOF'
    # 数据文件
    *.fastq
    *.fasta
    *.bam
    *.vcf
    
    # 结果文件
    results/*.txt
    results/*.pdf
    
    # 临时文件
    *.log
    *.tmp
    __pycache__/
    
    # 系统文件
    .DS_Store
    Thumbs.db
    EOF
    
    # 添加README
    cat > README.md << 'EOF'
    # Bioinformatics Analysis Project
    
    ## Description
    RNA-seq analysis pipeline
    
    ## Requirements
    - Python 3.9+
    - STAR
    - featureCounts
    
    ## Usage
    bash scripts/run_pipeline.sh
    
    ## Author
    Your Name
    EOF
    
    # 提交初始文件
    git add .gitignore README.md
    git commit -m "Initial commit: add project structure and README"
    
    # 添加分析脚本
    cat > scripts/quality_control.py << 'EOF'
    #!/usr/bin/env python3
    """Quality control for sequencing data"""
    
    import sys
    import argparse
    
    def main():
        parser = argparse.ArgumentParser(description='Quality control')
        parser.add_argument('-i', '--input', required=True, help='Input FASTQ file')
        parser.add_argument('-o', '--output', required=True, help='Output report')
        args = parser.parse_args()
        
        print(f"Processing {args.input}")
        # QC code here
        
    if __name__ == '__main__':
        main()
    EOF
    
    git add scripts/
    git commit -m "Add quality control script"
    
    # 创建开发分支
    git checkout -b develop
    
    # 推送到远程
    git remote add origin https://github.com/user/bioinfo_analysis.git
    git push -u origin master
    git push origin develop

容器化技术
----------

容器化技术是一种轻量级的虚拟化技术，它将应用程序及其所有依赖打包到一个独立的容器中运行。与传统的虚拟机相比，容器更加轻量、启动更快、资源占用更少。在生物信息学中，容器化技术解决了软件依赖和环境配置的问题，确保分析流程在不同计算环境中的一致性和可重复性。

Docker基础
~~~~~~~~~~

Docker是最流行的容器化平台，它使用客户端-服务器架构。Docker客户端与Docker守护进程通信，守护进程负责构建、运行和分发容器。Docker镜像是一个只读模板，包含创建容器的指令。Docker容器是镜像的运行实例，可以被启动、停止、移动和删除。

镜像使用分层存储，每个层代表Dockerfile中的一条指令。这种设计使得镜像构建高效，共享层可以节省存储空间。镜像可以从Docker Hub等仓库拉取，也可以自己构建。构建镜像使用Dockerfile，它是一个文本文件，包含构建镜像的所有指令。

Dockerfile是构建镜像的蓝图，从基础镜像开始，逐步添加软件和配置。FROM指令指定基础镜像，RUN指令执行命令，COPY和ADD指令添加文件，WORKDIR设置工作目录，EXPOSE声明端口，CMD和ENTRYPOINT指定容器启动命令。编写Dockerfile需要考虑镜像大小、构建速度和安全性。

容器是镜像的运行实例。使用docker run命令创建和启动容器。容器有独立的文件系统、网络配置和进程空间。容器可以交互式运行，也可以在后台运行。容器的数据默认存储在容器内部，容器删除后数据会丢失，可以使用数据卷持久化数据。

Docker Compose
~~~~~~~~~~~~~~

Docker Compose是一个用于定义和运行多容器应用的工具。使用YAML文件配置应用的服务、网络和卷。Docker Compose适合开发、测试和CI/CD环境，可以一键启动复杂的多服务应用。

Compose文件定义了服务、网络和卷。服务是容器的配置，包括镜像、端口映射、环境变量、依赖关系等。网络定义了服务之间的通信方式。卷定义了数据持久化的方式。Compose文件版本指定了文件格式的版本。

常用命令包括docker-compose up启动服务，docker-compose down停止服务，docker-compose ps查看服务状态，docker-compose logs查看日志。使用-d参数在后台运行服务。使用--build参数在启动前重新构建镜像。

在生物信息学中，Docker Compose可以用于搭建复杂的数据分析平台。例如，可以配置数据库服务、Web服务、分析服务等，一键启动整个分析环境。这大大简化了环境部署的复杂度。

Singularity容器
~~~~~~~~~~~~~~~

Singularity是另一种容器技术，特别适合在高性能计算环境中使用。与Docker需要root权限不同，Singularity允许普通用户运行容器，这使得它在共享的计算集群中更受欢迎。

Singularity使用SIF格式的镜像文件，这是一个单独的文件，便于传输和存档。Singularity可以从Docker Hub拉取镜像并转换为SIF格式，也可以使用Singularity定义文件构建镜像。Singularity容器默认只读，可以使用overlay或绑定挂载实现数据写入。

Singularity的命令与Docker类似。singularity pull拉取镜像，singularity build构建镜像，singularity run运行容器，singularity shell进入容器交互模式。Singularity还支持执行容器中的特定命令。

在生物信息学集群环境中，Singularity是运行容器化分析流程的首选。许多生物信息学工具都提供了Singularity镜像，可以直接使用。Singularity与工作流管理系统（如Nextflow、Snakemake）集成良好。

容器化生物信息学流程
~~~~~~~~~~~~~~~~~~~~

容器化生物信息学流程可以确保分析的可重复性。将分析流程的每个步骤打包成容器，可以保证在不同环境中运行结果一致。工作流管理系统可以协调多个容器的执行，管理数据流和依赖关系。

常用的生物信息学容器镜像可以在Biocontainers项目找到。Biocontainers为Bioconda中的软件提供了对应的容器镜像，覆盖了数千个生物信息学工具。使用这些镜像可以避免复杂的软件安装过程。

构建自定义镜像时，应该选择合适的基础镜像。对于生物信息学应用，通常选择包含常用工具的基础镜像。Dockerfile应该简洁明了，每条指令的目的清晰。应该清理不必要的文件，减小镜像大小。

数据管理是容器化流程的重要考虑。容器内部的数据在容器删除后会丢失。应该使用数据卷或绑定挂载将数据存储在容器外部。在集群环境中，通常将数据存储在共享文件系统中，通过绑定挂载访问。

集群作业调度
------------

集群作业调度系统管理计算集群的资源分配和任务执行。在高性能计算环境中，用户提交作业到调度系统，调度系统根据资源可用性和调度策略分配计算节点执行作业。掌握集群作业调度是进行大规模生物信息学分析的必备技能。

作业调度系统概述
~~~~~~~~~~~~~~~~

作业调度系统是高性能计算集群的核心组件，它负责接收用户提交的作业、分配计算资源、监控作业执行、管理作业队列。调度系统确保集群资源的公平共享和高效利用，防止单个用户独占资源。

常见的作业调度系统包括Slurm、PBS、SGE等。Slurm是目前最流行的开源调度系统，被许多超算中心和科研机构采用。PBS是传统的调度系统，有OpenPBS和PBS Pro两个版本。SGE是Sun Grid Engine的简称，现在有Open Grid Scheduler等开源版本。

作业调度涉及几个核心概念。作业是用户提交的计算任务，包含资源需求、执行命令等信息。队列是作业的等待列表，作业在队列中等待调度。节点是计算资源的单位，可以是物理服务器或虚拟机。核心是CPU的计算单元，是资源分配的基本单位。

资源请求是作业描述的重要组成部分。用户需要指定作业需要的CPU核心数、内存大小、运行时间等。调度系统根据资源请求和集群资源状况决定何时何地运行作业。合理的资源请求可以提高作业的调度效率和成功率。

Slurm调度系统
~~~~~~~~~~~~~

Slurm是Simple Linux Utility for Resource Management的缩写，是一个开源的、高度可扩展的集群管理和作业调度系统。Slurm具有高可用性、可扩展性和可移植性，支持多种调度策略和资源管理方式。

Slurm的基本命令包括sbatch提交批处理作业、srun提交交互式作业、squeue查看作业队列、scancel取消作业、sinfo查看节点信息、sacct查看作业历史。这些命令是日常使用Slurm的基础。

sbatch命令用于提交批处理作业脚本。脚本中包含Slurm指令，以#SBATCH开头，指定作业名称、资源需求、输出文件等。脚本主体是要执行的命令。sbatch命令返回作业ID，用户可以通过作业ID查询和管理作业。

srun命令用于提交交互式作业或并行作业。srun会在分配的节点上执行命令，命令的输出直接显示在终端。srun常用于测试和调试，也可以在脚本中使用srun启动并行任务。

PBS调度系统
~~~~~~~~~~~

PBS是Portable Batch System的缩写，是一个历史悠久的作业调度系统。PBS使用qsub提交作业、qstat查看作业状态、qdel取消作业。PBS的作业脚本使用#PBS指令指定资源需求。

PBS的作业脚本与Slurm类似，但指令语法略有不同。PBS使用-l参数指定资源，如-l nodes=1:ppn=8指定1个节点8个核心。PBS使用-N指定作业名称，-o和-e指定输出和错误文件。

PBS支持作业数组，可以一次提交多个相似的作业。作业数组适合处理大量相似的任务，如处理多个样本文件。使用-t参数指定数组索引范围，PBS会为数组中的每个索引创建一个子作业。

PBS的作业状态包括Q（排队）、R（运行）、C（完成）、E（退出）等。使用qstat命令查看作业状态，可以指定作业ID或用户名过滤输出。PBS还提供pbsnodes命令查看节点状态。

作业脚本编写
~~~~~~~~~~~~

作业脚本是提交到调度系统的文本文件，包含资源请求和执行命令。脚本通常以shell脚本格式编写，开头是调度系统的指令，后面是要执行的命令。编写好的作业脚本可以提高作业的成功率和效率。

资源请求应该根据实际需求合理设置。CPU核心数应该与程序的并行度匹配，过多的核心请求会延长等待时间。内存请求应该留有余量，避免因内存不足导致作业失败。运行时间应该根据经验估计，过短会导致作业被终止，过长会延长等待时间。

作业脚本应该包含错误处理和日志记录。检查输入文件是否存在，检查上一阶段是否成功完成。将标准输出和错误输出重定向到文件，便于调试。使用作业名称区分不同的作业。

生物信息学作业通常需要处理大量数据。应该使用共享文件系统存储数据，避免数据复制。使用环境变量或参数传递文件路径，使脚本更加灵活。对于多步骤分析，可以将每个步骤写成独立的作业，使用作业依赖关系串联。

生物信息学集群应用
~~~~~~~~~~~~~~~~~~

在生物信息学中，集群计算用于处理大规模数据分析任务。常见的应用包括基因组比对、变异检测、转录组定量、宏基因组分析等。这些任务通常需要大量的计算资源和内存。

基因组比对是计算密集型任务，可以使用多线程加速。BWA、STAR等比对工具支持多线程运行。作业脚本应该请求足够的CPU核心和内存。对于大型基因组，可能需要请求大内存节点。

转录组定量可以使用并行处理加速。对于多样本研究，可以将样本分配到不同的计算节点并行处理。使用作业数组可以简化多样本作业的提交。定量结果汇总可以在单独的作业中完成。

变异数据处理包括变异检测、注释、过滤等步骤。GATK等工具支持并行处理。对于全基因组数据，可以按染色体分割任务并行处理。作业依赖关系可以确保处理步骤的正确顺序。

工作流管理系统如Nextflow和Snakemake可以简化集群作业的管理。这些系统自动处理作业提交、依赖管理、错误恢复等。工作流定义文件描述分析步骤和数据流，系统负责执行和监控。

数据格式与处理
--------------

生物信息学数据具有特定的格式和结构，理解这些格式是进行数据分析的基础。常见的数据格式包括序列格式、比对格式、变异格式、注释格式等。掌握这些格式的结构和处理方法是生物信息学的基本功。

序列文件格式
~~~~~~~~~~~~

FASTA格式是最基本的序列格式，用于存储DNA、RNA或蛋白质序列。FASTA文件由多个记录组成，每个记录以大于号开头，后跟序列标识符和描述信息，接下来是序列数据。序列可以跨越多行，通常每行60或80个字符。

FASTQ格式是存储测序数据的标准格式，包含序列和质量信息。每个FASTQ记录包含四行：第一行以@开头，包含序列标识符；第二行是序列；第三行以+开头，可以重复标识符；第四行是质量字符串，每个字符对应序列中一个碱基的质量值。

质量值使用Phred质量评分表示，衡量碱基测序的准确性。Phred质量分数Q表示错误概率P，Q = -10 * log10(P)。质量分数20表示错误概率1%，质量分数30表示错误概率0.1%。质量字符串使用ASCII字符编码，不同的编码方案有不同的偏移值。

GenBank格式是GenBank数据库使用的序列格式，包含序列、注释和文献信息。GenBank记录以LOCUS行开始，包含序列长度、分子类型等基本信息。FEATURES部分包含基因、CDS等特征的注释。ORIGIN部分是序列数据。

比对文件格式
~~~~~~~~~~~~

SAM格式是序列比对的标准格式，存储测序reads与参考基因组的比对信息。SAM文件由头部和比对记录组成。头部行以@开头，包含参考序列信息、程序信息等。比对记录包含11个必需字段和可选字段。

BAM格式是SAM的二进制版本，占用空间更小，读取速度更快。BAM文件通常需要排序和建立索引，以便快速访问。samtools是处理SAM/BAM文件的主要工具，可以查看、排序、索引、过滤SAM/BAM文件。

比对记录的必需字段包括查询名称、标志位、参考序列名、比对位置、比对质量、CIGAR字符串、配对参考序列名、配对位置、模板长度、序列和质量值。标志位使用位编码表示比对的属性，如是否配对、是否比对上等。

CIGAR字符串描述比对的结构，使用操作符表示匹配、插入、缺失等。M表示匹配或错配，I表示插入，D表示缺失，N表示跳过，S表示软裁剪，H表示硬裁剪。例如，100M表示100个匹配，50M5I45M表示50个匹配、5个插入、45个匹配。

变异数据格式
~~~~~~~~~~~~

VCF格式是存储基因变异信息的标准格式，包括SNP、插入缺失、结构变异等。VCF文件由头部和变异记录组成。头部以##开头，包含格式定义、信息字段描述等。列标题行以#CHROM开头，列出各列的名称。

变异记录包含8个必需字段：染色体、位置、变异ID、参考碱基、变异碱基、质量分数、过滤标志、信息字段。对于多样本数据，还有每个样本的基因型字段。位置是变异的第一个碱基的位置，使用1-based坐标系统。

BCF格式是VCF的二进制版本，与BAM类似，占用空间更小，处理速度更快。BCF文件也需要索引才能快速访问。bcftools是处理VCF/BCF文件的主要工具，可以查看、过滤、注释、合并变异数据。

变异注释是变异数据处理的重要步骤。注释添加变异的功能影响、频率信息、临床意义等。常用的注释工具包括ANNOVAR、SnpEff、VEP等。注释结果可以添加到VCF的INFO字段或单独的文件中。

注释文件格式
~~~~~~~~~~~~

GFF格式是基因组注释的标准格式，描述基因、转录本、外显子等特征的位置和属性。GFF文件每行一个特征，包含9个字段：序列ID、来源、类型、起始位置、终止位置、分数、链、相位、属性。GFF有多个版本，GFF2、GFF3是常用的版本。

GTF格式是GFF的一个变体，主要用于基因注释。GTF与GFF3的主要区别在于属性字段的格式。GTF的属性字段使用键值对，以分号分隔。GFF3的属性字段使用键=值格式，也以分号分隔。

BED格式是简单的区间格式，用于描述基因组区域。BED文件每行一个区间，至少包含染色体、起始位置、终止位置三个字段。BED使用0-based坐标系统，起始位置是第一个碱基的索引，终止位置是最后一个碱基的索引加1。

BED格式支持扩展字段，可以添加名称、分数、链等信息。BED12格式可以描述区块特征，如外显子。BED文件可以用于定义分析区域、可视化注释等。UCSC基因组浏览器和其他工具广泛支持BED格式。

.. code-block:: docker

    # Dockerfile示例：生物信息学分析环境
    
    # 使用官方Python镜像作为基础
    FROM python:3.9-slim
    
    # 设置维护者信息
    LABEL maintainer="your.email@example.com"
    
    # 设置工作目录
    WORKDIR /app
    
    # 安装系统依赖
    RUN apt-get update && apt-get install -y \
        wget \
        bzip2 \
        libncurses5-dev \
        zlib1g-dev \
        && rm -rf /var/lib/apt/lists/*
    
    # 安装Miniconda
    RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh \
        && bash Miniconda3-latest-Linux-x86_64.sh -b -p /opt/conda \
        && rm Miniconda3-latest-Linux-x86_64.sh
    
    # 添加conda到PATH
    ENV PATH="/opt/conda/bin:${PATH}"
    
    # 配置conda频道
    RUN conda config --add channels defaults \
        && conda config --add channels bioconda \
        && conda config --add channels conda-forge
    
    # 安装生物信息学工具
    RUN conda install -y \
        fastqc \
        bwa \
        samtools \
        bcftools \
        && conda clean -afy
    
    # 安装Python包
    COPY requirements.txt .
    RUN pip install --no-cache-dir -r requirements.txt
    
    # 复制分析脚本
    COPY scripts/ /app/scripts/
    
    # 设置入口点
    ENTRYPOINT ["python", "/app/scripts/pipeline.py"]

.. code-block:: bash

    # Docker基本操作示例
    
    # 构建镜像
    docker build -t bioinfo-analysis:1.0 .
    
    # 查看本地镜像
    docker images
    
    # 运行容器
    docker run -it bioinfo-analysis:1.0
    
    # 挂载数据目录运行
    docker run -v /path/to/data:/data bioinfo-analysis:1.0
    
    # 后台运行容器
    docker run -d --name analysis bioinfo-analysis:1.0
    
    # 查看运行中的容器
    docker ps
    
    # 进入运行中的容器
    docker exec -it analysis bash
    
    # 查看容器日志
    docker logs analysis
    
    # 停止容器
    docker stop analysis
    
    # 删除容器
    docker rm analysis
    
    # 删除镜像
    docker rmi bioinfo-analysis:1.0
    
    # 从Docker Hub拉取镜像
    docker pull biocontainers/fastqc:v0.11.9
    
    # 推送镜像到Docker Hub
    docker login
    docker tag bioinfo-analysis:1.0 username/bioinfo-analysis:1.0
    docker push username/bioinfo-analysis:1.0

.. code-block:: yaml

    # Docker Compose示例：多服务生物信息学平台
    
    version: '3.8'
    
    services:
      # 数据库服务
      postgres:
        image: postgres:13
        environment:
          POSTGRES_DB: bioinfo_db
          POSTGRES_USER: bioinfo
          POSTGRES_PASSWORD: password
        volumes:
          - postgres_data:/var/lib/postgresql/data
        ports:
          - "5432:5432"
    
      # Web服务
      web:
        build: ./web
        ports:
          - "8000:8000"
        depends_on:
          - postgres
        environment:
          DATABASE_URL: postgresql://bioinfo:password@postgres:5432/bioinfo_db
    
      # 分析服务
      analysis:
        build: ./analysis
        volumes:
          - ./data:/data
          - ./results:/results
        depends_on:
          - postgres
    
      # Jupyter服务
      jupyter:
        image: jupyter/datascience-notebook
        ports:
          - "8888:8888"
        volumes:
          - ./notebooks:/home/jovyan/work
          - ./data:/home/jovyan/data
    
    volumes:
      postgres_data:

.. code-block:: bash

    # Singularity操作示例
    
    # 从Docker Hub拉取镜像并转换为SIF格式
    singularity pull docker://biocontainers/fastqc:v0.11.9
    
    # 运行容器
    singularity exec fastqc_v0.11.9.sif fastqc sample.fastq
    
    # 交互式运行
    singularity shell fastqc_v0.11.9.sif
    
    # 绑定挂载数据目录
    singularity exec -B /path/to/data:/data fastqc_v0.11.9.sif fastqc /data/sample.fastq
    
    # 构建Singularity镜像
    cat > fastqc.def << 'EOF'
    Bootstrap: docker
    From: biocontainers/fastqc:v0.11.9
    
    %post
        echo "Building FastQC container"
    
    %runscript
        fastqc $@
    EOF
    
    singularity build fastqc.sif fastqc.def

.. code-block:: bash

    # Slurm作业脚本示例
    
    # 基本Slurm作业脚本
    cat > slurm_job.sh << 'EOF'
    #!/bin/bash
    #SBATCH --job-name=rna_seq_analysis
    #SBATCH --output=output_%j.txt
    #SBATCH --error=error_%j.txt
    #SBATCH --partition=normal
    #SBATCH --nodes=1
    #SBATCH --ntasks=8
    #SBATCH --mem=32G
    #SBATCH --time=24:00:00
    
    # 加载模块
    module load STAR/2.7.10a
    module load samtools/1.15
    
    # 设置变量
    INPUT_DIR="/data/fastq"
    OUTPUT_DIR="/results/alignment"
    REFERENCE="/reference/genome/STAR_index"
    
    # 运行STAR比对
    STAR --runThreadN 8 \
         --genomeDir $REFERENCE \
         --readFilesIn $INPUT_DIR/sample_R1.fastq.gz $INPUT_DIR/sample_R2.fastq.gz \
         --readFilesCommand zcat \
         --outFileNamePrefix $OUTPUT_DIR/sample_ \
         --outSAMtype BAM SortedByCoordinate
    
    echo "Alignment completed"
    EOF
    
    # 提交作业
    sbatch slurm_job.sh
    
    # 查看作业状态
    squeue -u $USER
    
    # 查看作业详情
    sacct -j <job_id> --format=JobID,JobName,Elapsed,State,MaxRSS
    
    # 取消作业
    scancel <job_id>
    
    # 查看节点信息
    sinfo
    
    # 交互式作业
    srun --pty bash

.. code-block:: bash

    # PBS作业脚本示例
    
    cat > pbs_job.sh << 'EOF'
    #!/bin/bash
    #PBS -N rna_seq_analysis
    #PBS -o output.txt
    #PBS -e error.txt
    #PBS -q normal
    #PBS -l nodes=1:ppn=8
    #PBS -l mem=32gb
    #PBS -l walltime=24:00:00
    
    # 加载模块
    module load STAR/2.7.10a
    module load samtools/1.15
    
    # 切换到工作目录
    cd $PBS_O_WORKDIR
    
    # 运行分析
    STAR --runThreadN 8 \
         --genomeDir /reference/genome/STAR_index \
         --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
         --readFilesCommand zcat \
         --outFileNamePrefix sample_ \
         --outSAMtype BAM SortedByCoordinate
    EOF
    
    # 提交作业
    qsub pbs_job.sh
    
    # 查看作业状态
    qstat -u $USER
    
    # 取消作业
    qdel <job_id>
    
    # 作业数组示例
    cat > pbs_array.sh << 'EOF'
    #!/bin/bash
    #PBS -N process_samples
    #PBS -t 1-10
    #PBS -l nodes=1:ppn=4,mem=16gb,walltime=12:00:00
    
    # 获取样本列表
    SAMPLE=$(sed -n "${PBS_ARRAYID}p" samples.txt)
    
    # 处理样本
    fastqc ${SAMPLE}.fastq
    EOF
    
    qsub pbs_array.sh

.. code-block:: python

    # 数据格式处理示例
    
    # FASTA文件处理
    def read_fasta(file_path):
        sequences = {}
        current_id = None
        current_seq = []
        
        with open(file_path, 'r') as f:
            for line in f:
                line = line.strip()
                if line.startswith('>'):
                    if current_id:
                        sequences[current_id] = ''.join(current_seq)
                    current_id = line[1:].split()[0]
                    current_seq = []
                else:
                    current_seq.append(line)
            
            if current_id:
                sequences[current_id] = ''.join(current_seq)
        
        return sequences
    
    # FASTQ文件处理
    def read_fastq(file_path):
        sequences = []
        with open(file_path, 'r') as f:
            while True:
                header = f.readline().strip()
                if not header:
                    break
                seq = f.readline().strip()
                plus = f.readline().strip()
                qual = f.readline().strip()
                
                sequences.append({
                    'id': header[1:],
                    'sequence': seq,
                    'quality': qual
                })
        
        return sequences
    
    # 解析质量值
    def parse_quality(qual_string, offset=33):
        return [ord(c) - offset for c in qual_string]
    
    # SAM/BAM文件处理
    import pysam
    
    # 读取BAM文件
    bamfile = pysam.AlignmentFile("aligned.bam", "rb")
    
    for read in bamfile:
        print(f"Query: {read.query_name}")
        print(f"Reference: {read.reference_name}")
        print(f"Position: {read.reference_start}")
        print(f"CIGAR: {read.cigarstring}")
        print(f"Flag: {read.flag}")
    
    bamfile.close()
    
    # VCF文件处理
    import vcf
    
    vcf_reader = vcf.Reader(open("variants.vcf", 'r'))
    
    for record in vcf_reader:
        print(f"CHROM: {record.CHROM}")
        print(f"POS: {record.POS}")
        print(f"REF: {record.REF}")
        print(f"ALT: {record.ALT}")
        
        for sample in record.samples:
            print(f"Sample: {sample.sample}")
            print(f"GT: {sample['GT']}")
    
    # GFF文件处理
    def read_gff(file_path):
        features = []
        with open(file_path, 'r') as f:
            for line in f:
                if line.startswith('#'):
                    continue
                
                parts = line.strip().split('\t')
                if len(parts) < 9:
                    continue
                
                feature = {
                    'seqid': parts[0],
                    'source': parts[1],
                    'type': parts[2],
                    'start': int(parts[3]),
                    'end': int(parts[4]),
                    'score': parts[5],
                    'strand': parts[6],
                    'phase': parts[7],
                    'attributes': parts[8]
                }
                
                features.append(feature)
        
        return features

数据库基础
----------

数据库是存储和管理数据的系统，它提供了数据的持久化存储、高效查询和并发访问等功能。在生物信息学中，数据库用于存储基因组数据、表达数据、变异数据等大规模生物数据。理解数据库的基本原理和使用方法是生物信息学数据处理的重要技能。

关系型数据库
~~~~~~~~~~~~

关系型数据库是最常用的数据库类型，它使用表格存储数据，表格之间通过关系连接。关系型数据库使用SQL语言进行数据操作，支持事务、索引、约束等特性，保证数据的一致性和完整性。

SQL是结构化查询语言，用于操作关系型数据库。SQL包括数据定义语言（DDL）、数据操作语言（DML）、数据查询语言（DQL）和数据控制语言（DCL）。常用的SQL语句包括CREATE创建表、INSERT插入数据、SELECT查询数据、UPDATE更新数据、DELETE删除数据。

表是关系型数据库的基本存储单元，由行和列组成。每列有特定的数据类型，如整数、浮点数、字符串、日期等。每行代表一条记录。表可以有主键，唯一标识每行记录。表之间可以通过外键建立关系。

索引是提高查询性能的重要手段。索引类似于书的目录，可以快速定位到需要的数据。索引建立在表的列上，可以是单列索引或组合索引。索引加速查询，但会降低插入、更新、删除的速度，因为索引也需要维护。

SQLite与MySQL
~~~~~~~~~~~~~

SQLite是一个轻量级的嵌入式数据库，它不需要服务器，整个数据库存储在单个文件中。SQLite适合小型应用和单机程序，在生物信息学中常用于存储分析结果和中间数据。Python标准库提供了sqlite3模块，可以直接操作SQLite数据库。

MySQL是最流行的开源关系型数据库，采用客户端-服务器架构。MySQL适合大型应用和多用户环境，支持大规模数据存储和高并发访问。MySQL提供了丰富的数据类型、存储引擎和优化选项。在生物信息学中，MySQL常用于构建数据库应用和网站后端。

连接数据库需要提供连接参数，如主机地址、端口、用户名、密码、数据库名等。Python可以使用mysql-connector-python、PyMySQL等库连接MySQL数据库。连接后可以创建游标对象，执行SQL语句，获取结果。

数据库设计应该遵循规范化原则，避免数据冗余和更新异常。常用的规范化形式包括第一范式、第二范式、第三范式等。但过度规范化会影响查询性能，需要在规范化和性能之间找到平衡。

NoSQL数据库
~~~~~~~~~~~

NoSQL数据库是非关系型数据库的统称，它们不使用表格存储数据，而是使用键值对、文档、列族、图等数据模型。NoSQL数据库适合处理大规模非结构化数据，具有良好的可扩展性和灵活性。

键值数据库是最简单的NoSQL数据库，数据以键值对形式存储。Redis是最流行的键值数据库，支持丰富的数据类型和操作。键值数据库适合缓存、会话管理等场景。

文档数据库存储文档形式的数据，通常是JSON或BSON格式。MongoDB是最流行的文档数据库，支持灵活的模式和强大的查询功能。文档数据库适合存储半结构化数据，如基因注释、样本信息等。

列族数据库按列族存储数据，适合大规模数据分析。HBase是基于Hadoop的列族数据库，可以存储PB级别的数据。列族数据库适合时序数据、日志数据等。

图数据库存储图结构数据，节点和边可以存储属性。Neo4j是最流行的图数据库，支持图查询语言Cypher。图数据库适合存储关系网络，如蛋白质相互作用网络、代谢通路等。

生物信息学数据库应用
~~~~~~~~~~~~~~~~~~~~

生物信息学数据库存储各种生物数据，包括基因组序列、基因注释、表达数据、变异数据等。这些数据库为生物信息学研究提供了重要的数据资源。

基因组数据库存储基因组序列和注释信息。Ensembl、UCSC Genome Browser、NCBI等提供了基因组数据库的访问接口。可以使用SQL或API查询基因位置、序列、注释等信息。

表达数据库存储基因表达数据。GEO、ArrayExpress等数据库存储了大量的表达谱数据。可以使用数据库查询特定基因在不同条件下的表达情况。

变异数据库存储遗传变异信息。dbSNP存储已知的SNP变异，ClinVar存储临床相关变异，gnomAD存储人群频率数据。这些数据库对变异解读和疾病研究非常重要。

构建本地数据库可以提高数据访问效率。可以下载公共数据库的数据，导入本地数据库系统。使用索引和优化查询可以快速检索需要的信息。本地数据库还可以整合多个数据源，构建定制化的数据资源。

统计学基础
----------

统计学是收集、分析、解释和展示数据的科学。在生物信息学中，统计学方法用于分析实验数据、识别显著差异、评估结果的可靠性。掌握统计学基础是进行生物信息学分析的必要条件。

描述性统计
~~~~~~~~~~

描述性统计是对数据的基本特征进行总结和描述。常用的描述性统计量包括集中趋势指标和离散程度指标。集中趋势指标包括均值、中位数、众数，它们描述数据的中心位置。离散程度指标包括方差、标准差、范围、四分位数范围，它们描述数据的分散程度。

均值是所有数值的平均，是最常用的集中趋势指标。均值对极端值敏感，当数据存在极端值时，均值可能不能很好地代表数据的中心位置。中位数是将数据排序后位于中间位置的值，对极端值不敏感。众数是出现频率最高的值。

方差是每个数据点与均值之差的平方的平均值，标准差是方差的平方根。标准差与原始数据有相同的单位，更直观地表示数据的离散程度。范围是最大值与最小值之差，四分位数范围是第三四分位数与第一四分位数之差。

数据可视化是描述性统计的重要组成部分。直方图显示数据的分布，箱线图显示数据的五数概括，散点图显示两个变量之间的关系。可视化可以帮助理解数据的特征和发现异常值。

概率分布
~~~~~~~~

概率分布描述随机变量取各个值的概率。离散型随机变量使用概率质量函数，连续型随机变量使用概率密度函数。理解常见的概率分布是进行统计推断的基础。

正态分布是最重要的连续型概率分布，也称为高斯分布。正态分布由均值和标准差两个参数决定，其概率密度函数呈钟形曲线。许多自然现象近似服从正态分布，如身高、体重、测量误差等。中心极限定理表明，大量独立随机变量的和近似服从正态分布。

二项分布描述n次独立试验中成功次数的分布，每次试验成功的概率为p。二项分布是离散型分布，常用于计数数据，如测序reads中某个碱基出现的次数。

泊松分布描述单位时间或空间内事件发生次数的分布。泊松分布常用于稀有事件的计数，如基因突变的数量。当二项分布的n很大、p很小时，二项分布近似泊松分布。

t分布用于小样本的统计推断，当样本量较小时，样本均值的分布服从t分布而不是正态分布。t分布比正态分布有更厚的尾部，随着自由度增加，t分布趋近于正态分布。

假设检验
~~~~~~~~

假设检验是统计推断的核心方法，用于根据样本数据判断关于总体的假设是否成立。假设检验包括提出假设、选择检验方法、计算检验统计量、确定P值、做出结论等步骤。

零假设是需要检验的假设，通常表示没有差异或没有效果。备择假设是与零假设对立的假设，表示存在差异或有效果。检验的结果要么拒绝零假设，要么不能拒绝零假设。不能拒绝零假设并不意味着零假设为真，只是证据不足以拒绝它。

P值是在零假设为真的条件下，观察到当前结果或更极端结果的概率。P值越小，拒绝零假设的证据越强。通常使用显著性水平作为阈值，如果P值小于显著性水平，则拒绝零假设。常用的显著性水平是0.05和0.01。

第一类错误是零假设为真时拒绝零假设，也称为假阳性。第二类错误是零假设为假时不拒绝零假设，也称为假阴性。显著性水平控制第一类错误的概率。检验的功效是正确拒绝错误零假设的概率，等于1减去第二类错误的概率。

常用的假设检验方法包括t检验、卡方检验、方差分析等。t检验用于比较两组均值。卡方检验用于分类数据的独立性检验和拟合优度检验。方差分析用于比较多组均值。

多重检验校正
~~~~~~~~~~~~

在生物信息学分析中，经常需要同时进行大量假设检验，如差异表达分析中检验数千个基因。多重检验会增加假阳性的概率，需要进行校正。

族错误率是至少出现一次第一类错误的概率。随着检验次数增加，族错误率会迅速增加。控制族错误率的方法包括Bonferroni校正、Holm方法等。Bonferroni校正将显著性水平除以检验次数，是最保守的校正方法。

错误发现率是被拒绝的假设中错误拒绝的比例。控制错误发现率的方法比控制族错误率的方法更宽松，可以发现更多的真实差异。Benjamini-Hochberg方法是常用的控制错误发现率的方法。

选择合适的校正方法取决于研究目的。如果需要严格控制假阳性，应该使用控制族错误率的方法。如果希望发现更多潜在差异，可以接受一定的假阳性，可以使用控制错误发现率的方法。

在差异表达分析中，通常使用错误发现率校正。常用的阈值是错误发现率小于0.05或0.1。这意味着在报告的差异表达基因中，假阳性的比例控制在5%或10%以内。

机器学习基础
------------

机器学习是让计算机从数据中学习规律和模式的方法。在生物信息学中，机器学习用于分类、聚类、降维、预测等任务。掌握机器学习的基本原理和方法是进行高级生物信息学分析的重要技能。

机器学习概述
~~~~~~~~~~~~

机器学习算法根据训练数据是否包含标签，分为监督学习、无监督学习和半监督学习。监督学习使用带标签的数据训练模型，用于分类和回归任务。无监督学习使用无标签的数据，用于聚类和降维任务。半监督学习使用部分带标签的数据。

监督学习的目标是学习输入到输出的映射关系。分类任务预测离散的类别标签，如疾病分类、细胞类型识别。回归任务预测连续的数值，如基因表达量预测、药物反应预测。常见的监督学习算法包括逻辑回归、支持向量机、决策树、随机森林等。

无监督学习的目标是发现数据中的结构和模式。聚类任务将相似的样本分组，如基因共表达聚类、患者分层。降维任务将高维数据映射到低维空间，如主成分分析、t-SNE可视化。常见的无监督学习算法包括K均值聚类、层次聚类、主成分分析等。

模型评估是机器学习的重要环节。对于分类任务，常用的评估指标包括准确率、精确率、召回率、F1分数、ROC曲线和AUC值。对于回归任务，常用的评估指标包括均方误差、均方根误差、平均绝对误差、决定系数等。交叉验证是评估模型泛化能力的常用方法。

特征选择与降维
~~~~~~~~~~~~~~

特征选择是从原始特征中选择最相关的特征子集。特征选择可以降低维度、提高模型性能、加速训练、增强可解释性。在生物信息学中，特征选择常用于从大量基因中筛选关键基因。

过滤方法根据特征的统计特性进行选择，独立于模型。常用的过滤方法包括方差阈值、相关系数、卡方检验、互信息等。过滤方法计算速度快，但可能忽略特征之间的交互作用。

包装方法根据模型的性能选择特征子集。常用的包装方法包括递归特征消除、前向选择、后向消除等。包装方法考虑特征之间的交互作用，但计算成本高，容易过拟合。

嵌入方法在模型训练过程中进行特征选择。LASSO回归通过L1正则化实现特征选择，决策树和随机森林可以输出特征重要性。嵌入方法结合了过滤方法和包装方法的优点。

降维是将高维数据映射到低维空间的技术。主成分分析是最常用的线性降维方法，通过正交变换找到方差最大的方向。t-SNE和UMAP是非线性降维方法，适合可视化高维数据。降维可以减少计算复杂度，发现数据的潜在结构。

常用算法介绍
~~~~~~~~~~~~

逻辑回归是经典的分类算法，使用逻辑函数将线性组合映射到0到1之间的概率。逻辑回归简单高效，可以输出概率，适合二分类问题。在生物信息学中，逻辑回归常用于疾病风险预测、生物标志物筛选等。

支持向量机是强大的分类算法，通过寻找最大间隔超平面进行分类。支持向量机可以使用核函数处理非线性问题，常用的核函数包括线性核、多项式核、径向基核等。支持向量机在小样本高维数据上表现良好，常用于基因表达数据分类。

决策树是基于树结构进行决策的算法，每个内部节点表示一个特征测试，每个叶节点表示一个类别。决策树易于理解和解释，可以处理分类和回归问题。但决策树容易过拟合，可以通过剪枝和集成学习改善。

随机森林是集成学习方法，通过构建多棵决策树并综合它们的结果进行预测。随机森林对过拟合有较强的抵抗力，可以输出特征重要性，适合高维数据。在生物信息学中，随机森林广泛用于分类和特征选择任务。

K均值聚类是最常用的聚类算法，通过迭代更新聚类中心，将样本分配到最近的聚类。K均值聚类简单高效，但需要预先指定聚类数，对初始值敏感。层次聚类通过构建层次树进行聚类，不需要预先指定聚类数。

主成分分析是最常用的降维方法，通过正交变换找到方差最大的方向。主成分分析可以用于数据可视化、降噪、特征提取。在生物信息学中，主成分分析常用于表达数据的可视化和批次效应检测。

生物信息学应用
~~~~~~~~~~~~~~

机器学习在生物信息学中有广泛的应用。在基因组学中，机器学习用于基因预测、功能注释、调控元件识别等。在转录组学中，机器学习用于差异表达分析、细胞类型识别、轨迹推断等。在蛋白质组学中，机器学习用于蛋白质结构预测、功能预测、相互作用预测等。

疾病分类是机器学习在生物信息学中的典型应用。使用基因表达数据或其他组学数据，可以构建分类模型区分疾病亚型或预测疾病风险。特征选择可以筛选与疾病相关的关键基因，为疾病机制研究提供线索。

药物发现是另一个重要应用领域。机器学习可以预测药物与靶点的相互作用、药物的活性和毒性。虚拟筛选使用机器学习模型从化合物库中筛选候选药物，大大加速药物研发过程。

单细胞数据分析是机器学习的新兴应用领域。降维和聚类用于识别细胞类型，轨迹推断用于重建细胞发育过程，批次校正用于整合多个数据集。机器学习方法在单细胞数据分析中发挥着越来越重要的作用。

.. code-block:: python

    # SQLite数据库操作示例
    import sqlite3
    
    # 连接数据库（如果不存在则创建）
    conn = sqlite3.connect('bioinfo.db')
    cursor = conn.cursor()
    
    # 创建基因表
    cursor.execute('''
    CREATE TABLE IF NOT EXISTS genes (
        gene_id TEXT PRIMARY KEY,
        gene_name TEXT NOT NULL,
        chromosome TEXT,
        start_pos INTEGER,
        end_pos INTEGER,
        strand TEXT
    )
    ''')
    
    # 插入数据
    cursor.execute('''
    INSERT INTO genes (gene_id, gene_name, chromosome, start_pos, end_pos, strand)
    VALUES (?, ?, ?, ?, ?, ?)
    ''', ('ENSG00000141510', 'TP53', 'chr17', 7668421, 7687490, '-'))
    
    # 批量插入
    genes_data = [
        ('ENSG00000012048', 'BRCA1', 'chr17', 43044295, 43125483, '-'),
        ('ENSG00000146648', 'EGFR', 'chr7', 55019021, 55211628, '+'),
        ('ENSG00000136997', 'MYC', 'chr8', 127735434, 127742951, '+')
    ]
    
    cursor.executemany('''
    INSERT INTO genes (gene_id, gene_name, chromosome, start_pos, end_pos, strand)
    VALUES (?, ?, ?, ?, ?, ?)
    ''', genes_data)
    
    # 查询数据
    cursor.execute("SELECT * FROM genes WHERE chromosome = 'chr17'")
    for row in cursor.fetchall():
        print(row)
    
    # 更新数据
    cursor.execute('''
    UPDATE genes SET strand = '-' WHERE gene_name = 'TP53'
    ''')
    
    # 删除数据
    cursor.execute("DELETE FROM genes WHERE gene_name = 'MYC'")
    
    # 提交事务
    conn.commit()
    
    # 关闭连接
    conn.close()

.. code-block:: python

    # MySQL数据库操作示例
    import pymysql
    
    # 连接MySQL数据库
    connection = pymysql.connect(
        host='localhost',
        user='bioinfo',
        password='password',
        database='bioinfo_db',
        charset='utf8mb4',
        cursorclass=pymysql.cursors.DictCursor
    )
    
    try:
        with connection.cursor() as cursor:
            # 创建表达数据表
            cursor.execute('''
            CREATE TABLE IF NOT EXISTS expression (
                id INT AUTO_INCREMENT PRIMARY KEY,
                sample_id VARCHAR(50),
                gene_id VARCHAR(50),
                expression_value FLOAT,
                FOREIGN KEY (gene_id) REFERENCES genes(gene_id)
            )
            ''')
            
            # 插入表达数据
            cursor.execute('''
            INSERT INTO expression (sample_id, gene_id, expression_value)
            VALUES (%s, %s, %s)
            ''', ('sample1', 'ENSG00000141510', 15.5))
            
            # 查询特定基因的表达
            cursor.execute('''
            SELECT sample_id, expression_value 
            FROM expression 
            WHERE gene_id = %s
            ''', ('ENSG00000141510',))
            
            results = cursor.fetchall()
            for row in results:
                print(f"Sample: {row['sample_id']}, Expression: {row['expression_value']}")
            
            # 连接查询
            cursor.execute('''
            SELECT g.gene_name, e.sample_id, e.expression_value
            FROM genes g
            JOIN expression e ON g.gene_id = e.gene_id
            WHERE e.expression_value > 10
            ''')
            
            high_expr = cursor.fetchall()
            for row in high_expr:
                print(f"{row['gene_name']}: {row['expression_value']}")
        
        connection.commit()
    
    finally:
        connection.close()

.. code-block:: python

    # 描述性统计示例
    import numpy as np
    import pandas as pd
    from scipy import stats
    
    # 生成示例数据
    expression_data = np.random.normal(10, 2, 1000)
    
    # 计算描述性统计量
    mean_val = np.mean(expression_data)
    median_val = np.median(expression_data)
    std_val = np.std(expression_data)
    var_val = np.var(expression_data)
    
    print(f"均值: {mean_val:.2f}")
    print(f"中位数: {median_val:.2f}")
    print(f"标准差: {std_val:.2f}")
    print(f"方差: {var_val:.2f}")
    
    # 使用pandas计算
    df = pd.DataFrame({'expression': expression_data})
    print(df.describe())
    
    # 分组统计
    df['group'] = np.random.choice(['control', 'treatment'], 1000)
    grouped = df.groupby('group')['expression'].agg(['mean', 'std', 'count'])
    print(grouped)
    
    # 可视化
    import matplotlib.pyplot as plt
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    # 直方图
    axes[0].hist(expression_data, bins=30, edgecolor='black')
    axes[0].set_xlabel('Expression Level')
    axes[0].set_ylabel('Frequency')
    axes[0].set_title('Histogram')
    
    # 箱线图
    df.boxplot(column='expression', by='group', ax=axes[1])
    axes[1].set_title('Boxplot by Group')
    
    # Q-Q图
    stats.probplot(expression_data, dist="norm", plot=axes[2])
    axes[2].set_title('Q-Q Plot')
    
    plt.tight_layout()
    plt.savefig('descriptive_stats.png', dpi=300)
    plt.close()

.. code-block:: python

    # 假设检验示例
    from scipy import stats
    import numpy as np
    
    # 生成示例数据
    control = np.random.normal(10, 2, 50)
    treatment = np.random.normal(12, 2, 50)
    
    # t检验（独立样本）
    t_stat, p_value = stats.ttest_ind(control, treatment)
    print(f"独立样本t检验:")
    print(f"t统计量: {t_stat:.4f}")
    print(f"P值: {p_value:.4f}")
    
    # 配对t检验
    before = np.random.normal(10, 2, 30)
    after = before + np.random.normal(1, 0.5, 30)
    
    t_stat, p_value = stats.ttest_rel(before, after)
    print(f"\n配对t检验:")
    print(f"t统计量: {t_stat:.4f}")
    print(f"P值: {p_value:.4f}")
    
    # 方差分析（ANOVA）
    group1 = np.random.normal(10, 2, 30)
    group2 = np.random.normal(12, 2, 30)
    group3 = np.random.normal(11, 2, 30)
    
    f_stat, p_value = stats.f_oneway(group1, group2, group3)
    print(f"\n单因素方差分析:")
    print(f"F统计量: {f_stat:.4f}")
    print(f"P值: {p_value:.4f}")
    
    # 卡方检验
    observed = np.array([[50, 30], [20, 40]])
    chi2, p_value, dof, expected = stats.chi2_contingency(observed)
    print(f"\n卡方检验:")
    print(f"卡方值: {chi2:.4f}")
    print(f"P值: {p_value:.4f}")
    
    # 非参数检验
    u_stat, p_value = stats.mannwhitneyu(control, treatment)
    print(f"\nMann-Whitney U检验:")
    print(f"U统计量: {u_stat:.4f}")
    print(f"P值: {p_value:.4f}")

.. code-block:: python

    # 多重检验校正示例
    import numpy as np
    from scipy import stats
    from statsmodels.stats.multitest import multipletests
    
    # 模拟差异表达分析
    n_genes = 1000
    n_samples = 10
    
    # 生成表达数据
    control = np.random.normal(10, 2, (n_genes, n_samples))
    treatment = np.random.normal(10, 2, (n_genes, n_samples))
    
    # 为部分基因添加差异
    treatment[:50, :] += 2
    
    # 对每个基因进行t检验
    p_values = []
    for i in range(n_genes):
        t_stat, p_val = stats.ttest_ind(control[i, :], treatment[i, :])
        p_values.append(p_val)
    
    p_values = np.array(p_values)
    
    # 多重检验校正
    # Bonferroni校正
    bonf_corrected = multipletests(p_values, method='bonferroni')
    
    # Benjamini-Hochberg校正（FDR）
    fdr_corrected = multipletests(p_values, method='fdr_bh')
    
    print(f"原始显著基因数 (p < 0.05): {sum(p_values < 0.05)}")
    print(f"Bonferroni校正后显著基因数: {sum(bonf_corrected[0])}")
    print(f"FDR校正后显著基因数: {sum(fdr_corrected[0])}")
    
    # 可视化P值分布
    import matplotlib.pyplot as plt
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    axes[0].hist(p_values, bins=50, edgecolor='black')
    axes[0].set_xlabel('P-value')
    axes[0].set_title('Original P-values')
    axes[0].axvline(0.05, color='red', linestyle='--')
    
    axes[1].hist(bonf_corrected[1], bins=50, edgecolor='black')
    axes[1].set_xlabel('Adjusted P-value')
    axes[1].set_title('Bonferroni Correction')
    axes[1].axvline(0.05, color='red', linestyle='--')
    
    axes[2].hist(fdr_corrected[1], bins=50, edgecolor='black')
    axes[2].set_xlabel('Adjusted P-value')
    axes[2].set_title('FDR Correction')
    axes[2].axvline(0.05, color='red', linestyle='--')
    
    plt.tight_layout()
    plt.savefig('multiple_testing.png', dpi=300)
    plt.close()

.. code-block:: python

    # 机器学习示例：分类任务
    import numpy as np
    import pandas as pd
    from sklearn.model_selection import train_test_split, cross_val_score
    from sklearn.preprocessing import StandardScaler
    from sklearn.linear_model import LogisticRegression
    from sklearn.svm import SVC
    from sklearn.ensemble import RandomForestClassifier
    from sklearn.metrics import classification_report, confusion_matrix, roc_auc_score
    from sklearn.feature_selection import SelectKBest, f_classif
    
    # 生成模拟基因表达数据
    np.random.seed(42)
    n_samples = 200
    n_genes = 1000
    
    # 生成表达矩阵
    X = np.random.randn(n_samples, n_genes)
    
    # 为前50个基因添加差异（这些是与疾病相关的基因）
    X[:100, :50] += 1.5
    
    # 创建标签（0=健康，1=疾病）
    y = np.array([0] * 100 + [1] * 100)
    
    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=42, stratify=y
    )
    
    # 特征选择
    selector = SelectKBest(f_classif, k=100)
    X_train_selected = selector.fit_transform(X_train, y_train)
    X_test_selected = selector.transform(X_test)
    
    # 标准化
    scaler = StandardScaler()
    X_train_scaled = scaler.fit_transform(X_train_selected)
    X_test_scaled = scaler.transform(X_test_selected)
    
    # 训练多个模型
    models = {
        'Logistic Regression': LogisticRegression(random_state=42),
        'SVM': SVC(kernel='rbf', random_state=42),
        'Random Forest': RandomForestClassifier(n_estimators=100, random_state=42)
    }
    
    for name, model in models.items():
        # 交叉验证
        cv_scores = cross_val_score(model, X_train_scaled, y_train, cv=5)
        print(f"\n{name}:")
        print(f"交叉验证准确率: {cv_scores.mean():.4f} (+/- {cv_scores.std():.4f})")
        
        # 训练和预测
        model.fit(X_train_scaled, y_train)
        y_pred = model.predict(X_test_scaled)
        
        # 评估
        print(f"测试集准确率: {model.score(X_test_scaled, y_test):.4f}")
        print("\n分类报告:")
        print(classification_report(y_test, y_pred))

.. code-block:: python

    # 机器学习示例：聚类和降维
    import numpy as np
    import pandas as pd
    from sklearn.cluster import KMeans
    from sklearn.decomposition import PCA
    from sklearn.manifold import TSNE
    import matplotlib.pyplot as plt
    
    # 生成模拟单细胞表达数据
    np.random.seed(42)
    n_cells = 500
    n_genes = 2000
    
    # 生成3个细胞类型
    cell_type1 = np.random.normal(5, 1, (150, n_genes))
    cell_type2 = np.random.normal(8, 1.5, (200, n_genes))
    cell_type3 = np.random.normal(3, 0.8, (150, n_genes))
    
    X = np.vstack([cell_type1, cell_type2, cell_type3])
    true_labels = np.array([0] * 150 + [1] * 200 + [2] * 150)
    
    # PCA降维
    pca = PCA(n_components=50)
    X_pca = pca.fit_transform(X)
    
    print(f"前50个主成分解释的方差比例: {sum(pca.explained_variance_ratio_):.4f}")
    
    # K-means聚类
    kmeans = KMeans(n_clusters=3, random_state=42)
    cluster_labels = kmeans.fit_predict(X_pca)
    
    # t-SNE可视化
    tsne = TSNE(n_components=2, random_state=42)
    X_tsne = tsne.fit_transform(X_pca)
    
    # 可视化结果
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    # PCA结果
    scatter1 = axes[0].scatter(X_pca[:, 0], X_pca[:, 1], c=true_labels, cmap='viridis', alpha=0.6)
    axes[0].set_xlabel('PC1')
    axes[0].set_ylabel('PC2')
    axes[0].set_title('PCA - True Labels')
    plt.colorbar(scatter1, ax=axes[0])
    
    # t-SNE结果（真实标签）
    scatter2 = axes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], c=true_labels, cmap='viridis', alpha=0.6)
    axes[1].set_xlabel('t-SNE 1')
    axes[1].set_ylabel('t-SNE 2')
    axes[1].set_title('t-SNE - True Labels')
    plt.colorbar(scatter2, ax=axes[1])
    
    # t-SNE结果（聚类标签）
    scatter3 = axes[2].scatter(X_tsne[:, 0], X_tsne[:, 1], c=cluster_labels, cmap='viridis', alpha=0.6)
    axes[2].set_xlabel('t-SNE 1')
    axes[2].set_ylabel('t-SNE 2')
    axes[2].set_title('t-SNE - K-means Clusters')
    plt.colorbar(scatter3, ax=axes[2])
    
    plt.tight_layout()
    plt.savefig('clustering_visualization.png', dpi=300)
    plt.close()
    
    # 评估聚类质量
    from sklearn.metrics import adjusted_rand_score, silhouette_score
    
    ari = adjusted_rand_score(true_labels, cluster_labels)
    silhouette = silhouette_score(X_pca, cluster_labels)
    
    print(f"调整兰德指数: {ari:.4f}")
    print(f"轮廓系数: {silhouette:.4f}")

网络基础与API
-------------

网络通信是现代生物信息学分析的重要组成部分。理解网络基础和API使用方法，可以方便地获取公共数据库的数据、调用在线分析工具、构建网络应用等。掌握网络编程基础是生物信息学工作者的必备技能。

网络基础概念
~~~~~~~~~~~~

计算机网络是多台计算机通过通信线路连接起来，实现资源共享和信息传递的系统。网络通信遵循一定的协议，最常用的是TCP/IP协议族。TCP/IP协议分为四层：应用层、传输层、网络层和链路层。

IP地址是网络中设备的唯一标识。IPv4地址是32位地址，通常表示为四个十进制数，如192.168.1.1。IPv6地址是128位地址，解决了IPv4地址不足的问题。域名是IP地址的易记形式，通过DNS解析为IP地址。

端口是网络通信的逻辑端点，每个端口对应一个服务。常用的端口包括HTTP的80端口、HTTPS的443端口、SSH的22端口等。端口号范围是0到65535，其中0到1023是知名端口，1024到49151是注册端口，49152到65535是动态端口。

HTTP协议是应用层最重要的协议，用于Web通信。HTTP请求包括请求方法、请求头、请求体。常用的请求方法包括GET获取资源、POST提交数据、PUT更新资源、DELETE删除资源。HTTP响应包括状态码、响应头、响应体。常用的状态码包括200成功、301重定向、404未找到、500服务器错误。

RESTful API
~~~~~~~~~~~

REST是一种软件架构风格，RESTful API是符合REST原则的Web API。RESTful API使用HTTP方法表示操作，使用URL表示资源，使用JSON或XML传输数据。RESTful API设计简洁、易于理解、扩展性好。

URL是统一资源定位符，用于标识网络资源。RESTful API的URL应该表示资源而不是动作。例如，获取基因列表的URL应该是/genes而不是/getGenes。资源可以是单个资源或资源集合，如/genes表示基因集合，/genes/TP53表示特定基因。

HTTP方法对应CRUD操作。GET用于读取资源，POST用于创建资源，PUT用于更新资源，DELETE用于删除资源。GET和DELETE请求通常没有请求体，POST和PUT请求通常包含JSON格式的请求体。

API认证用于验证客户端身份。常用的认证方式包括基本认证、API密钥、OAuth等。基本认证使用用户名和密码，API密钥使用预分配的密钥，OAuth使用令牌进行授权。认证信息通常放在请求头中。

Python网络请求
~~~~~~~~~~~~~~

Python提供了多个库用于发送HTTP请求。requests库是最常用的HTTP客户端库，提供了简洁的API。urllib是Python标准库，功能完整但使用稍复杂。httpx是新一代HTTP客户端，支持HTTP/2和异步请求。

发送GET请求使用requests.get方法，可以传递参数、请求头、认证信息等。响应对象包含状态码、响应头、响应体等信息。响应体可以使用text属性获取文本，使用json方法解析JSON，使用content属性获取字节。

发送POST请求使用requests.post方法，需要传递请求体。请求体可以是字典、JSON字符串、文件等。requests库会自动处理编码和内容类型。上传文件使用files参数，可以上传单个文件或多个文件。

处理响应时需要检查状态码，判断请求是否成功。可以使用raise_for_status方法，如果状态码表示错误则抛出异常。对于JSON响应，使用json方法解析。对于大文件下载，使用流式传输避免内存溢出。

生物信息学API应用
~~~~~~~~~~~~~~~~~

许多生物信息学数据库和工具提供了API接口。NCBI提供了E-utilities API，可以查询PubMed、Gene、Protein等数据库。Ensembl提供了REST API，可以查询基因组注释、变异数据等。UniProt提供了API，可以查询蛋白质序列和功能注释。

使用API可以自动化数据获取流程。例如，批量查询基因信息、下载序列数据、获取注释信息等。API通常有访问频率限制，需要合理控制请求频率。可以使用缓存减少重复请求，使用异步请求提高效率。

NCBI E-utilities包括多个工具：ESearch搜索记录、EFetch获取记录、ESummary获取摘要、ELink获取链接。使用E-utilities需要提供数据库类型、查询词、返回格式等参数。NCBI要求提供API密钥，使用密钥可以获得更高的访问频率。

Ensembl REST API提供了丰富的端点，包括基因组序列、基因注释、变异信息、比较基因组学等。API支持多种返回格式，包括JSON、XML、BED等。Ensembl还提供了Perl和Python客户端库，简化API调用。

数据可视化
----------

数据可视化是将数据转换为图形或图像的过程，帮助人们理解数据、发现模式、传达信息。在生物信息学中，数据可视化用于展示基因表达、基因组特征、进化关系、蛋白质结构等。掌握数据可视化技术是生物信息学分析的重要技能。

可视化基础
~~~~~~~~~~

数据可视化遵循一定的设计原则。好的可视化应该准确表达数据、清晰传达信息、美观易于理解。可视化类型应该根据数据类型和分析目的选择。常用的可视化类型包括柱状图、折线图、散点图、饼图、热图、树状图等。

柱状图用于比较分类数据的数值。折线图用于展示连续数据的变化趋势。散点图用于展示两个变量之间的关系。饼图用于展示各部分占总体的比例。热图用于展示矩阵数据的模式。树状图用于展示层次关系。

颜色是可视化的重要元素。颜色可以区分类别、表示数值、突出重点。选择颜色时需要考虑色盲友好、打印友好、文化含义等因素。常用的配色方案包括分类配色、顺序配色、发散配色等。

坐标轴和标签帮助读者理解图形。坐标轴应该有清晰的标题和刻度。标签应该简洁明了，避免冗余信息。图例应该放在合适的位置，不遮挡数据。标题应该概括图形的主要内容。

Matplotlib基础
~~~~~~~~~~~~~~

Matplotlib是Python最基础的绘图库，提供了类似MATLAB的绘图接口。Matplotlib可以创建各种类型的图形，支持高度自定义。大多数Python可视化库都基于Matplotlib构建。

Figure是整个图形对象，可以包含多个子图。Axes是单个子图，包含坐标轴、数据、标签等元素。Artist是图形中的所有可见元素，包括线条、文本、图例等。理解这些概念有助于灵活使用Matplotlib。

基本绘图使用pyplot模块。plot函数绘制折线图，scatter函数绘制散点图，bar函数绘制柱状图，hist函数绘制直方图。每个函数都有多个参数控制图形外观，如颜色、线型、标记、大小等。

子图布局可以使用subplot函数创建规则网格，使用gridspec创建复杂布局，使用add_axes手动定位。调整子图间距使用subplots_adjust函数。保存图形使用savefig函数，支持多种格式和分辨率。

Seaborn高级可视化
~~~~~~~~~~~~~~~~~

Seaborn是基于Matplotlib的高级可视化库，提供了更美观的默认样式和更简单的API。Seaborn特别适合统计数据可视化，内置了多种统计图形类型。

分布图用于展示数据的分布。distplot绘制直方图和核密度估计，kdeplot绘制核密度估计，rugplot绘制边际分布。这些函数可以组合使用，全面展示数据分布。

关系图用于展示变量之间的关系。scatterplot绘制散点图，lineplot绘制折线图，relplot是关系图的通用接口。可以添加颜色、大小、样式等维度表示更多变量。

分类图用于展示分类数据的分布。boxplot绘制箱线图，violinplot绘制小提琴图，stripplot绘制散点图，swarmplot绘制蜂群图。这些图形可以组合使用，如violinplot和stripplot组合。

热图和聚类图在生物信息学中广泛使用。heatmap绘制热图，clustermap绘制带聚类的热图。热图适合展示表达矩阵、相关性矩阵等数据。聚类可以发现数据的模式和结构。

基因组可视化
~~~~~~~~~~~~

基因组可视化展示基因组特征和测序数据。常用的基因组可视化工具包括IGV、UCSC Genome Browser、JBrowse等。Python中可以使用pyGenomeTracks、Dalliance等库进行基因组可视化。

染色体视图展示整个染色体的特征分布。可以使用Ideogram绘制染色体带型图，使用条形图展示基因密度、GC含量等特征。染色体视图帮助理解基因组的宏观结构。

基因视图展示特定区域的详细注释。可以绘制基因结构、转录本、外显子、内含子等特征。不同类型的特征使用不同的形状和颜色。基因视图常用于展示候选基因区域的注释。

测序数据可视化展示比对结果和覆盖度。可以绘制reads堆叠图、覆盖度曲线、变异位点等。测序数据可视化帮助评估数据质量、发现变异、验证结果。

生物信息学常用工具
------------------

生物信息学分析依赖大量的专业工具和软件。这些工具涵盖了序列比对、变异检测、表达分析、功能注释等各个方面。了解常用工具的原理和使用方法是进行生物信息学分析的基础。

序列比对工具
~~~~~~~~~~~~

序列比对是生物信息学的基础操作，用于寻找序列之间的相似性。序列比对可以分为全局比对、局部比对、成对比对、多序列比对等类型。不同的比对工具适用于不同的场景。

BLAST是最常用的序列比对工具，用于在数据库中搜索相似序列。BLAST包括多个程序：blastn用于核酸序列，blastp用于蛋白质序列，blastx将核酸翻译后比对蛋白质数据库，tblastn用蛋白质序列比对核酸翻译数据库。BLAST使用启发式算法，速度快但可能遗漏一些相似性较低的序列。

BWA是常用的短序列比对工具，用于将测序reads比对到参考基因组。BWA-MEM适合较长的reads，BWA-ALN适合较短的reads。BWA使用Burrows-Wheeler变换，内存效率高，比对速度快。

Bowtie2是另一个流行的短序列比对工具，特别适合超长reads的比对。Bowtie2支持间隙比对，可以处理插入缺失变异。Bowtie2提供多种比对模式，平衡速度和灵敏度。

STAR是RNA-seq专用的比对工具，使用未压缩的后缀数组进行快速比对。STAR可以发现剪接点，适合真核生物的转录组分析。STAR比对速度快，但内存需求较大。

变异检测工具
~~~~~~~~~~~~

变异检测是从测序数据中识别遗传变异的过程。常见的变异类型包括单核苷酸多态性（SNP）、插入缺失（Indel）、结构变异（SV）、拷贝数变异（CNV）等。不同的变异类型需要不同的检测工具。

GATK是最广泛使用的变异检测工具集，提供了完整的变异检测流程。GATK Best Practices定义了标准的分析流程，包括数据预处理、变异检测、变异过滤等步骤。GATK使用HaplotypeCaller进行变异检测，可以同时检测SNP和Indel。

FreeBayes是基于单倍型的变异检测工具，使用贝叶斯方法进行变异识别。FreeBayes不需要复杂的参数设置，适合多种实验设计。FreeBayes可以处理多倍体和混合样本。

Samtools和BCFtools提供了简单的变异检测方法。mpileup命令生成pileup格式文件，call命令进行变异检测。这些工具适合快速分析和小规模数据。

结构变异检测需要专门的工具。Delly使用配对reads和分裂reads检测结构变异。Manta提供快速的SV检测流程。CNVnator使用读深度方法检测拷贝数变异。

表达分析工具
~~~~~~~~~~~~

基因表达分析是转录组研究的核心内容。表达分析包括表达定量、差异表达分析、功能富集分析等步骤。不同的测序平台和分析目的需要不同的工具。

HISAT2是新一代的转录组比对工具，使用分层索引实现快速比对。HISAT2内存需求小，比对速度快，可以处理大规模转录组数据。HISAT2可以替代TopHat2进行RNA-seq分析。

featureCounts是快速的表达定量工具，使用比对结果和基因注释计算表达量。featureCounts支持多种特征类型和重叠处理策略。featureCounts速度快，内存效率高。

Salmon和Kallisto是无比对的表达定量工具，使用伪比对技术直接从reads估计表达量。这些工具速度极快，可以在几分钟内完成定量。但需要注意参考转录本的完整性。

DESeq2和edgeR是常用的差异表达分析工具，使用负二项分布模型处理计数数据。这些工具考虑了生物学变异和技术变异，可以准确识别差异表达基因。limma-voom是另一种选择，使用线性模型进行分析。

功能注释工具
~~~~~~~~~~~~

功能注释是将生物学意义赋予基因或蛋白质的过程。功能注释包括基因本体论注释、通路注释、蛋白质结构域注释等。功能注释帮助理解基因的功能和调控机制。

BLAST可以用于功能注释，通过序列相似性推断功能。但需要注意同源基因可能具有不同的功能。保守的结构域通常与保守的功能相关。

InterProScan整合了多个蛋白质特征数据库，可以识别蛋白质结构域、家族、位点等。InterProScan提供全面的蛋白质功能注释，是蛋白质功能分析的标准工具。

GO注释使用基因本体论描述基因功能。GO分为三个本体：分子功能、生物过程、细胞组分。GO富集分析可以识别差异表达基因的功能特征。

KEGG和Reactome是常用的通路数据库。通路富集分析可以识别基因参与的生物学通路。GSEA是基因集富集分析方法，不需要预先定义差异基因。

流程管理与自动化
----------------

生物信息学分析通常涉及多个步骤和多个工具的组合。流程管理系统可以帮助组织分析步骤、跟踪数据流向、实现自动化执行。掌握流程管理技术可以提高分析效率和可重复性。

流程设计原则
~~~~~~~~~~~~

好的分析流程应该具有可重复性、可扩展性、可维护性。可重复性意味着相同的输入和参数应该产生相同的输出。可扩展性意味着流程可以处理不同规模的数据。可维护性意味着流程易于理解和修改。

模块化设计将复杂流程分解为独立的模块。每个模块完成特定的功能，有明确的输入和输出。模块之间通过数据文件连接，可以独立测试和调试。模块化提高了流程的灵活性和可维护性。

版本控制是流程管理的重要组成部分。使用Git管理流程代码，记录每次修改的原因和内容。使用版本标签标记重要的里程碑。使用分支开发新功能，避免影响主流程。

文档是流程可维护性的关键。流程应该有清晰的README文档，说明流程的目的、依赖、使用方法。每个步骤应该有注释，说明参数的选择和处理逻辑。文档帮助其他研究者理解和使用流程。

Snakemake流程管理
~~~~~~~~~~~~~~~~~

Snakemake是流行的流程管理系统，使用Python语法定义工作流。Snakemake基于规则定义分析步骤，自动解析依赖关系，支持并行执行和集群提交。Snakemake适合复杂的生物信息学分析流程。

规则是Snakemake的基本单元，定义了输入、输出和执行命令。规则之间通过文件名连接，Snakemake自动构建依赖图。当输入文件更新时，Snakemake自动重新执行相关的规则。

配置文件将参数与流程代码分离。可以使用YAML或JSON格式的配置文件存储参数。配置文件使得流程可以适应不同的数据集和分析需求，而不需要修改代码。

Snakemake支持多种执行模式。本地执行适合小规模数据，多线程执行利用多核CPU。集群提交适合大规模数据，可以自动提交作业到Slurm、PBS等调度系统。云端执行适合弹性计算需求。

Nextflow流程管理
~~~~~~~~~~~~~~~~

Nextflow是另一个流行的流程管理系统，使用Groovy语法定义工作流。Nextflow支持多种执行平台，包括本地、集群、云平台。Nextflow与容器技术集成良好，可以保证流程的可移植性。

Nextflow使用进程定义分析步骤，使用通道传递数据。进程是独立的执行单元，可以并行运行。通道是异步队列，连接进程的输入和输出。这种模型适合数据驱动的分析流程。

Nextflow支持多种容器技术。Docker容器提供了完整的运行环境，Singularity容器适合HPC环境。使用容器可以避免软件依赖问题，保证流程在不同系统上的一致性。

Nextflow提供了丰富的执行器。本地执行器适合开发和测试，Slurm、PBS、LSF等执行器适合集群环境，AWS、Google Cloud、Azure等执行器适合云平台。Nextflow可以无缝切换执行平台。

科学计算与高性能计算
--------------------

生物信息学分析经常需要处理大规模数据和复杂计算。科学计算和高性能计算技术可以加速分析过程，处理更大的数据集。了解并行计算和分布式计算的原理是进行高级生物信息学分析的必要条件。

并行计算基础
~~~~~~~~~~~~

并行计算是同时使用多个计算资源解决问题的方法。并行计算可以分为数据并行和任务并行。数据并行将数据分割到多个处理器，每个处理器执行相同的操作。任务并行将不同的任务分配给不同的处理器。

并行计算面临的主要挑战是通信开销和负载均衡。处理器之间需要交换数据，通信开销可能抵消并行带来的加速。负载不均衡会导致部分处理器空闲，降低整体效率。好的并行算法需要最小化通信，均衡负载。

多线程是最简单的并行方式，在单个进程内创建多个线程。Python的threading模块提供了多线程支持，但由于全局解释器锁（GIL），多线程不适合CPU密集型任务。multiprocessing模块使用多进程绕过GIL限制。

多进程使用多个进程并行执行，每个进程有独立的内存空间。Python的multiprocessing模块提供了进程池、共享内存、进程间通信等功能。多进程适合CPU密集型任务，但进程创建和通信开销较大。

GPU计算
~~~~~~~

GPU是图形处理器，具有大量的计算核心，适合并行计算。GPU计算在深度学习、分子动力学、基因组比对等领域有广泛应用。使用GPU可以大幅加速计算密集型任务。

CUDA是NVIDIA的GPU编程平台，提供了C/C++编程接口。Python可以通过PyCUDA、CuPy等库使用GPU。这些库提供了类似NumPy的接口，可以方便地将CPU代码迁移到GPU。

深度学习框架如TensorFlow、PyTorch提供了GPU加速。这些框架自动将计算任务分配到GPU，用户不需要编写GPU代码。GPU加速可以显著缩短模型训练时间。

GPU计算需要考虑数据传输开销。CPU和GPU之间的数据传输是瓶颈，应该尽量减少传输次数。将数据一次性传输到GPU，完成所有计算后再传输回CPU，可以最大化GPU的优势。

分布式计算
~~~~~~~~~~

分布式计算使用多台计算机协同解决问题。分布式计算可以处理单机无法容纳的数据，提供更高的计算能力。Hadoop和Spark是常用的分布式计算框架。

Hadoop使用HDFS存储数据，使用MapReduce进行计算。HDFS将数据分块存储在多个节点，提供高容错性。MapReduce将计算分为Map和Reduce两个阶段，适合批处理任务。

Spark是新一代分布式计算框架，使用内存计算提高性能。Spark提供了丰富的API，包括RDD、DataFrame、Dataset等。Spark支持SQL查询、流处理、机器学习、图计算等多种计算模式。

Dask是Python的分布式计算库，提供了类似NumPy和Pandas的接口。Dask可以处理大于内存的数据集，利用多核和集群进行并行计算。Dask适合数据科学和科学计算应用。

.. code-block:: python

    # 网络请求示例：获取生物信息学数据
    import requests
    import json
    
    # NCBI E-utilities API示例
    def search_pubmed(query, max_results=10):
        base_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"
        
        # 搜索PubMed
        search_url = f"{base_url}/esearch.fcgi"
        params = {
            'db': 'pubmed',
            'term': query,
            'retmax': max_results,
            'retmode': 'json'
        }
        
        response = requests.get(search_url, params=params)
        data = response.json()
        
        return data['esearchresult']['idlist']
    
    # 获取PubMed记录详情
    def fetch_pubmed(pmids):
        base_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"
        
        fetch_url = f"{base_url}/efetch.fcgi"
        params = {
            'db': 'pubmed',
            'id': ','.join(pmids),
            'rettype': 'abstract',
            'retmode': 'text'
        }
        
        response = requests.get(fetch_url, params=params)
        return response.text
    
    # Ensembl REST API示例
    def get_gene_info(gene_symbol, species='human'):
        server = "https://rest.ensembl.org"
        endpoint = f"/lookup/symbol/{species}/{gene_symbol}"
        
        response = requests.get(
            server + endpoint,
            headers={"Content-Type": "application/json"}
        )
        
        if response.ok:
            return response.json()
        else:
            return None
    
    # 获取序列
    def get_sequence(region, species='human'):
        server = "https://rest.ensembl.org"
        endpoint = f"/sequence/region/{species}/{region}"
        
        response = requests.get(
            server + endpoint,
            headers={"Content-Type": "application/json"}
        )
        
        if response.ok:
            return response.json()
        else:
            return None
    
    # 使用示例
    if __name__ == "__main__":
        # 搜索TP53相关文献
        pmids = search_pubmed("TP53 cancer", max_results=5)
        print(f"找到 {len(pmids)} 篇文献")
        
        # 获取基因信息
        gene_info = get_gene_info("TP53")
        if gene_info:
            print(f"基因ID: {gene_info['id']}")
            print(f"描述: {gene_info['description']}")
            print(f"位置: {gene_info['seq_region_name']}:{gene_info['start']}-{gene_info['end']}")
        
        # 获取序列
        seq_data = get_sequence("chr17:7668421-7687490")
        if seq_data:
            print(f"序列长度: {len(seq_data['seq'])} bp")

.. code-block:: python

    # 数据可视化示例
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import seaborn as sns
    
    # 设置样式
    sns.set_style("whitegrid")
    plt.rcParams['figure.figsize'] = (10, 6)
    
    # 生成模拟表达数据
    np.random.seed(42)
    genes = [f'Gene_{i}' for i in range(1, 51)]
    samples = [f'Sample_{i}' for i in range(1, 11)]
    
    expression_data = np.random.randn(50, 10) * 2 + 10
    expression_df = pd.DataFrame(expression_data, index=genes, columns=samples)
    
    # 添加分组信息
    groups = ['Control'] * 5 + ['Treatment'] * 5
    sample_info = pd.DataFrame({'group': groups}, index=samples)
    
    # 热图
    fig, ax = plt.subplots(figsize=(12, 10))
    sns.heatmap(
        expression_df.iloc[:20],
        cmap='RdBu_r',
        center=0,
        annot=False,
        cbar_kws={'label': 'Expression Level'}
    )
    ax.set_title('Gene Expression Heatmap')
    ax.set_xlabel('Samples')
    ax.set_ylabel('Genes')
    plt.tight_layout()
    plt.savefig('heatmap.png', dpi=300)
    plt.close()
    
    # 箱线图
    fig, ax = plt.subplots(figsize=(12, 6))
    melted_df = expression_df.T.melt(var_name='Gene', value_name='Expression')
    melted_df = melted_df[melted_df['Gene'].isin(genes[:10])]
    
    sns.boxplot(data=melted_df, x='Gene', y='Expression', ax=ax)
    ax.set_title('Expression Distribution by Gene')
    ax.set_xlabel('Gene')
    ax.set_ylabel('Expression Level')
    plt.xticks(rotation=45)
    plt.tight_layout()
    plt.savefig('boxplot.png', dpi=300)
    plt.close()
    
    # 火山图
    log2fc = np.random.randn(1000) * 2
    neg_log10_pval = -np.log10(np.random.uniform(0.0001, 1, 1000))
    
    volcano_df = pd.DataFrame({
        'log2FC': log2fc,
        'neg_log10_pval': neg_log10_pval
    })
    
    volcano_df['significant'] = (
        (volcano_df['log2FC'].abs() > 1) & 
        (volcano_df['neg_log10_pval'] > -np.log10(0.05))
    )
    
    fig, ax = plt.subplots(figsize=(10, 8))
    sns.scatterplot(
        data=volcano_df,
        x='log2FC',
        y='neg_log10_pval',
        hue='significant',
        palette={True: 'red', False: 'gray'},
        alpha=0.6,
        ax=ax
    )
    ax.axhline(-np.log10(0.05), color='blue', linestyle='--', label='p=0.05')
    ax.axvline(-1, color='green', linestyle='--')
    ax.axvline(1, color='green', linestyle='--')
    ax.set_title('Volcano Plot')
    ax.set_xlabel('log2 Fold Change')
    ax.set_ylabel('-log10(p-value)')
    plt.legend()
    plt.tight_layout()
    plt.savefig('volcano.png', dpi=300)
    plt.close()
    
    # 主成分分析可视化
    from sklearn.decomposition import PCA
    from sklearn.preprocessing import StandardScaler
    
    # 标准化数据
    scaler = StandardScaler()
    scaled_data = scaler.fit_transform(expression_df.T)
    
    # PCA
    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(scaled_data)
    
    pca_df = pd.DataFrame(
        pca_result,
        columns=['PC1', 'PC2'],
        index=samples
    )
    pca_df['group'] = groups
    
    # 绘制PCA图
    fig, ax = plt.subplots(figsize=(10, 8))
    sns.scatterplot(
        data=pca_df,
        x='PC1',
        y='PC2',
        hue='group',
        s=100,
        ax=ax
    )
    ax.set_title(f'PCA Plot (PC1: {pca.explained_variance_ratio_[0]:.2%}, PC2: {pca.explained_variance_ratio_[1]:.2%})')
    ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.2%})')
    ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.2%})')
    plt.legend()
    plt.tight_layout()
    plt.savefig('pca.png', dpi=300)
    plt.close()

.. code-block:: python

    # Snakemake流程示例
    # 文件: Snakefile
    
    # 配置文件
    configfile: "config.yaml"
    
    # 规则：质量控制
    rule fastqc:
        input:
            "data/{sample}.fastq.gz"
        output:
            html="results/fastqc/{sample}_fastqc.html",
            zip="results/fastqc/{sample}_fastqc.zip"
        threads: 2
        shell:
            """
            fastqc -t {threads} -o results/fastqc {input}
            """
    
    # 规则：比对
    rule align:
        input:
            reads="data/{sample}.fastq.gz",
            index=config["reference"] + ".bwt"
        output:
            "results/bam/{sample}.bam"
        threads: 8
        params:
            reference=config["reference"]
        shell:
            """
            bwa mem -t {threads} {params.reference} {input.reads} | \
            samtools view -Sb - > {output}
            """
    
    # 规则：排序
    rule sort:
        input:
            "results/bam/{sample}.bam"
        output:
            "results/bam/{sample}.sorted.bam"
        threads: 4
        shell:
            """
            samtools sort -@ {threads} -o {output} {input}
            """
    
    # 规则：变异检测
    rule call_variants:
        input:
            bam="results/bam/{sample}.sorted.bam",
            reference=config["reference"]
        output:
            "results/vcf/{sample}.vcf"
        threads: 4
        shell:
            """
            bcftools mpileup -f {input.reference} {input.bam} | \
            bcftools call -mv -o {output}
            """
    
    # 规则：生成报告
    rule report:
        input:
            expand("results/vcf/{sample}.vcf", sample=config["samples"])
        output:
            "results/report.html"
        shell:
            """
            python generate_report.py {input} > {output}
            """

.. code-block:: python

    # 并行计算示例
    import numpy as np
    from multiprocessing import Pool, cpu_count
    import time
    
    # 模拟计算密集型任务
    def process_gene(gene_id):
        # 模拟复杂计算
        time.sleep(0.1)
        return {
            'gene': gene_id,
            'score': np.random.randn(),
            'p_value': np.random.uniform(0, 1)
        }
    
    # 串行处理
    def serial_processing(gene_list):
        start_time = time.time()
        results = [process_gene(gene) for gene in gene_list]
        end_time = time.time()
        print(f"串行处理时间: {end_time - start_time:.2f}秒")
        return results
    
    # 并行处理
    def parallel_processing(gene_list, n_processes=None):
        if n_processes is None:
            n_processes = cpu_count()
        
        start_time = time.time()
        with Pool(n_processes) as pool:
            results = pool.map(process_gene, gene_list)
        end_time = time.time()
        print(f"并行处理时间 ({n_processes}进程): {end_time - start_time:.2f}秒")
        return results
    
    # 使用示例
    if __name__ == "__main__":
        gene_list = [f'Gene_{i}' for i in range(100)]
        
        print("串行处理:")
        serial_results = serial_processing(gene_list)
        
        print("\n并行处理:")
        parallel_results = parallel_processing(gene_list)
        
        print(f"\n加速比: {len(gene_list) * 0.1 / (len(gene_list) * 0.1 / cpu_count()):.2f}x")

.. code-block:: python

    # Dask并行计算示例
    import dask.array as da
    import dask.dataframe as dd
    import dask.bag as db
    from dask.distributed import Client
    import numpy as np
    
    # 启动本地集群
    client = Client(n_workers=4)
    
    # Dask Array示例：处理大型矩阵
    # 创建大型随机数组（不会立即计算）
    large_array = da.random.random((100000, 10000), chunks=(1000, 1000))
    
    # 计算统计量（延迟计算）
    mean_val = large_array.mean()
    std_val = large_array.std()
    
    # 执行计算
    print(f"数组均值: {mean_val.compute():.4f}")
    print(f"数组标准差: {std_val.compute():.4f}")
    
    # 矩阵运算
    result = da.dot(large_array[:100, :], large_array[:100, :].T)
    print(f"矩阵乘法结果形状: {result.compute().shape}")
    
    # Dask DataFrame示例：处理大型表格
    # 创建大型DataFrame
    n_rows = 1000000
    df = dd.from_pandas(
        pd.DataFrame({
            'gene': np.random.choice(['TP53', 'BRCA1', 'EGFR', 'MYC', 'KRAS'], n_rows),
            'sample': [f'Sample_{i}' for i in range(n_rows)],
            'expression': np.random.randn(n_rows) * 10 + 50,
            'group': np.random.choice(['Control', 'Treatment'], n_rows)
        }),
        npartitions=10
    )
    
    # 分组统计
    grouped = df.groupby('gene')['expression'].mean()
    print("\n基因平均表达量:")
    print(grouped.compute())
    
    # 筛选和聚合
    high_expr = df[df['expression'] > 60].groupby('group').count()
    print("\n高表达基因数量（按组）:")
    print(high_expr.compute())
    
    # Dask Bag示例：处理非结构化数据
    # 模拟FASTQ文件处理
    fastq_records = db.from_sequence([
        f"@SEQ_{i}\nATCGATCGATCG\n+\nIIIIIIIIIIII"
        for i in range(10000)
    ])
    
    # 解析FASTQ记录
    def parse_fastq(record):
        lines = record.split('\n')
        return {
            'id': lines[0][1:],
            'sequence': lines[1],
            'quality': lines[3]
        }
    
    parsed = fastq_records.map(parse_fastq)
    
    # 计算序列长度分布
    seq_lengths = parsed.map(lambda x: len(x['sequence']))
    print(f"\n平均序列长度: {seq_lengths.mean().compute():.2f}")
    
    # 关闭客户端
    client.close()

.. code-block:: python

    # 生物信息学工具调用示例
    import subprocess
    import os
    
    # 运行FastQC
    def run_fastqc(input_file, output_dir):
        cmd = [
            'fastqc',
            '-o', output_dir,
            '-t', '4',
            input_file
        ]
        
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True
        )
        
        if result.returncode == 0:
            print(f"FastQC完成: {input_file}")
        else:
            print(f"FastQC错误: {result.stderr}")
        
        return result.returncode
    
    # 运行BWA比对
    def run_bwa_mem(reference, reads, output_bam, threads=8):
        # 比对
        bwa_cmd = [
            'bwa', 'mem',
            '-t', str(threads),
            reference,
            reads
        ]
        
        # 排序
        samtools_cmd = [
            'samtools', 'sort',
            '-@', str(threads),
            '-o', output_bam
        ]
        
        # 管道连接
        bwa_process = subprocess.Popen(
            bwa_cmd,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE
        )
        
        samtools_process = subprocess.Popen(
            samtools_cmd,
            stdin=bwa_process.stdout,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE
        )
        
        bwa_process.stdout.close()
        stdout, stderr = samtools_process.communicate()
        
        if samtools_process.returncode == 0:
            print(f"比对完成: {output_bam}")
        else:
            print(f"比对错误: {stderr.decode()}")
        
        return samtools_process.returncode
    
    # 运行GATK变异检测
    def run_gatk_haplotypecaller(reference, input_bam, output_vcf):
        cmd = [
            'gatk', 'HaplotypeCaller',
            '-R', reference,
            '-I', input_bam,
            '-O', output_vcf
        ]
        
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True
        )
        
        if result.returncode == 0:
            print(f"变异检测完成: {output_vcf}")
        else:
            print(f"GATK错误: {result.stderr}")
        
        return result.returncode
    
    # 批量处理样本
    def batch_process(samples, reference, output_dir):
        for sample in samples:
            print(f"\n处理样本: {sample}")
            
            # 质量控制
            run_fastqc(
                f"data/{sample}.fastq.gz",
                f"{output_dir}/fastqc"
            )
            
            # 比对
            run_bwa_mem(
                reference,
                f"data/{sample}.fastq.gz",
                f"{output_dir}/bam/{sample}.bam"
            )
            
            # 变异检测
            run_gatk_haplotypecaller(
                reference,
                f"{output_dir}/bam/{sample}.bam",
                f"{output_dir}/vcf/{sample}.vcf"
            )
    
    # 使用示例
    if __name__ == "__main__":
        samples = ['sample1', 'sample2', 'sample3']
        reference = 'reference/hg38.fa'
        output_dir = 'results'
        
        # 创建输出目录
        os.makedirs(f"{output_dir}/fastqc", exist_ok=True)
        os.makedirs(f"{output_dir}/bam", exist_ok=True)
        os.makedirs(f"{output_dir}/vcf", exist_ok=True)
        
        batch_process(samples, reference, output_dir)