Python 生信学习路线¶

# === 变量和数据类型 ===
gene_name = "BRCA1"           # 字符串：基因名
gene_length = 81189           # 整数：基因长度（碱基对）
gc_content = 0.423            # 浮点数：GC含量
is_coding = True              # 布尔值：是否编码蛋白

# === 列表（最常用的数据结构）===
samples = ["WT_1", "WT_2", "KO_1", "KO_2"]  # 样本列表
expression = [12.3, 15.6, 3.2, 2.8]          # 表达量列表

# 列表推导式（生信中超常用，一行代替 for 循环）
# 白话：从 expression 里挑出大于 10 的值
high_expr = [x for x in expression if x > 10]  # 结果：[12.3, 15.6]

# === 字典（存储键值对，像查字典一样查数据）===
gene_info = {
    "name": "BRCA1",           # 基因名
    "chr": "chr17",            # 染色体位置
    "length": 81189,           # 长度
    "function": "DNA repair"   # 功能
}
print(gene_info["chr"])        # 输出：chr17

# === 函数（把重复的操作封装起来）===
def calc_gc_content(sequence):
    """计算DNA序列的GC含量"""
    sequence = sequence.upper()                    # 统一转大写
    gc_count = sequence.count("G") + sequence.count("C")  # 数G和C的个数
    return gc_count / len(sequence)                # GC个数除以总长度

# 使用函数
seq = "ATCGATCGGGCC"
print(f"GC含量: {calc_gc_content(seq):.2%}")  # 输出：GC含量: 58.33%

# === 读取 FASTA 文件（最基础的序列文件格式）===
def read_fasta(filepath):
    """读取FASTA文件，返回 {序列名: 序列} 的字典"""
    sequences = {}            # 空字典，准备存结果
    current_name = ""         # 当前序列名
    current_seq = []          # 当前序列的片段列表

    with open(filepath, "r") as f:          # 打开文件
        for line in f:                       # 逐行读取
            line = line.strip()              # 去掉行尾换行符
            if line.startswith(">"):         # 以 > 开头 = 序列名行
                if current_name:             # 如果之前有序列，先存起来
                    sequences[current_name] = "".join(current_seq)
                current_name = line[1:]      # 去掉 > 号，取序列名
                current_seq = []             # 重置序列片段
            else:
                current_seq.append(line)     # 拼接序列片段
        if current_name:                     # 存最后一条序列
            sequences[current_name] = "".join(current_seq)

    return sequences

# === 读取 TSV/CSV 表格（基因表达矩阵等）===
# 方法1：纯 Python（不依赖第三方库）
with open("gene_counts.tsv", "r") as f:
    header = f.readline().strip().split("\t")   # 读表头
    for line in f:                               # 逐行读数据
        fields = line.strip().split("\t")        # 按tab分割
        gene = fields[0]                         # 第一列是基因名
        counts = [int(x) for x in fields[1:]]   # 后面的列是计数

import numpy as np  # 导入numpy，约定俗成缩写为np

# === 创建数组（比列表快100倍的数值容器）===
expr_matrix = np.array([
    [12.3, 15.6, 3.2, 2.8],   # 基因1在4个样本的表达量
    [45.1, 42.3, 44.8, 43.9],  # 基因2
    [0.5,  0.3,  28.7, 31.2],  # 基因3
])

# === 常用操作 ===
print(expr_matrix.shape)       # 查看矩阵大小：(3, 4) = 3个基因x4个样本
print(expr_matrix.mean(axis=1))  # 每个基因的平均表达量（axis=1 按行算）
print(expr_matrix.mean(axis=0))  # 每个样本的平均表达量（axis=0 按列算）

# log2转换（基因表达分析的标准操作）
log2_expr = np.log2(expr_matrix + 1)  # +1 是为了避免 log2(0) 报错

# 标准化（Z-score，让不同基因可以比较）
mean_per_gene = expr_matrix.mean(axis=1, keepdims=True)  # 每基因均值
std_per_gene = expr_matrix.std(axis=1, keepdims=True)    # 每基因标准差
z_scores = (expr_matrix - mean_per_gene) / std_per_gene  # Z-score标准化

import pandas as pd  # 导入pandas，约定俗成缩写为pd

# === 读取基因表达矩阵 ===
# 假设文件格式：第一列基因名，后面是各样本的表达量
df = pd.read_csv("gene_expression.tsv", sep="\t", index_col=0)
# sep="\t" 表示tab分隔
# index_col=0 表示第一列作为行名（基因名）

# === 基本查看 ===
print(df.shape)        # 多少基因 x 多少样本
print(df.head())       # 看前5行
print(df.describe())   # 基本统计：均值、标准差、最大最小等

# === 数据筛选（面试高频）===
# 筛选表达量大于10的基因（所有样本的均值 > 10）
high_expr_genes = df[df.mean(axis=1) > 10]  # axis=1 按行（每个基因）算均值

# 筛选特定样本
wt_samples = df[["WT_1", "WT_2"]]   # 只取WT组的列
ko_samples = df[["KO_1", "KO_2"]]   # 只取KO组的列

# === 计算差异倍数（Fold Change）===
wt_mean = df[["WT_1", "WT_2"]].mean(axis=1)   # WT组均值
ko_mean = df[["KO_1", "KO_2"]].mean(axis=1)   # KO组均值
log2fc = np.log2(ko_mean / wt_mean)            # log2 Fold Change

# === 合并数据 ===
# 把基因注释信息和表达矩阵合并
annotation = pd.read_csv("gene_annotation.tsv", sep="\t")  # 读注释文件
merged = df.merge(annotation, left_index=True, right_on="gene_id")  # 按基因ID合并

# === 分组统计 ===
# 按染色体分组，统计每条染色体有多少基因
gene_per_chr = annotation.groupby("chromosome").size()  # 按chromosome列分组计数

import matplotlib.pyplot as plt  # 画图基础库
import seaborn as sns            # 更好看的统计图库

# === 火山图（Volcano Plot，差异分析必画）===
fig, ax = plt.subplots(figsize=(8, 6))  # 创建画布，8x6英寸
ax.scatter(
    log2fc,                              # x轴：log2 Fold Change
    -np.log10(pvalues),                  # y轴：-log10(p-value)
    c=["red" if abs(fc) > 1 and p < 0.05 else "gray"  # 显著的标红
        for fc, p in zip(log2fc, pvalues)],
    alpha=0.5,                           # 透明度
    s=5                                  # 点的大小
)
ax.set_xlabel("log2 Fold Change")        # x轴标签
ax.set_ylabel("-log10(P-value)")         # y轴标签
ax.axhline(-np.log10(0.05), ls="--", color="gray")  # p=0.05的阈值线
ax.axvline(-1, ls="--", color="gray")    # FC=-2的阈值线
ax.axvline(1, ls="--", color="gray")     # FC=2的阈值线
plt.tight_layout()                       # 自动调整布局
plt.savefig("volcano_plot.pdf", dpi=300) # 保存为PDF，300dpi
plt.close()                              # 关闭画布释放内存

# === 热图（Heatmap，聚类分析必画）===
sns.clustermap(
    df.head(50),               # 取前50个基因画
    cmap="RdBu_r",             # 红蓝配色（经典）
    z_score=0,                 # 按行做Z-score标准化
    figsize=(10, 12),          # 图片大小
    method="ward",             # 聚类方法
    metric="euclidean"         # 距离度量
)
plt.savefig("heatmap.pdf", dpi=300)
plt.close()

# === 箱线图（比较分组差异）===
fig, ax = plt.subplots(figsize=(6, 4))
data = pd.DataFrame({
    "Expression": list(wt_mean) + list(ko_mean),   # 合并两组数据
    "Group": ["WT"]*len(wt_mean) + ["KO"]*len(ko_mean)  # 分组标签
})
sns.boxplot(data=data, x="Group", y="Expression", ax=ax)
ax.set_title("Gene Expression Comparison")
plt.tight_layout()
plt.savefig("boxplot.pdf", dpi=300)
plt.close()

from Bio import SeqIO        # 序列文件读写模块
from Bio.Seq import Seq       # 序列对象
from Bio import Entrez        # NCBI数据库查询

# === 序列基本操作 ===
dna = Seq("ATGCGATCGATCG")   # 创建DNA序列对象
print(dna.complement())       # 互补链：TACGCTAGCTAGC
print(dna.reverse_complement())  # 反向互补链
print(dna.transcribe())       # 转录为RNA：AUGCGAUCGAUCG
print(dna.translate())        # 翻译为蛋白质

# === 读取 FASTA 文件（比自己写的版本更专业）===
for record in SeqIO.parse("genes.fasta", "fasta"):  # 逐条读取
    print(f"ID: {record.id}")                # 序列ID
    print(f"长度: {len(record.seq)}")         # 序列长度
    print(f"描述: {record.description}")      # 描述信息
    gc = (record.seq.count("G") + record.seq.count("C")) / len(record.seq)
    print(f"GC含量: {gc:.2%}")               # GC含量

# === 读取 FASTQ 文件（测序原始数据）===
quality_scores = []
for record in SeqIO.parse("reads.fastq", "fastq"):
    quals = record.letter_annotations["phred_quality"]  # 质量分数列表
    avg_qual = sum(quals) / len(quals)                   # 平均质量
    quality_scores.append(avg_qual)

print(f"平均质量分数: {sum(quality_scores)/len(quality_scores):.1f}")

# === 查询 NCBI 数据库 ===
Entrez.email = "your@email.com"   # NCBI要求提供邮箱
handle = Entrez.efetch(
    db="nucleotide",               # 查核酸数据库
    id="NM_007294",                # BRCA1的RefSeq ID
    rettype="fasta",               # 返回FASTA格式
    retmode="text"                 # 文本模式
)
record = SeqIO.read(handle, "fasta")  # 读取结果
print(f"BRCA1序列长度: {len(record.seq)}")
handle.close()

import pysam  # 处理BAM/SAM/VCF等文件的库

# === 读取 BAM 文件 ===
bamfile = pysam.AlignmentFile("sample.bam", "rb")  # rb = read binary（二进制读取）

# 统计比对到 chr17 区域的 reads
count = 0
for read in bamfile.fetch("chr17", 41196312, 41277500):  # BRCA1区域
    if not read.is_unmapped:     # 排除未比对的reads
        if read.mapping_quality >= 20:  # 只要比对质量>=20的
            count += 1
print(f"BRCA1区域的高质量reads数: {count}")

# === 统计比对率 ===
stats = pysam.idxstats(bamfile.filename.decode())  # 获取索引统计
total_mapped = 0       # 总比对reads数
total_unmapped = 0     # 总未比对reads数
for line in stats.strip().split("\n"):
    fields = line.split("\t")
    total_mapped += int(fields[2])      # 第3列：比对数
    total_unmapped += int(fields[3])    # 第4列：未比对数
mapping_rate = total_mapped / (total_mapped + total_unmapped)
print(f"比对率: {mapping_rate:.2%}")

bamfile.close()  # 用完要关闭

# === 读取 VCF 文件（变异信息）===
vcf = pysam.VariantFile("variants.vcf.gz")  # 读取VCF文件
for record in vcf:
    chrom = record.chrom         # 染色体
    pos = record.pos             # 位置
    ref = record.ref             # 参考碱基
    alt = record.alts[0]         # 突变碱基
    qual = record.qual           # 质量分数
    print(f"{chrom}:{pos} {ref}>{alt} (QUAL={qual})")

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score  # 数据分割和交叉验证
from sklearn.preprocessing import StandardScaler     # 标准化
from sklearn.ensemble import RandomForestClassifier  # 随机森林
from sklearn.metrics import (
    classification_report,    # 分类报告
    roc_auc_score,           # AUC分数
    confusion_matrix         # 混淆矩阵
)

# === 加载数据（以疾病分类为例）===
df = pd.read_csv("microbiome_features.tsv", sep="\t")  # 微生物组特征
X = df.drop(columns=["sample_id", "group"])  # 特征矩阵（去掉ID和标签列）
y = df["group"]  # 标签：T2D=1, Control=0

# === 数据分割 ===
X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.2,       # 20%作为测试集
    random_state=42,     # 随机种子（保证可重复）
    stratify=y           # 分层抽样（保证正负样本比例一致）
)

# === 标准化（让所有特征在同一数量级）===
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)   # 在训练集上拟合+转换
X_test_scaled = scaler.transform(X_test)          # 在测试集上只转换（不能再fit！）

# === 训练随机森林 ===
rf = RandomForestClassifier(
    n_estimators=100,     # 100棵树
    max_depth=10,         # 最大深度10
    random_state=42,      # 随机种子
    n_jobs=-1             # 用所有CPU核心
)
rf.fit(X_train_scaled, y_train)  # 训练模型

# === 评估 ===
y_pred = rf.predict(X_test_scaled)        # 预测
y_prob = rf.predict_proba(X_test_scaled)[:, 1]  # 预测概率（取正类概率）

print(classification_report(y_test, y_pred))     # 精确率、召回率、F1
print(f"AUC: {roc_auc_score(y_test, y_prob):.3f}")  # AUC分数

# === 交叉验证（更靠谱的评估方式）===
cv_scores = cross_val_score(rf, X_train_scaled, y_train, cv=5, scoring="roc_auc")
print(f"5折交叉验证 AUC: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# === 特征重要性（面试爱问）===
feature_importance = pd.Series(
    rf.feature_importances_,   # 每个特征的重要性分数
    index=X.columns            # 特征名
).sort_values(ascending=False)

print("Top 10 重要特征:")
print(feature_importance.head(10))

import torch                          # PyTorch核心
import torch.nn as nn                  # 神经网络模块
from torch.utils.data import DataLoader, TensorDataset  # 数据加载

# === 定义简单的分类网络 ===
class BioClassifier(nn.Module):
    """生信数据分类器：3层全连接网络"""
    def __init__(self, input_dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, 128),   # 输入层 -> 128个神经元
            nn.ReLU(),                    # 激活函数（给网络加非线性能力）
            nn.Dropout(0.3),              # 随机关闭30%的神经元（防过拟合）
            nn.Linear(128, 64),           # 128 -> 64
            nn.ReLU(),
            nn.Dropout(0.3),
            nn.Linear(64, 1),             # 64 -> 1（二分类输出）
            nn.Sigmoid()                  # 压缩到0-1之间（概率）
        )

    def forward(self, x):
        return self.net(x)               # 前向传播

# === 准备数据 ===
X_tensor = torch.FloatTensor(X_train_scaled)  # numpy转pytorch张量
y_tensor = torch.FloatTensor(y_train.values).unsqueeze(1)  # 标签也转张量
dataset = TensorDataset(X_tensor, y_tensor)   # 组合为数据集
loader = DataLoader(dataset, batch_size=32, shuffle=True)  # 批量加载器

# === 训练 ===
model = BioClassifier(X_train_scaled.shape[1])  # 实例化模型
criterion = nn.BCELoss()            # 二分类交叉熵损失
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)  # Adam优化器

for epoch in range(50):             # 训练50轮
    model.train()                    # 训练模式
    total_loss = 0
    for batch_X, batch_y in loader:  # 每批32个样本
        pred = model(batch_X)        # 前向传播
        loss = criterion(pred, batch_y)  # 计算损失
        optimizer.zero_grad()        # 清空梯度
        loss.backward()              # 反向传播
        optimizer.step()             # 更新参数
        total_loss += loss.item()
    if (epoch + 1) % 10 == 0:       # 每10轮打印一次
        print(f"Epoch {epoch+1}, Loss: {total_loss/len(loader):.4f}")

# Snakefile（Snakemake的流程定义文件，用Python语法写）

# 定义所有样本名
SAMPLES = ["WT_1", "WT_2", "KO_1", "KO_2"]

# 最终目标：生成多样本质控报告
rule all:
    input:
        "results/multiqc_report.html"   # 最终要生成的文件

# 步骤1：质控
rule fastqc:
    input:
        "data/{sample}.fastq.gz"        # 输入：原始测序数据
    output:
        "results/fastqc/{sample}_fastqc.html"  # 输出：质控报告
    shell:
        "fastqc {input} -o results/fastqc/"    # 运行fastqc命令

# 步骤2：汇总质控报告
rule multiqc:
    input:
        expand("results/fastqc/{sample}_fastqc.html", sample=SAMPLES)  # 所有样本的质控
    output:
        "results/multiqc_report.html"
    shell:
        "multiqc results/fastqc/ -o results/"