R 语言生信学习路线¶

# === 变量和数据类型 ===
gene_name <- "BRCA1"          # 字符串（R用 <- 赋值，也可以用 =）
gene_length <- 81189           # 数值
is_coding <- TRUE              # 逻辑值
gc_content <- 0.423            # 浮点数

# === 向量（R的基本数据结构，类似Python的列表）===
samples <- c("WT_1", "WT_2", "KO_1", "KO_2")  # 字符向量
expression <- c(12.3, 15.6, 3.2, 2.8)          # 数值向量
names(expression) <- samples                     # 给向量元素命名

# 向量化操作（R的精髓：不用写for循环！）
log2_expr <- log2(expression + 1)     # 对所有元素同时做log2转换
high_expr <- expression[expression > 10]  # 筛选大于10的元素

# === 数据框（Data Frame，相当于Python的pandas DataFrame）===
gene_df <- data.frame(
  gene_id = c("BRCA1", "TP53", "EGFR"),      # 基因ID列
  length = c(81189, 19149, 188307),            # 长度列
  chr = c("chr17", "chr17", "chr7"),           # 染色体列
  is_oncogene = c(FALSE, FALSE, TRUE),         # 是否为癌基因
  stringsAsFactors = FALSE                     # 字符串不自动转因子
)

# 数据框操作
gene_df$gc_content <- c(0.42, 0.38, 0.55)     # 添加新列
gene_df[gene_df$chr == "chr17", ]              # 筛选chr17的行
nrow(gene_df)                                   # 行数
ncol(gene_df)                                   # 列数

# === 函数 ===
calc_gc <- function(sequence) {
  # 计算DNA序列的GC含量
  seq_upper <- toupper(sequence)               # 转大写
  gc <- nchar(gsub("[^GCgc]", "", seq_upper))  # 数G和C的个数
  total <- nchar(seq_upper)                     # 总长度
  return(gc / total)                            # 返回GC比例
}

calc_gc("ATCGATCGGGCC")  # 结果：0.5833

# === 读取表格文件 ===
# 读取TSV文件（tab分隔）
expr_df <- read.delim(
  "gene_expression.tsv",    # 文件路径
  header = TRUE,             # 第一行是表头
  row.names = 1,             # 第一列作为行名
  check.names = FALSE        # 不自动修改列名中的特殊字符
)

# 读取CSV文件
metadata <- read.csv("sample_info.csv", header = TRUE)

# === 写入文件 ===
write.csv(gene_df, "output.csv", row.names = FALSE)        # 写CSV
write.table(gene_df, "output.tsv", sep = "\t", quote = FALSE)  # 写TSV

# === 读取大文件（用data.table，比read.csv快10倍）===
library(data.table)                        # 加载data.table包
big_df <- fread("huge_matrix.tsv")         # 自动识别分隔符，极速读取

library(tidyverse)  # 一次性加载 dplyr, ggplot2, tidyr 等

# === 读取示例数据 ===
expr_df <- read_tsv("gene_expression.tsv")  # tidyverse 风格的读取函数

# === 五大动词（面试必背）===

# 1. filter() —— 筛选行（白话：挑出我要的行）
sig_genes <- expr_df %>%                    # %>% 管道符，把左边的结果传给右边
  filter(padj < 0.05, abs(log2FC) > 1)      # 筛选：调整p值<0.05 且 |FC|>1

# 2. select() —— 选择列（白话：只要这几列）
subset_df <- expr_df %>%
  select(gene_id, log2FC, padj, gene_name)  # 只保留这4列

# 3. mutate() —— 新增/修改列（白话：加一列新数据）
expr_df <- expr_df %>%
  mutate(
    significance = case_when(                # 类似Excel的IF嵌套
      padj < 0.05 & log2FC > 1  ~ "Up",     # 上调
      padj < 0.05 & log2FC < -1 ~ "Down",   # 下调
      TRUE                      ~ "NS"      # 不显著
    ),
    neg_log10_pval = -log10(padj)            # 计算-log10(p值)
  )

# 4. arrange() —— 排序（白话：按某列排个序）
top_genes <- expr_df %>%
  arrange(padj) %>%                          # 按p值从小到大排
  head(20)                                   # 取前20个

# 5. summarise() + group_by() —— 分组统计
summary_stats <- expr_df %>%
  group_by(significance) %>%                 # 按显著性分组
  summarise(
    n_genes = n(),                           # 每组基因数
    mean_fc = mean(log2FC),                  # 每组平均FC
    median_pval = median(padj)               # 每组中位p值
  )

# === 宽表转长表（wide -> long）===
# 生信中经常需要：表达矩阵（宽表）-> ggplot画图需要的长表
expr_wide <- data.frame(
  gene = c("BRCA1", "TP53"),
  WT_1 = c(12.3, 45.1),      # 宽表：每个样本一列
  WT_2 = c(15.6, 42.3),
  KO_1 = c(3.2, 44.8),
  KO_2 = c(2.8, 43.9)
)

expr_long <- expr_wide %>%
  pivot_longer(
    cols = -gene,              # 除了gene列，其他都转成长格式
    names_to = "sample",       # 列名变成sample列的值
    values_to = "expression"   # 值变成expression列
  )
# 结果：gene | sample | expression
#        BRCA1 | WT_1   | 12.3
#        BRCA1 | WT_2   | 15.6  ... 以此类推

# === 长表转宽表（long -> wide）===
expr_back <- expr_long %>%
  pivot_wider(
    names_from = sample,       # sample列的值变成新的列名
    values_from = expression   # expression列的值填充
  )

# === 分离和合并列 ===
expr_long <- expr_long %>%
  separate(sample, into = c("group", "replicate"), sep = "_")  # WT_1 -> group=WT, replicate=1

library(ggplot2)

# === 火山图（生信面试必会）===
ggplot(expr_df, aes(x = log2FC, y = neg_log10_pval, color = significance)) +
  geom_point(alpha = 0.5, size = 1) +              # 画散点
  scale_color_manual(values = c(                     # 自定义颜色
    "Up" = "red",
    "Down" = "blue",
    "NS" = "gray60"
  )) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +  # p值阈值线
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +      # FC阈值线
  labs(
    x = "log2 Fold Change",
    y = "-log10(Adjusted P-value)",
    title = "Volcano Plot: KO vs WT"
  ) +
  theme_bw() +                                       # 白色背景主题
  theme(legend.position = "top")                      # 图例放顶部

ggsave("volcano.pdf", width = 8, height = 6, dpi = 300)  # 保存PDF

# === 箱线图（分组比较）===
ggplot(expr_long, aes(x = group, y = expression, fill = group)) +
  geom_boxplot(outlier.shape = NA) +                  # 箱线图（隐藏离群点）
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.5) + # 加抖动的散点
  scale_fill_manual(values = c("WT" = "#2196F3", "KO" = "#F44336")) +
  labs(x = "Group", y = "Expression (TPM)") +
  theme_classic() +                                    # 经典主题
  stat_compare_means(method = "t.test")                # 加统计检验p值

# === 热图（用pheatmap包，比ggplot更方便）===
library(pheatmap)
pheatmap(
  as.matrix(expr_wide[, -1]),     # 去掉基因名列，转矩阵
  scale = "row",                   # 按行标准化（Z-score）
  clustering_method = "ward.D2",   # 聚类方法
  color = colorRampPalette(c("blue", "white", "red"))(100),  # 蓝白红渐变
  show_rownames = TRUE,            # 显示行名（基因名）
  fontsize_row = 8,                # 行名字体大小
  main = "Gene Expression Heatmap" # 标题
)

# === 安装（只需一次）===
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("DESeq2")

library(DESeq2)

# === 准备数据 ===
# count_matrix: 行=基因，列=样本，值=原始counts（整数）
# col_data: 样本信息表，包含分组等

count_matrix <- as.matrix(read.delim("raw_counts.tsv", row.names = 1))
col_data <- data.frame(
  condition = factor(c("WT", "WT", "WT", "KO", "KO", "KO")),  # 分组信息
  row.names = colnames(count_matrix)                             # 样本名
)

# === 创建 DESeq2 对象 ===
dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,    # 原始计数矩阵
  colData = col_data,          # 样本信息
  design = ~ condition         # 实验设计公式：按condition分组比较
)

# === 预过滤（去掉低表达基因，减少多重检验负担）===
keep <- rowSums(counts(dds) >= 10) >= 3  # 至少3个样本中counts>=10
dds <- dds[keep, ]

# === 运行差异分析（一行搞定！）===
dds <- DESeq(dds)    # 内部做了：估计size factor -> 估计离散度 -> 统计检验

# === 提取结果 ===
res <- results(dds,
  contrast = c("condition", "KO", "WT"),  # 比较：KO vs WT
  alpha = 0.05                             # FDR阈值
)

# 查看结果摘要
summary(res)

# 转为数据框并排序
res_df <- as.data.frame(res) %>%
  rownames_to_column("gene_id") %>%        # 行名变成gene_id列
  arrange(padj) %>%                         # 按调整p值排序
  filter(!is.na(padj))                      # 去掉NA

# 显著差异基因
sig_genes <- res_df %>%
  filter(padj < 0.05, abs(log2FoldChange) > 1)
cat("显著差异基因数:", nrow(sig_genes), "\n")
cat("上调:", sum(sig_genes$log2FoldChange > 0), "\n")
cat("下调:", sum(sig_genes$log2FoldChange < 0), "\n")

# === 可视化 ===
# MA图
plotMA(res, ylim = c(-5, 5), main = "MA Plot")

# PCA图
vsd <- vst(dds, blind = FALSE)  # 方差稳定转换（比log2更适合PCA）
plotPCA(vsd, intgroup = "condition")

library(edgeR)

# === 创建对象 ===
group <- factor(c("WT", "WT", "WT", "KO", "KO", "KO"))
dge <- DGEList(counts = count_matrix, group = group)  # 创建edgeR对象

# === 过滤低表达基因 ===
keep <- filterByExpr(dge)        # edgeR自带的过滤函数
dge <- dge[keep, , keep.lib.sizes = FALSE]

# === 标准化 ===
dge <- calcNormFactors(dge)      # TMM标准化（edgeR默认方法）

# === 估计离散度 ===
design <- model.matrix(~ group)  # 设计矩阵
dge <- estimateDisp(dge, design) # 估计离散度

# === 差异检验 ===
fit <- glmQLFit(dge, design)     # 拟合广义线性模型
qlf <- glmQLFTest(fit)           # 似然比检验

# === 提取结果 ===
topTags(qlf, n = 20)             # 前20个差异基因

library(clusterProfiler)
library(org.Hs.eg.db)   # 人类基因注释数据库

# === 基因ID转换（通常要把gene symbol转成ENTREZ ID）===
gene_list <- sig_genes$gene_id   # 显著差异基因列表
gene_entrez <- bitr(
  gene_list,
  fromType = "SYMBOL",           # 从基因symbol
  toType = "ENTREZID",           # 转成ENTREZ ID
  OrgDb = org.Hs.eg.db           # 用人类注释库
)

# === GO富集分析 ===
go_result <- enrichGO(
  gene = gene_entrez$ENTREZID,   # 差异基因的ENTREZ ID
  OrgDb = org.Hs.eg.db,          # 注释数据库
  ont = "BP",                     # BP=生物过程, MF=分子功能, CC=细胞组分
  pAdjustMethod = "BH",          # BH校正（即FDR）
  pvalueCutoff = 0.05,           # p值阈值
  qvalueCutoff = 0.05            # q值阈值
)

# 可视化
dotplot(go_result, showCategory = 20)   # 气泡图（最常用）
barplot(go_result, showCategory = 15)   # 条形图
cnetplot(go_result, showCategory = 5)   # 基因-通路网络图

# === KEGG通路富集 ===
kegg_result <- enrichKEGG(
  gene = gene_entrez$ENTREZID,
  organism = "hsa",              # hsa = Homo sapiens（人类）
  pvalueCutoff = 0.05
)

dotplot(kegg_result, showCategory = 15)

# === GSEA分析（基因集富集分析）===
# 需要全部基因的排序列表（不是只取显著的！）
gene_rank <- res_df$log2FoldChange       # 用log2FC作为排序指标
names(gene_rank) <- res_df$gene_id       # 以基因名命名
gene_rank <- sort(gene_rank, decreasing = TRUE)  # 从大到小排序

gsea_result <- gseGO(
  geneList = gene_rank,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pvalueCutoff = 0.05
)
gseaplot2(gsea_result, geneSetID = 1)    # 画GSEA图

library(ComplexHeatmap)
library(circlize)   # 颜色映射

# === 高级热图 ===
# 准备颜色
col_fun <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))  # Z-score颜色映射

# 准备注释
ha_top <- HeatmapAnnotation(
  Group = col_data$condition,                         # 样本分组注释
  col = list(Group = c("WT" = "#2196F3", "KO" = "#F44336"))  # 分组颜色
)

# 画热图
Heatmap(
  as.matrix(sig_genes_expr),    # 表达矩阵
  name = "Z-score",              # 图例名
  col = col_fun,                 # 颜色方案
  top_annotation = ha_top,       # 顶部注释
  show_row_names = FALSE,        # 基因太多就不显示行名
  cluster_columns = TRUE,        # 列聚类
  cluster_rows = TRUE,           # 行聚类
  row_km = 3                     # 把基因分成3个cluster
)