生物信息学实验基础

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点

一个典型的生物信息学分析项目通常遵循以下流程：
第一阶段：项目规划
1. 明确科学问题
2. 确定所需数据类型和来源
3. 选择合适的分析方法和工具
4. 评估计算资源需求
5. 制定数据管理计划
第二阶段：数据获取与预处理
1. 数据采集：从公共数据库下载（如NCBI SRA、ENA）或从合作者处获取原始数据
2. 质量控制（QC）：评估数据质量，识别并处理低质量数据。对于测序数据，使用FastQC等工具检查碱基质量分布、GC含量、接头污染、重复序列等
3. 数据清洗：去除低质量序列、接头序列、PCR重复等
4. 格式转换：将数据转换为分析所需的标准格式
第三阶段：核心分析
根据研究问题选择相应的分析方法，例如：
- 序列比对：将测序reads定位到参考基因组
- 变异检测：识别SNP、Indel等遗传变异
- 表达量计算：定量基因或转录本的表达水平
- 功能注释：将分析结果与生物学功能关联
第四阶段：结果整合与解释
1. 多来源结果整合
2. 统计显著性评估和多重检验校正
3. 功能富集分析（GO、KEGG等）
4. 可视化展示
第五阶段：报告与共享
1. 撰写分析报告
2. 整理代码和文档
3. 数据和代码归档
4. 论文发表和结果共享

掌握标准数据格式是进行生物信息学分析的基础：
| 格式 | 用途 | 说明 |
|------|------|------|
| FASTA | 序列存储 | >header开头的文本格式，用于存储DNA、RNA、蛋白质序列 |
| FASTQ | 测序数据 | FASTA+质量值，每4行一个read（header、序列、+、质量） |
| SAM/BAM | 比对结果 | SAM为文本格式，BAM为二进制压缩格式，存储reads比对信息 |
| VCF | 变异信息 | 存储SNP、Indel等变异的位置、基因型和质量信息 |
| GFF/GTF | 基因组注释 | 存储基因、外显子、CDS等特征的位置和属性信息 |
| BED | 基因组区间 | 简单的三列格式（chr, start, end），用于表示基因组区域 |
| GCT/TPM | 表达矩阵 | 存储基因表达量矩阵，行是基因，列是样本 |

FAIR原则：科学数据管理应遵循FAIR原则——可发现（Findable）、可访问（Accessible）、可互操作（Interoperable）、可重用（Reusable）。
可重复性的层次：
1. 结果可重复（Repeatability）：同一研究者使用相同数据和方法重复实验，获得相同结果
2. 可复现（Reproducibility）：不同研究者使用相同数据和方法，获得相同结果
3. 可再实现（Replicability）：在不同数据集上应用相同方法，获得一致结论
实现可重复性的最佳实践：
- 使用版本控制系统（Git）管理代码
- 记录软件名称和版本号
- 保存完整的分析参数
- 使用虚拟环境或容器隔离运行环境
- 编写清晰的文档和README文件
- 使用工作流管理系统自动化分析流程

假设我们要比较某种处理条件下（Treatment）和对照（Control）的基因表达差异，每组3个生物学重复。
项目结构：

rnaseq_project/
├── data/
│   ├── raw/              # 原始测序数据（FASTQ）
│   ├── reference/        # 参考基因组和注释文件
│   └── processed/        # 处理后的数据
├── scripts/              # 分析脚本
│   ├── 01_qc.sh
│   ├── 02_alignment.sh
│   ├── 03_quantification.sh
│   └── 04_de_analysis.R
├── results/              # 分析结果
│   ├── qc_reports/
│   ├── bam_files/
│   ├── counts/
│   └── de_results/
├── envs/                 # 环境配置文件
│   └── rnaseq_env.yml
├── README.md             # 项目说明文档
└── Snakefile             # Snakemake工作流文件

分析流程概览：

# 第一步：质量控制
fastqc data/raw/*.fastq.gz -o results/qc_reports/
multiqc results/qc_reports/ -o results/qc_reports/summary/
# 第二步：序列比对
hisat2 -p 8 -x reference/genome -1 sample_R1.fq.gz -2 sample_R2.fq.gz | \
    samtools sort -@ 4 -o results/bam_files/sample.bam
# 第三步：表达量定量
featureCounts -T 8 -a reference/genes.gtf -o results/counts/counts.txt \
    results/bam_files/*.bam
# 第四步：差异表达分析（在R中完成）
# DESeq2分析

13.2.1 Linux操作系统

变量和数据类型：

# 数字
integer = 42
floating = 3.14159
# 字符串
dna = "ATGCGCTAGCTA"
protein = 'MPLK'
# 字符串操作
print(len(dna))              # 长度
print(dna[0:3])              # 切片（前3个碱基）
print(dna.count("GC"))       # 计数
print(dna.replace("T", "U"))  # 替换（DNA转RNA）
print(dna.find("ATG"))       # 查找子串位置
# 列表（有序可变）
sequences = ["ATCG", "GCTA", "TAGC"]
sequences.append("CGAT")     # 添加元素
print(sequences[0])          # 访问
print(len(sequences))        # 长度
# 字典（键值对）
gene_expr = {
    "GAPDH": 25.3,
    "ACTB": 18.7,
    "TP53": 3.2
}
print(gene_expr["GAPDH"])    # 访问
gene_expr["BRCA1"] = 7.5     # 添加
# 元组（有序不可变）
coordinates = (100, 200)

控制流：

# 条件语句
gc_content = 0.55
if gc_content > 0.6:
    print("High GC")
elif gc_content > 0.4:
    print("Moderate GC")
else:
    print("Low GC")
# for循环
for seq in sequences:
    gc = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"{seq}: GC={gc:.2%}")
# while循环
i = 0
while i < len(sequences):
    print(sequences[i])
    i += 1
# 列表推导式（Pythonic写法）
gc_values = [(s.count("G") + s.count("C")) / len(s) for s in sequences]

函数：

def calculate_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return gc_count / len(sequence)
# 调用
gc = calculate_gc_content("ATGCGCTAGCTA")
print(f"GC content: {gc:.2%}")
# 默认参数
def reverse_complement(seq, rna=False):
    """计算反向互补序列"""
    complement = {"A": "T", "T": "A", "G": "C", "C": "G",
                  "a": "t", "t": "a", "g": "c", "c": "g"}
    if rna:
        complement["A"] = "U"
        complement["a"] = "u"
    rc = "".join(complement.get(base, base) for base in reversed(seq))
    return rc
print(reverse_complement("ATGC"))       # GCAT
print(reverse_complement("ATGC", rna=True))  # GCAU

文件操作：

# 读取文件
with open("sequences.fasta", "r") as f:
    content = f.read()           # 读取全部
    lines = f.readlines()        # 读取为列表
# 写入文件
with open("output.txt", "w") as f:
    f.write("Hello\n")
# 逐行读取（推荐大文件使用）
with open("data.txt", "r") as f:
    for line in f:
        line = line.strip()      # 去除末尾换行符
        if line.startswith(">"):
            print(f"Header: {line}")

安装：

pip install biopython
# 或
conda install -c conda-forge biopython

序列处理（Seq对象）：

from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
# 创建序列对象
dna_seq = Seq("ATGCGCTAGCTA")
# 序列属性
print(f"长度: {len(dna_seq)}")
print(f"GC含量: {GC(dna_seq):.1f}%")
print(f"分子量: {molecular_weight(dna_seq):.1f}")
# 序列操作
print(dna_seq.complement())        # 互补序列
print(dna_seq.reverse_complement()) # 反向互补
print(dna_seq.transcribe())        # DNA转RNA
print(dna_seq.translate())         # 翻译为蛋白质
# 转录和翻译
mrna = dna_seq.transcribe()
protein = mrna.translate()
print(f"Protein: {protein}")

解析FASTA/FASTQ文件：

from Bio import SeqIO
# 解析FASTA文件
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Description: {record.description}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")
    print()
# 解析FASTQ文件（测序数据）
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality'][:10]}")
# 将记录写入文件
records = []
for record in SeqIO.parse("input.fasta", "fasta"):
    if len(record.seq) > 100:  # 只保留长序列
        records.append(record)
SeqIO.write(records, "filtered.fasta", "fasta")
# 转换格式（FASTQ to FASTA）
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")

访问NCBI数据库：

from Bio import Entrez, SeqIO
# 设置邮箱（NCBI要求）
Entrez.email = "your.email@example.com"
# 搜索PubMed
handle = Entrez.esearch(db="pubmed", term="CRISPR[Title] AND 2023[PDAT]", retmax=10)
record = Entrez.read(handle)
print(f"Found {record['Count']} articles")
print(f"IDs: {record['IdList']}")
# 获取序列
handle = Entrez.efetch(db="nucleotide", id="NM_001301717", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
print(f"Sequence: {record.seq[:100]}...")

序列比对：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment
# 全局比对（Needleman-Wunsch）
alignments = pairwise2.align.globalxx("ATCG", "ATG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 局部比对（Smith-Waterman）
alignments = pairwise2.align.localxx("ATCGGCTA", "CGG")
for alignment in alignments:
    print(format_alignment(*alignment))
# 带参数的比对（匹配+1，错配-1，空位开启-2，空位延伸-1）
alignments = pairwise2.align.globalms("ATCG", "ATG", 1, -1, -2, -1)

BLAST解析：

from Bio.Blast import NCBIXML
# 解析BLAST XML结果
with open("blast_result.xml") as result_handle:
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        print(f"Hit: {alignment.title}")
        for hsp in alignment.hsps:
            print(f"  E-value: {hsp.expect}")
            print(f"  Identity: {hsp.identities}/{hsp.align_length}")
            print(f"  Query: {hsp.query[:50]}...")
            print(f"  Match: {hsp.match[:50]}...")
            print(f"  Sbjct: {hsp.sbjct[:50]}...")

变量和数据类型：

# 赋值（推荐使用 <-）
x <- 42
y <- 3.14
name <- "GeneA"
is_active <- TRUE
# 向量（最基本的数据结构）
expr_values <- c(10.5, 23.4, 5.2, 18.9, 12.1)
genes <- c("BRCA1", "TP53", "EGFR", "MYC")
# 向量操作
length(expr_values)
mean(expr_values)
sd(expr_values)
summary(expr_values)
# 数据框（类似Excel表格）
df <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR", "MYC"),
    expression = c(25.3, 18.7, 12.4, 8.9),
    group = c("treatment", "treatment", "control", "control")
)
# 访问数据
df$gene                    # 提取列
df$expression
df[1, ]                    # 第一行
df[, "expression"]         # 按名列
df[df$group == "treatment", ]  # 条件筛选
# 列表（可包含不同类型）
my_list <- list(
    name = "sample1",
    counts = c(100, 200, 300),
    metadata = data.frame(key = c("A", "B"), val = c(1, 2))
)
my_list$name
my_list[["counts"]]

控制流：

# 条件语句
x <- 15
if (x > 10) {
    print("x is greater than 10")
} else if (x > 5) {
    print("x is between 5 and 10")
} else {
    print("x is 5 or less")
}
# ifelse向量化条件
scores <- c(85, 92, 78, 65, 88)
grades <- ifelse(scores >= 90, "A",
                 ifelse(scores >= 80, "B",
                        ifelse(scores >= 70, "C", "D")))
# for循环
for (gene in genes) {
    print(paste("Processing:", gene))
}
# apply族函数（向量化操作，避免显式循环）
# apply用于矩阵/数组
# lapply用于列表，返回列表
# sapply用于列表，返回向量/矩阵
# tapply用于分组计算
# 示例：对数据框的数值列计算均值
numeric_cols <- sapply(df, is.numeric)
lapply(df[, numeric_cols], mean)

函数：

# 定义函数
calculate_fold_change <- function(treatment, control) {
    """计算差异倍数（log2 fold change）"""
    fc <- treatment - control  # 假设已经是log2转换的值
    return(fc)
}
# 使用
fc <- calculate_fold_change(10.5, 8.2)
print(paste("Fold change:", round(fc, 2)))
# 默认参数
normalize <- function(values, method = "zscore") {
    if (method == "zscore") {
        return((values - mean(values)) / sd(values))
    } else if (method == "minmax") {
        return((values - min(values)) / (max(values) - min(values)))
    } else {
        stop("Unknown normalization method")
    }
}

library(ggplot2)
# 创建示例数据
data <- data.frame(
    gene = rep(c("GeneA", "GeneB", "GeneC"), each = 10),
    expression = c(rnorm(10, 10, 2), rnorm(10, 15, 3), rnorm(10, 8, 1)),
    group = rep(c("Control", "Treatment"), each = 5, times = 3)
)
# 散点图
p1 <- ggplot(data, aes(x = group, y = expression, color = group)) +
    geom_point(position = position_jitter(width = 0.2), size = 3) +
    geom_boxplot(alpha = 0.3, outlier.shape = NA) +
    facet_wrap(~gene) +
    labs(title = "Gene Expression Comparison",
         x = "Condition",
         y = "Expression Level") +
    theme_minimal()
print(p1)
# 热图（使用pheatmap包）
library(pheatmap)
expr_matrix <- matrix(rnorm(100), nrow = 10)
rownames(expr_matrix) <- paste0("Gene", 1:10)
colnames(expr_matrix) <- paste0("Sample", 1:10)
pheatmap(expr_matrix, 
         scale = "row",
         clustering_method = "ward.D2",
         color = colorRampPalette(c("navy", "white", "firebrick"))(50))
# 火山图（差异表达结果）
de_results <- data.frame(
    gene = paste0("Gene", 1:1000),
    log2FC = rnorm(1000, 0, 2),
    pvalue = runif(1000)
)
de_results$padj <- p.adjust(de_results$pvalue, method = "BH")
ggplot(de_results, aes(x = log2FC, y = -log10(padj))) +
    geom_point(aes(color = abs(log2FC) > 1 & padj < 0.05), alpha = 0.5) +
    scale_color_manual(values = c("grey", "red")) +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
    labs(x = "log2 Fold Change", y = "-log10 adjusted p-value") +
    theme_bw()

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
# 安装Bioconductor包
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
BiocManager::install("Biostrings")
# Biostrings：序列处理
library(Biostrings)
dna <- DNAString("ATGCGCTAGCTA")
complement(dna)
reverseComplement(dna)
translate(dna)
# 计算GC含量
letterFrequency(dna, "GC") / length(dna)
# 读取FASTA文件
fasta_file <- readDNAStringSet("sequences.fasta")
width(fasta_file)  # 序列长度
alphabetFrequency(fasta_file)  # 碱基频率
# DESeq2：差异表达分析（简要示例）
library(DESeq2)
# countData: 基因计数矩阵（行是基因，列是样本）
# colData: 样本信息
# dds <- DESeqDataSetFromMatrix(countData = counts, 
#                               colData = sample_info,
#                               design = ~ condition)
# dds <- DESeq(dds)
# results <- results(dds)

核心概念：
- 环境（Environment）：独立的软件安装空间，不同环境之间互不干扰
- 通道（Channel）：软件包的来源仓库。Bioconda是生物信息学专用通道
- 环境文件（environment.yml）：记录环境中所有软件及其版本的配置文件
基本操作：

# 查看现有环境
conda env list
# 创建新环境
conda create -n rnaseq python=3.10
# 激活环境
conda activate rnaseq
# 安装软件
conda install -c bioconda star
conda install -c bioconda samtools
conda install -c bioconda featurecounts
conda install -c bioconda fastqc
conda install -c bioconda multiqc
# 一次性安装多个包
conda install -c bioconda star samtools featurecounts fastqc multiqc
# 查看已安装的包
conda list
# 导出环境配置
conda env export > environment.yml
# 从配置文件创建环境
conda env create -f environment.yml
# 删除环境
conda remove -n rnaseq --all
# 退出当前环境
conda deactivate

示例environment.yml文件：

name: rnaseq
channels:
  - bioconda
  - conda-forge
  - defaults
dependencies:
  - python=3.10
  - star=2.7.10b
  - samtools=1.16
  - featurecounts=2.0
  - fastqc=0.11.9
  - multiqc=1.14
  - trimmomatic=0.39
  - picard=2.27
  - r-base=4.2
  - bioconductor-deseq2=1.38
  - pip
  - pip:
    - salmon==1.9.0

基本语法：
Snakemake的核心是定义规则（rule），每个规则包含输入（input）、输出（output）和命令（shell/run/script）。

# Snakefile
# 定义样本列表
SAMPLES = ["sample1", "sample2", "sample3"]
# 目标规则：定义最终需要生成的文件
rule all:
    input:
        "results/counts_matrix.txt",
        "results/multiqc_report.html"
# 质量控制规则
rule fastqc:
    input:
        "data/{sample}_R1.fastq.gz",
        "data/{sample}_R2.fastq.gz"
    output:
        html_r1="qc/{sample}_R1_fastqc.html",
        zip_r1="qc/{sample}_R1_fastqc.zip",
        html_r2="qc/{sample}_R2_fastqc.html",
        zip_r2="qc/{sample}_R2_fastqc.zip"
    params:
        outdir="qc"
    threads: 2
    shell:
        "fastqc -t {threads} -o {params.outdir} {input}"
# 序列比对规则
rule star_align:
    input:
        r1="data/{sample}_R1.fastq.gz",
        r2="data/{sample}_R2.fastq.gz",
        index="reference/STAR_index/Genome"
    output:
        bam="aligned/{sample}.sorted.bam",
        bai="aligned/{sample}.sorted.bam.bai"
    params:
        prefix="aligned/{sample}_",
        index_dir="reference/STAR_index"
    threads: 8
    shell:
        """
        STAR --runThreadN {threads} \
             --genomeDir {params.index_dir} \
             --readFilesIn {input.r1} {input.r2} \
             --readFilesCommand zcat \
             --outFileNamePrefix {params.prefix} \
             --outSAMtype BAM SortedByCoordinate
        mv {params.prefix}Aligned.sortedByCoord.out.bam {output.bam}
        samtools index {output.bam}
        """
# 表达量定量
rule featurecounts:
    input:
        bams=expand("aligned/{sample}.sorted.bam", sample=SAMPLES),
        gtf="reference/genes.gtf"
    output:
        counts="results/counts_matrix.txt"
    threads: 4
    shell:
        "featureCounts -T {threads} -a {input.gtf} -o {output.counts} {input.bams}"
# MultiQC汇总
rule multiqc:
    input:
        expand("qc/{sample}_{read}_fastqc.zip", sample=SAMPLES, read=["R1", "R2"])
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc qc/ -o results/"

运行Snakemake：

# 预览（dry-run）
snakemake -n
# 本地运行（使用4个核）
snakemake --cores 4
# 强制重新运行
snakemake --cores 4 --forceall
# 只运行特定规则
snakemake --cores 4 featurecounts
# 使用Conda环境（自动为每个规则创建环境）
snakemake --cores 4 --use-conda
# 集群提交（SLURM）
snakemake --cluster "sbatch --time={resources.time} --mem={resources.mem} --cpus-per-task={threads}" \
          --jobs 10

Nextflow采用数据流编程模型，更适合复杂的并行计算：

// main.nf
params.reads = "data/*_{R1,R2}.fastq.gz"
params.genome = "reference/genome.fa"
params.gtf = "reference/genes.gtf"
// 定义进程
process FASTQC {
    tag "$sample_id"
    input:
    tuple val(sample_id), path(reads)
    output:
    path "fastqc_*"
    script:
    """
    fastqc -t 2 -o . ${reads}
    """
}
process STAR_ALIGN {
    tag "$sample_id"
    cpus 8
    input:
    tuple val(sample_id), path(reads)
    path index
    output:
    tuple val(sample_id), path("*.sorted.bam")
    script:
    """
    STAR --runThreadN 8 \
         --genomeDir $index \
         --readFilesIn ${reads[0]} ${reads[1]} \
         --readFilesCommand zcat \
         --outSAMtype BAM SortedByCoordinate \
         --outFileNamePrefix ${sample_id}_
    mv ${sample_id}_Aligned.sortedByCoord.out.bam ${sample_id}.sorted.bam
    """
}
// 工作流
workflow {
    Channel
        .fromFilePairs(params.reads, checkIfExists: true)
        .set { read_pairs }
    FASTQC(read_pairs)
    STAR_ALIGN(read_pairs, params.genome)
}

FastQC从多个维度评估测序数据质量：
- 基本统计：总reads数、序列长度、GC含量
- 每个碱基的质量：箱线图展示read各位置的质量分布
- 每个序列的质量：所有reads平均质量的分布
- 每个碱基的序列含量：各位置A/T/G/C的比例（应接近一致）
- GC含量：样本GC分布与理论正态分布的比较
- 接头序列：检测已知接头序列的污染
- 序列重复：评估PCR重复和过度表达的序列

BWA-MEM（Burrows-Wheeler Aligner's Maximal Exact Matches）是目前最常用的DNA序列比对算法：
1. 索引构建：将参考基因组转换为FM-index（基于Burrows-Wheeler变换），支持快速精确匹配查找。
2. 种子搜索：在参考基因组中查找read的精确匹配子串（maximal exact matches）。
3. 链式扩展：将种子连接起来，构建可能的比对位置。
4. 局部比对：使用Smith-Waterman算法对候选位置进行精细比对，处理gap（插入缺失）。
5. 比对质量计算：根据最佳和次佳比对的得分差异，计算MAPQ（mapping quality）。

SAM文件包含11个必需字段：
| 列 | 字段名 | 说明 |
|----|--------|------|
| 1 | QNAME | Read名称 |
| 2 | FLAG | 比对标志（位掩码） |
| 3 | RNAME | 参考序列名称 |
| 4 | POS | 比对起始位置（1-based） |
| 5 | MAPQ | 比对质量（Phred标度） |
| 6 | CIGAR | 比对运算字符串 |
| 7 | RNEXT | 配对read的参考序列 |
| 8 | PNEXT | 配对read的位置 |
| 9 | TLEN | 模板长度 |
| 10 | SEQ | 序列 |
| 11 | QUAL | 质量值 |
FLAG字段解读（关键值）：
| FLAG | 含义 |
|------|------|
| 1 | Read有配对 |
| 2 | Read在配对中正确比对 |
| 4 | Read未比对 |
| 8 | 配对read未比对 |
| 16 | Read比对到反向互补链 |
| 64 | 是第一个read |
| 128 | 是第二个read |
| 256 | 非主要比对 |
| 512 | 未通过QC |
| 1024 | PCR重复或光学重复 |

GATK（Genome Analysis Toolkit）是Broad Institute开发的行业标准变异检测软件包，其核心流程包括：
1. 标记重复（MarkDuplicates）：识别PCR重复reads，只保留一个拷贝用于变异检测。
2. 碱基质量校正（Base Quality Score Recalibration, BQSR）：系统性的测序误差会导致碱基质量值不准确。BQSR通过比较观察到的错配率和期望的错配率，构建校正模型，输出更准确的碱基质量值。
3. 变异检测（HaplotypeCaller）：GATK的核心变异检测引擎：
- 在感兴趣区域进行局部de novo组装
- 识别潜在的单倍型
- 使用PairHMM算法将reads与单倍型比对
- 输出基因型似然值和变异质量
4. 变异质量校正（Variant Quality Score Recalibration, VQSR）：利用已知变异位点（如dbSNP、HapMap）训练高斯混合模型，区分真实的变异和假阳性。

# === 第一步：FastQC质控 ===
mkdir -p qc_results
# 对单个样本运行FastQC
fastqc -t 4 -o qc_results sample_R1.fastq.gz sample_R2.fastq.gz
# 批量处理
for r1 in raw/*_R1.fastq.gz; do
    r2="${r1/_R1/_R2}"
    fastqc -t 4 -o qc_results "$r1" "$r2"
done
# MultiQC汇总所有质控报告
multiqc qc_results/ -o qc_results/summary/
# === 第二步：Trimmomatic质量修剪 ===
# 假设FastQC显示：
# 1. 3'端质量下降（需要SLIDINGWINDOW修剪）
# 2. 存在adapter污染（需要ILLUMINACLIP）
mkdir -p trimmed
# PE模式（双端）
trimmomatic PE -threads 8 \
    raw/sample_R1.fastq.gz raw/sample_R2.fastq.gz \
    trimmed/sample_R1_paired.fq.gz trimmed/sample_R1_unpaired.fq.gz \
    trimmed/sample_R2_paired.fq.gz trimmed/sample_R2_unpaired.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
# 参数说明：
# ILLUMINACLIP:adapters.fa:2:30:10
#   adapters.fa = adapter序列文件
#   2 = seed匹配时允许的最大错配数
#   30 = palindrome模式下的阈值
#   10 = simple模式下的阈值
# LEADING:3 - 从read起始切除质量值<3的碱基
# TRAILING:3 - 从read末尾切除质量值<3的碱基
# SLIDINGWINDOW:4:15 - 4bp窗口平均质量<15时切除
# MINLEN:36 - 丢弃长度<36bp的read

# === 第一步：建立BWA索引 ===
bwa index reference.fa
# 或使用samtools建立FASTA索引（用于IGV等工具）
samtools faidx reference.fa
# === 第二步：序列比对 ===
# BWA-MEM适合70bp-1Mbp的reads
bwa mem -t 16 \
    -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa \
    sample_R1_paired.fq.gz sample_R2_paired.fq.gz \
    > sample.sam
# -t 16: 使用16个线程
# -R: 添加read group信息（GATK必需）
# === 第三步：SAM转BAM、排序、索引 ===
# 方法1：管道一步完成
bwa mem -t 16 -R "@RG\tID:sample1\tSM:sample1\tLB:lib1\tPL:ILLUMINA" \
    reference.fa sample_R1.fq.gz sample_R2.fq.gz | \
    samtools sort -@ 4 -o sample.sorted.bam -
samtools index sample.sorted.bam
# 方法2：分步处理
samtools view -bS sample.sam > sample.bam          # SAM转BAM
samtools sort sample.bam -o sample.sorted.bam      # 排序
samtools index sample.sorted.bam                    # 建立索引
rm sample.sam sample.bam                           # 删除中间文件
# === 第四步：比对质量统计 ===
samtools flagstat sample.sorted.bam > sample.flagstat
samtools idxstats sample.sorted.bam > sample.idxstats
# === 常用SAMtools命令 ===
# 查看BAM头部
samtools view -H sample.sorted.bam
# 查看特定区域的比对
samtools view sample.sorted.bam chr1:1000000-2000000 | head
# 提取特定FLAG的reads（例如：只提取properly paired reads）
samtools view -f 2 -b sample.sorted.bam > sample.proper_pair.bam
# 过滤掉未比对reads（FLAG 4）
samtools view -F 4 -b sample.sorted.bam > sample.mapped_only.bam
# BAM转FASTQ（用于重新比对）
samtools bam2fq sample.sorted.bam > sample.fastq
# 合并多个BAM
samtools merge merged.bam sample1.sorted.bam sample2.sorted.bam

# === GATK最佳实践流程 ===
# 1. 标记重复（MarkDuplicates）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.dedup.bam \
    -M sample.metrics.txt \
    --REMOVE_DUPLICATES false  # 标记但不删除，GATK推荐
samtools index sample.dedup.bam
# 2. 碱基质量校正（BQSR）
# 需要已知变异位点（如dbSNP）作为训练集
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R reference.fa \
    --known-sites dbsnp.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R reference.fa \
    -I sample.dedup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.recal.bam
# 3. 变异检测（HaplotypeCaller）
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample.recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF  # 输出gVCF格式（适合多样本联合分析）
# 4. 多样本联合基因型（如果有多个样本）
# GenomicsDBImport导入gVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path cohort_db \
    -L intervals.list
# GenotypeGVCFs输出最终VCF
gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://cohort_db \
    -O cohort.vcf.gz
# 5. 变异质控（VQSR，可选）
gatk VariantRecalibrator \
    -R reference.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum \
    -mode SNP \
    -O cohort.recal \
    --tranches-file cohort.tranches
gatk ApplyVQSR \
    -R reference.fa \
    -V cohort.vcf.gz \
    --recal-file cohort.recal \
    --tranches-file cohort.tranches \
    --truth-sensitivity-filter-level 99.5 \
    -mode SNP \
    -O cohort.vqsr.vcf.gz
# 6. 硬过滤（如果VQSR不适用，如小样本）
gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0" \
    --filter-name "basic_filters"

# IGV是一个Java桌面应用程序，无需命令行操作
# 但它可以配合命令行工具准备输入文件
# 1. 生成IGV兼容的TDF文件（用于大BAM文件的快速浏览）
igvtools count sample.sorted.bam sample.tdf reference.fa
# 2. 生成BEDGRAPH覆盖度文件
bedtools genomecov -ibam sample.sorted.bam -bg > sample.coverage.bedgraph
# 3. 生成VCF索引（tabix）
bgzip cohort.vcf.gz
tabix -p vcf cohort.vcf.gz

IGV使用步骤：
1. 启动IGV（需要Java环境）
2. 加载参考基因组：Genomes → Load Genome from File → 选择reference.fa
3. 加载比对文件：File → Load from File → 选择sample.sorted.bam
4. 加载变异文件：File → Load from File → 选择cohort.vcf.gz
5. 导航到感兴趣的基因区域（如输入 TP53 或坐标 chr17:7,661,778-7,687,538）
6. 观察reads的比对情况、覆盖深度和变异位点