实验1：R基础编程与生物医学数据分析

掌握R语言核心语法与生信数据处理

实验目标

完成本实验后，你将能够：

  1. 熟练创建和操作R的各种数据结构
  2. 使用向量化操作处理基因表达数据
  3. 编写自定义函数封装生信分析流程
  4. 进行基础的数据处理与统计分析
  5. 安装和管理Bioconductor生信包

实验时长：90分钟
实验类型：个人完成

实验准备

环境要求

  • R >= 4.3.0
  • RStudio >= 2023.06

创建新项目

  1. 打开RStudio
  2. 菜单：File → New Project → New Directory → New Project
  3. 目录名：lab1-r-basics
  4. 勾选 Create a git repository（如果有Git）

创建脚本文件

在RStudio中，创建新脚本：File → New File → R Script，保存为 lab1.R

第一部分：数据类型与结构 (25分钟)

练习1.1：基础类型与生物医学应用

在R脚本中完成以下操作：

# TODO: 创建以下变量（结合生信场景）
# 1. 一个数值型变量，存储TP53基因的表达量 8.5，命名为 tp53_expr
# 2. 一个整型变量，存储测序读数 10000，命名为 read_count  
# 3. 一个字符型变量，存储基因名 "BRCA1"，命名为 gene_name
# 4. 一个逻辑型变量，表示是否差异表达 TRUE，命名为 is_de

# 你的代码：





# 使用 class() 函数检查每个变量的类型

提示：使用 L 后缀创建整数，如 10000L

点击查看参考答案
代码 
# 创建变量
 tp53_expr <- 8.5
read_count <- 10000L
gene_name <- "BRCA1"
is_de <- TRUE

# 检查类型
cat("tp53_expr 类型:", class(tp53_expr), "- 基因表达量\n")
tp53_expr 类型: numeric - 基因表达量
代码 
cat("read_count 类型:", class(read_count), "- 测序读数\n")
read_count 类型: integer - 测序读数
代码 
cat("gene_name 类型:", class(gene_name), "- 基因名\n")
gene_name 类型: character - 基因名
代码 
cat("is_de 类型:", class(is_de), "- 差异表达标志\n")
is_de 类型: logical - 差异表达标志

练习1.2：基因表达向量操作

# TODO: 完成以下基因表达向量操作

# 1. 创建TP53基因在10个样本中的表达值向量
#    值：12.5, 15.2, 8.7, 18.3, 11.9, 6.2, 7.5, 5.8, 8.1, 6.9
#    命名为 tp53_expression
tp53_expression <- ____

# 2. 创建样本名称向量：Tumor_1到Tumor_5，Normal_1到Normal_5
sample_names <- ____

# 3. 计算TP53的平均表达量和标准差
mean_expr <- ____
sd_expr <- ____

# 4. 提取高表达样本（表达量 > 10）
high_expr <- ____

# 5. 计算高表达样本的数量
n_high <- ____

# 6. 找出表达量最高的样本位置
max_pos <- ____
点击查看参考答案
代码 
# 1. 创建表达向量
tp53_expression <- c(12.5, 15.2, 8.7, 18.3, 11.9, 6.2, 7.5, 5.8, 8.1, 6.9)
cat("TP53表达值:", tp53_expression, "\n\n")
TP53表达值: 12.5 15.2 8.7 18.3 11.9 6.2 7.5 5.8 8.1 6.9
代码 
# 2. 创建样本名
sample_names <- c(paste0("Tumor_", 1:5), paste0("Normal_", 1:5))
cat("样本名:", sample_names, "\n\n")
样本名: Tumor_1 Tumor_2 Tumor_3 Tumor_4 Tumor_5 Normal_1 Normal_2 Normal_3 Normal_4 Normal_5
代码 
# 3. 统计计算
mean_expr <- mean(tp53_expression)
sd_expr <- sd(tp53_expression)
cat("平均表达量:", round(mean_expr, 2), "\n")
平均表达量: 10.11
代码 
cat("标准差:", round(sd_expr, 2), "\n\n")
标准差: 4.2
代码 
# 4. 高表达样本
high_expr <- tp53_expression[tp53_expression > 10]
cat("高表达样本值:", high_expr, "\n\n")
高表达样本值: 12.5 15.2 18.3 11.9
代码 
# 5. 高表达样本数
n_high <- sum(tp53_expression > 10)
cat("高表达样本数:", n_high, "\n\n")
高表达样本数: 4
代码 
# 6. 最高表达位置
max_pos <- which.max(tp53_expression)
cat("最高表达样本位置:", max_pos, "-", sample_names[max_pos], "\n")
最高表达样本位置: 4 - Tumor_4

练习1.3：患者临床数据框

创建一个包含患者临床信息和基因表达的数据框：

# TODO: 创建如下数据框
# 患者ID: P001, P002, P003, P004, P005
# 年龄: 45, 52, 38, 61, 29
# 性别: M, F, F, M, F
# 肿瘤分期: II, III, I, IV, II
# TP53表达量: 12.5, 15.2, 8.7, 18.3, 11.9

patients <- data.frame(
  # 你的代码：
  
  
  
  
  
  stringsAsFactors = FALSE
)

# 查看数据框结构
str(patients)

# 查看摘要统计
summary(patients)
点击查看参考答案
代码 
# 创建数据框
patients <- data.frame(
  patient_id = c("P001", "P002", "P003", "P004", "P005"),
  age = c(45, 52, 38, 61, 29),
  gender = c("M", "F", "F", "M", "F"),
  tumor_stage = c("II", "III", "I", "IV", "II"),
  tp53_expr = c(12.5, 15.2, 8.7, 18.3, 11.9),
  stringsAsFactors = FALSE
)

# 查看结构
cat("数据框结构:\n")
数据框结构:
代码 
str(patients)
'data.frame':   5 obs. of  5 variables:
 $ patient_id : chr  "P001" "P002" "P003" "P004" ...
 $ age        : num  45 52 38 61 29
 $ gender     : chr  "M" "F" "F" "M" ...
 $ tumor_stage: chr  "II" "III" "I" "IV" ...
 $ tp53_expr  : num  12.5 15.2 8.7 18.3 11.9
代码 
cat("\n摘要统计:\n")

摘要统计:
代码 
summary(patients)
  patient_id             age        gender          tumor_stage       
 Length:5           Min.   :29   Length:5           Length:5          
 Class :character   1st Qu.:38   Class :character   Class :character  
 Mode  :character   Median :45   Mode  :character   Mode  :character  
                    Mean   :45                                        
                    3rd Qu.:52                                        
                    Max.   :61                                        
   tp53_expr    
 Min.   : 8.70  
 1st Qu.:11.90  
 Median :12.50  
 Mean   :13.32  
 3rd Qu.:15.20  
 Max.   :18.30
代码 
cat("\n数据框内容:\n")

数据框内容:
代码 
print(patients)
  patient_id age gender tumor_stage tp53_expr
1       P001  45      M          II      12.5
2       P002  52      F         III      15.2
3       P003  38      F           I       8.7
4       P004  61      M          IV      18.3
5       P005  29      F          II      11.9

练习1.4：临床数据筛选与分析

# TODO: 基于 patients 数据框进行分析

# 1. 提取所有女性患者的信息
female_patients <- ____

# 2. 提取年龄大于50岁的患者ID和TP53表达量
older_patients <- ____

# 3. 计算不同性别患者的平均TP53表达量
# 提示：使用 aggregate() 函数
avg_by_gender <- ____

# 4. 找出TP53表达量最高的患者信息
high_tp53 <- ____

# 5. 添加一列 "风险等级"：
#    TP53表达量 >= 15 为 "高风险"
#    TP53表达量 10-15 为 "中风险"
#    TP53表达量 < 10 为 "低风险"
patients$risk <- ifelse(____, "高风险",
                 ifelse(____, "中风险", "低风险"))

patients

提示： - 条件筛选使用 data[condition, ] - which.max() 可以找到最大值的位置 - aggregate() 用于分组统计

点击查看参考答案
代码 
# 重新创建数据框（确保存在）
patients <- data.frame(
  patient_id = c("P001", "P002", "P003", "P004", "P005"),
  age = c(45, 52, 38, 61, 29),
  gender = c("M", "F", "F", "M", "F"),
  tumor_stage = c("II", "III", "I", "IV", "II"),
  tp53_expr = c(12.5, 15.2, 8.7, 18.3, 11.9),
  stringsAsFactors = FALSE
)

# 1. 提取女性患者
female_patients <- patients[patients$gender == "F", ]
cat("女性患者:\n")
女性患者:
代码 
print(female_patients)
  patient_id age gender tumor_stage tp53_expr
2       P002  52      F         III      15.2
3       P003  38      F           I       8.7
5       P005  29      F          II      11.9
代码 
# 2. 提取年龄大于50岁的患者
older_patients <- patients[patients$age > 50, c("patient_id", "tp53_expr")]
cat("\n年龄大于50岁的患者:\n")

年龄大于50岁的患者:
代码 
print(older_patients)
  patient_id tp53_expr
2       P002      15.2
4       P004      18.3
代码 
# 3. 按性别计算平均TP53表达量
avg_by_gender <- aggregate(tp53_expr ~ gender, data = patients, FUN = mean)
cat("\n不同性别平均TP53表达量:\n")

不同性别平均TP53表达量:
代码 
print(avg_by_gender)
  gender tp53_expr
1      F  11.93333
2      M  15.40000
代码 
# 4. TP53表达量最高的患者
high_tp53 <- patients[which.max(patients$tp53_expr), ]
cat("\nTP53表达量最高的患者:\n")

TP53表达量最高的患者:
代码 
print(high_tp53)
  patient_id age gender tumor_stage tp53_expr
4       P004  61      M          IV      18.3
代码 
# 5. 添加风险等级
patients$risk <- ifelse(patients$tp53_expr >= 15, "高风险",
                 ifelse(patients$tp53_expr >= 10, "中风险", "低风险"))
cat("\n添加风险等级后的数据框:\n")

添加风险等级后的数据框:
代码 
print(patients)
  patient_id age gender tumor_stage tp53_expr   risk
1       P001  45      M          II      12.5 中风险
2       P002  52      F         III      15.2 高风险
3       P003  38      F           I       8.7 低风险
4       P004  61      M          IV      18.3 高风险
5       P005  29      F          II      11.9 中风险

第二部分：向量化操作与函数 (25分钟)

练习2.1：向量化 vs 循环性能对比

比较循环和向量化计算的性能差异：

# TODO: 比较循环和向量化计算10万个基因的标准化值

# 生成10万个基因的表达值
set.seed(123)
gene_expression <- rnorm(100000, mean = 10, sd = 2)

# 方法1：使用循环（慢）
result_loop <- numeric(length(gene_expression))

system.time({
  for (i in 1:length(gene_expression)) {
    result_loop[i] <- (gene_expression[i] - mean(gene_expression)) / sd(gene_expression)
  }
})

# 方法2：使用向量化（快）
system.time({
  result_vectorized <- ____
})

# 验证结果是否相同
all.equal(result_loop, result_vectorized)

# 计算速度提升倍数
# 提示：比较两个system.time的输出
点击查看参考答案
代码 
# 生成数据
set.seed(123)
gene_expression <- rnorm(100000, mean = 10, sd = 2)

# 方法1：循环
cat("=== 循环方法 ===\n")
=== 循环方法 ===
代码 
result_loop <- numeric(length(gene_expression))

time_loop <- system.time({
  for (i in 1:length(gene_expression)) {
    result_loop[i] <- (gene_expression[i] - mean(gene_expression)) / sd(gene_expression)
  }
})
print(time_loop)
  用户   系统   流逝 
54.556  0.187 54.868
代码 
# 方法2：向量化
cat("\n=== 向量化方法 ===\n")

=== 向量化方法 ===
代码 
time_vec <- system.time({
  result_vectorized <- (gene_expression - mean(gene_expression)) / sd(gene_expression)
})
print(time_vec)
 用户  系统  流逝 
0.001 0.000 0.001
代码 
# 验证
cat("\n结果是否相同:", all.equal(result_loop, result_vectorized), "\n")

结果是否相同: TRUE
代码 
# 计算速度提升
speedup <- time_loop["elapsed"] / time_vec["elapsed"]
cat("速度提升:", round(speedup, 1), "倍\n")
速度提升: 54868 倍

练习2.2：生信分析函数

完成差异表达分析相关函数：

# TODO: 完善以下生信分析函数

# 1. 计算Fold Change的函数
calculate_fc <- function(tumor_expr, normal_expr, log2 = TRUE) {
  # 计算均值
  mean_tumor <- ____
  mean_normal <- ____
  
  # 计算Fold Change
  fc <- ____
  
  # 如果log2=TRUE，返回log2 Fold Change
  if (log2) {
    fc <- ____
  }
  
  return(fc)
}

# 2. 基因表达统计函数
analyze_expression <- function(expression_values) {
  # 返回列表，包含：
  # - mean: 平均值
  # - sd: 标准差
  # - min: 最小值
  # - max: 最大值
  # - cv: 变异系数 (sd/mean)
  
  return(list(
    mean = ____,
    sd = ____,
    min = ____,
    max = ____,
    cv = ____
  ))
}

# 测试函数
tumor <- c(12.5, 15.2, 18.3, 14.1, 16.7)
normal <- c(6.2, 7.5, 5.8, 8.1, 6.9)

calculate_fc(tumor, normal)
analyze_expression(c(tumor, normal))
点击查看参考答案
代码 
# 1. Fold Change计算函数
calculate_fc <- function(tumor_expr, normal_expr, log2 = TRUE) {
  mean_tumor <- mean(tumor_expr)
  mean_normal <- mean(normal_expr)
  
  fc <- mean_tumor / mean_normal
  
  if (log2) {
    fc <- log2(fc)
  }
  
  return(fc)
}

# 2. 表达统计函数
analyze_expression <- function(expression_values) {
  return(list(
    mean = mean(expression_values),
    sd = sd(expression_values),
    min = min(expression_values),
    max = max(expression_values),
    cv = sd(expression_values) / mean(expression_values)
  ))
}

# 测试
tumor <- c(12.5, 15.2, 18.3, 14.1, 16.7)
normal <- c(6.2, 7.5, 5.8, 8.1, 6.9)

cat("Fold Change (log2):", round(calculate_fc(tumor, normal), 2), "\n\n")
Fold Change (log2): 1.15
代码 
cat("表达统计:\n")
表达统计:
代码 
stats <- analyze_expression(c(tumor, normal))
print(stats)
$mean
[1] 11.13

$sd
[1] 4.745302

$min
[1] 5.8

$max
[1] 18.3

$cv
[1] 0.4263524

练习2.3：BMI计算与分类

完成BMI计算函数（临床常用）：

# TODO: 完善BMI计算函数
calculate_bmi <- function(weight_kg, height_m) {
  # 1. 计算BMI
  bmi <- ____
  
  # 2. 根据BMI判断分类（中国标准）
  # < 18.5: 偏瘦
  # 18.5 - 23.9: 正常
  # 24 - 27.9: 超重
  # >= 28: 肥胖
  if (____) {
    category <- "偏瘦"
  } else if (____) {
    category <- "正常"
  } else if (____) {
    category <- "超重"
  } else {
    category <- "肥胖"
  }
  
  # 3. 返回列表
  return(list(
    bmi = round(bmi, 2),
    category = category
  ))
}

# 测试函数
calculate_bmi(70, 1.75)
calculate_bmi(55, 1.60)
calculate_bmi(85, 1.70)
点击查看参考答案
代码 
# BMI计算函数
calculate_bmi <- function(weight_kg, height_m) {
  bmi <- weight_kg / (height_m ^ 2)
  
  if (bmi < 18.5) {
    category <- "偏瘦"
  } else if (bmi < 24) {
    category <- "正常"
  } else if (bmi < 28) {
    category <- "超重"
  } else {
    category <- "肥胖"
  }
  
  return(list(
    bmi = round(bmi, 2),
    category = category
  ))
}

# 测试
cat("测试1 (70kg, 1.75m):\n")
测试1 (70kg, 1.75m):
代码 
result1 <- calculate_bmi(70, 1.75)
print(result1)
$bmi
[1] 22.86

$category
[1] "正常"
代码 
cat("\n测试2 (55kg, 1.60m):\n")

测试2 (55kg, 1.60m):
代码 
result2 <- calculate_bmi(55, 1.60)
print(result2)
$bmi
[1] 21.48

$category
[1] "正常"
代码 
cat("\n测试3 (85kg, 1.70m):\n")

测试3 (85kg, 1.70m):
代码 
result3 <- calculate_bmi(85, 1.70)
print(result3)
$bmi
[1] 29.41

$category
[1] "肥胖"

第三部分：综合生信分析 (30分钟)

练习3.1：TP53基因表达差异分析

完成完整的差异表达分析流程：

代码 
# 创建示例数据
samples <- c(paste0("Tumor_", 1:5), paste0("Normal_", 1:5))
groups <- c(rep("Tumor", 5), rep("Normal", 5))
tp53_expr <- c(12.5, 15.2, 8.7, 18.3, 11.9,
               6.2, 7.5, 5.8, 8.1, 6.9)

# 创建数据框
data <- data.frame(
  sample = samples,
  group = groups,
  expression = tp53_expr
)

print(data)
     sample  group expression
1   Tumor_1  Tumor       12.5
2   Tumor_2  Tumor       15.2
3   Tumor_3  Tumor        8.7
4   Tumor_4  Tumor       18.3
5   Tumor_5  Tumor       11.9
6  Normal_1 Normal        6.2
7  Normal_2 Normal        7.5
8  Normal_3 Normal        5.8
9  Normal_4 Normal        8.1
10 Normal_5 Normal        6.9

任务：

# TODO: 完成差异表达分析

# 1. 分组统计：计算Tumor和Normal组的均值和标准差
# 提示：使用 aggregate() 函数
group_stats <- ____
group_stats

# 2. 差异分析：使用t检验比较两组
# 提示：使用 t.test() 函数
tumor <- data$expression[data$group == "Tumor"]
normal <- data$expression[data$group == "Normal"]
t_test_result <- ____

# 3. 计算Fold Change（使用之前定义的函数）
fc <- calculate_fc(tumor, normal)

# 4. 结果汇总
results <- list(
  gene = "TP53",
  mean_tumor = ____,
  mean_normal = ____,
  fold_change = ____,
  p_value = ____,
  significant = ____  # p < 0.05?
)

# 打印结果
cat("=== TP53差异表达分析结果 ===\n")
cat("肿瘤组均值:", round(results$mean_tumor, 2), "\n")
cat("正常组均值:", round(results$mean_normal, 2), "\n")
cat("Log2 Fold Change:", round(results$fold_change, 2), "\n")
cat("P-value:", format(results$p_value, scientific = TRUE, digits = 3), "\n")
cat("显著性:", ifelse(results$significant, "显著", "不显著"), "\n")
点击查看参考答案
代码 
# 重新创建数据
samples <- c(paste0("Tumor_", 1:5), paste0("Normal_", 1:5))
groups <- c(rep("Tumor", 5), rep("Normal", 5))
tp53_expr <- c(12.5, 15.2, 8.7, 18.3, 11.9, 6.2, 7.5, 5.8, 8.1, 6.9)

data <- data.frame(
  sample = samples,
  group = groups,
  expression = tp53_expr
)

# 1. 分组统计
group_stats <- aggregate(expression ~ group, data = data, 
                        FUN = function(x) c(mean = mean(x), sd = sd(x)))
cat("分组统计:\n")
分组统计:
代码 
print(group_stats)
   group expression.mean expression.sd
1 Normal       6.9000000     0.9354143
2  Tumor      13.3200000     3.6182869
代码 
# 2. t检验
tumor <- data$expression[data$group == "Tumor"]
normal <- data$expression[data$group == "Normal"]
t_test_result <- t.test(tumor, normal)

# 3. Fold Change
fc <- calculate_fc(tumor, normal)

# 4. 结果汇总
results <- list(
  gene = "TP53",
  mean_tumor = mean(tumor),
  mean_normal = mean(normal),
  fold_change = fc,
  p_value = t_test_result$p.value,
  significant = t_test_result$p.value < 0.05
)

# 打印
cat("\n=== TP53差异表达分析结果 ===\n")

=== TP53差异表达分析结果 ===
代码 
cat("肿瘤组均值:", round(results$mean_tumor, 2), "\n")
肿瘤组均值: 13.32
代码 
cat("正常组均值:", round(results$mean_normal, 2), "\n")
正常组均值: 6.9
代码 
cat("Log2 Fold Change:", round(results$fold_change, 2), "\n")
Log2 Fold Change: 0.95
代码 
cat("P-value:", format(results$p_value, scientific = TRUE, digits = 3), "\n")
P-value: 1.46e-02
代码 
cat("显著性:", ifelse(results$significant, "显著", "不显著"), "\n")
显著性: 显著

练习3.2：多基因表达矩阵分析

代码 
# 基因表达矩阵（5个基因 × 4个样本）
expr_matrix <- matrix(
  c(12.5, 15.2, 8.7, 18.3,
    8.3, 9.1, 7.9, 8.7,
    15.2, 14.5, 16.1, 13.8,
    20.1, 19.3, 21.5, 18.2,
    6.7, 7.2, 6.1, 7.8),
  nrow = 5, byrow = TRUE
)

dimnames(expr_matrix) <- list(
  c("TP53", "BRCA1", "EGFR", "MYC", "KRAS"),
  c("S1", "S2", "S3", "S4")
)

print(expr_matrix)
        S1   S2   S3   S4
TP53  12.5 15.2  8.7 18.3
BRCA1  8.3  9.1  7.9  8.7
EGFR  15.2 14.5 16.1 13.8
MYC   20.1 19.3 21.5 18.2
KRAS   6.7  7.2  6.1  7.8

任务：

# TODO: 分析表达矩阵

# 1. 计算每个基因在所有样本中的平均表达量
gene_means <- ____

# 2. 计算每个样本中所有基因的平均表达量
sample_means <- ____

# 3. 找出平均表达量最高的基因
most_expressed_gene <- ____

# 4. 计算每个基因的表达标准差
gene_sds <- ____

# 5. 找出表达最稳定的基因（标准差最小）
most_stable_gene <- ____

# 6. 创建基因统计汇总表
gene_summary <- data.frame(
  gene = rownames(expr_matrix),
  mean = ____,
  sd = ____,
  cv = ____  # 变异系数 = sd/mean
)

# 按平均表达量排序
gene_summary <- gene_summary[order(____, decreasing = TRUE), ]
gene_summary
点击查看参考答案
代码 
# 重新创建矩阵
expr_matrix <- matrix(
  c(12.5, 15.2, 8.7, 18.3,
    8.3, 9.1, 7.9, 8.7,
    15.2, 14.5, 16.1, 13.8,
    20.1, 19.3, 21.5, 18.2,
    6.7, 7.2, 6.1, 7.8),
  nrow = 5, byrow = TRUE
)
dimnames(expr_matrix) <- list(
  c("TP53", "BRCA1", "EGFR", "MYC", "KRAS"),
  c("S1", "S2", "S3", "S4")
)

# 1. 每个基因的平均表达
gene_means <- rowMeans(expr_matrix)
cat("每个基因的平均表达:\n")
每个基因的平均表达:
代码 
print(round(gene_means, 2))
 TP53 BRCA1  EGFR   MYC  KRAS 
13.68  8.50 14.90 19.78  6.95
代码 
# 2. 每个样本的平均表达
sample_means <- colMeans(expr_matrix)
cat("\n每个样本的平均表达:\n")

每个样本的平均表达:
代码 
print(round(sample_means, 2))
   S1    S2    S3    S4 
12.56 13.06 12.06 13.36
代码 
# 3. 表达量最高的基因
most_expressed_gene <- names(which.max(gene_means))
cat("\n表达量最高的基因:", most_expressed_gene, "\n")

表达量最高的基因: MYC
代码 
# 4. 每个基因的标准差
gene_sds <- apply(expr_matrix, 1, sd)

# 5. 表达最稳定的基因
most_stable_gene <- names(which.min(gene_sds))
cat("表达最稳定的基因:", most_stable_gene, "\n")
表达最稳定的基因: BRCA1
代码 
# 6. 基因统计汇总
gene_summary <- data.frame(
  gene = rownames(expr_matrix),
  mean = gene_means,
  sd = gene_sds,
  cv = gene_sds / gene_means
)

gene_summary <- gene_summary[order(gene_summary$mean, decreasing = TRUE), ]
cat("\n基因表达统计汇总:\n")

基因表达统计汇总:
代码 
print(gene_summary)
       gene   mean        sd         cv
MYC     MYC 19.775 1.3889444 0.07023739
EGFR   EGFR 14.900 0.9831921 0.06598605
TP53   TP53 13.675 4.0762524 0.29808061
BRCA1 BRCA1  8.500 0.5163978 0.06075268
KRAS   KRAS  6.950 0.7234178 0.10408889

练习3.3：质控过滤

# TODO: 基于 gene_summary 进行质控过滤

# 过滤条件：
# 1. 平均表达量 >= 10
# 2. 变异系数 <= 0.15

# 1. 应用过滤条件
filtered_genes <- gene_summary[____, ]

# 2. 统计信息
cat("原始基因数:", nrow(gene_summary), "\n")
cat("过滤后基因数:", nrow(filtered_genes), "\n")
cat("过滤比例:", round((1 - nrow(filtered_genes)/nrow(gene_summary)) * 100, 1), "%\n")

# 3. 查看通过的基因
cat("\n通过质控的基因:\n")
print(filtered_genes)
点击查看参考答案
代码 
# 重新创建 gene_summary
gene_summary <- data.frame(
  gene = c("TP53", "BRCA1", "EGFR", "MYC", "KRAS"),
  mean = c(13.68, 8.50, 14.90, 19.78, 6.95),
  sd = c(4.03, 0.69, 1.12, 1.37, 0.72),
  cv = c(0.2945, 0.0812, 0.0752, 0.0693, 0.1036)
)

# 质控过滤
filtered_genes <- gene_summary[gene_summary$mean >= 10 & gene_summary$cv <= 0.15, ]

cat("通过质控的基因:\n")
通过质控的基因:
代码 
print(filtered_genes)
  gene  mean   sd     cv
3 EGFR 14.90 1.12 0.0752
4  MYC 19.78 1.37 0.0693
代码 
cat("\n原始基因数:", nrow(gene_summary), "\n")

原始基因数: 5
代码 
cat("过滤后基因数:", nrow(filtered_genes), "\n")
过滤后基因数: 2
代码 
cat("过滤比例:", round((1 - nrow(filtered_genes)/nrow(gene_summary)) * 100, 1), "%\n")
过滤比例: 60 %

第四部分：R包管理 (10分钟)

练习4.1：安装Bioconductor包

# TODO: 安装Bioconductor包

# 1. 安装BiocManager（如果尚未安装）
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# 2. 安装DESeq2包（RNA-seq差异表达分析）
# BiocManager::install("DESeq2")

# 3. 查看已安装的包
installed_packages <- ____
head(installed_packages[, c("Package", "Version")], 10)

# 4. 加载包并查看帮助
library(ggplot2)
help(____)  # 查看ggplot函数的帮助
点击查看参考答案
代码 
# 1. 安装BiocManager
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# 2. 安装DESeq2（取消注释运行）
# BiocManager::install("DESeq2")

# 3. 查看已安装包
installed_packages <- installed.packages()
head(installed_packages[, c("Package", "Version")], 10)

# 4. 加载和查看帮助
library(ggplot2)
help(ggplot)

实验提交

提交内容

  1. R脚本文件 (lab1.R)
    • 包含所有练习的代码
    • 添加必要注释
    • 确保代码可运行
  2. 实验报告 (PDF或HTML，使用R Markdown编写)
    • 包含关键代码和结果截图
    • 回答思考题
    • 总结学习心得

思考题

  1. 向量化 vs 循环
    • 为什么R中向量化操作通常比循环更快？
    • 在什么情况下必须使用循环？
  2. 数据类型选择
    • data.frame 和 matrix 的主要区别是什么？
    • 什么时候应该使用 list？
  3. 函数设计
    • 如何设计一个”好”的函数？
    • 参数命名有什么最佳实践？
  4. 生信应用
    • 在差异表达分析中，Fold Change和P值分别代表什么？
    • 为什么要进行质控过滤？

延伸阅读

实验结束，祝学习愉快！ 🎉