使用mapply批量绘制SNP-基因表达图时遇行数不匹配错误求助

问题描述 #

我可通过以下代码成功绘制单个SNP（如rs80002225）与特定基因（如Gene1）的小提琴图：

rs_total = read.csv("~/rs_130.csv", header = T, row.names = 1, check.names = FALSE )
rs_info = data.frame(t(rs_total))

xp_total = read.csv("~/xp_130.csv", header = T, row.names = 1,check.names = FALSE )
xp_info = data.frame(t(xp_total))

samp_id = colnames(rs_total, do.NULL = TRUE, prefix = "col")

rs_df = read.csv("~/Gene1_SNPs.csv", header = F,check.names = FALSE)

###############################################################################
  rs80002225 = data.frame(rs_info$rs80002225)
  Gene1 = data.frame(xp_info$Gene1)
  rs80002225 = cbind(rs80002225, samp_id)
  Gene1 = cbind(Gene1, samp_id)
  rs80002225_Gene1 = rs80002225 %>% inner_join( Gene1, by=c('samp_id'='samp_id'))
  colnames(rs80002225_Gene1)[1] = "rs80002225"
  colnames(rs80002225_Gene1)[3] = "Gene1"
  rs80002225_Gene1 %>%
    mutate(rs80002225 = factor(rs80002225,
                               levels=c("0",
                                        "1",
                                        "2")))
  zero_color = "#bf75d1"
  one_color =  "#778fe6"
  two_color = "#86b08e"
  
  rs_count = rs80002225_Gene1 %>%
    count(rs80002225)
  zero_n = rs_count %>%
    filter(rs80002225 == "0") %>%
    pull(n)
  one_n = rs_count %>%
    filter(rs80002225 == "1") %>%
    pull(n)
  two_n = rs_count %>%
    filter(rs80002225 == "2") %>%
    pull(n)
  
  
  rs80002225_Gene1 %>%
  ggplot(aes(x=rs80002225, y=Gene1, fill=factor(rs80002225))) +
    geom_violin(trim = FALSE)+
    # geom_boxplot(show.legend=FALSE, outlier.shape=NA, alpha=0.25, width=0.6,
    #              coef=0)+
    stat_summary(fun.data = median_hilow, fun.args=0.50, show.legend=FALSE,
                 geom="crossbar", alpha=0.25, width=0.6) +
    geom_jitter(show.legend=FALSE, width=0.25, shape=21, color="black") +
    scale_fill_manual(values = c(zero_color, one_color, two_color)) +
    labs(x = "rs80002225", y = "Gene1_Expression", fill = "")+
  theme_classic() +
  theme(axis.text.x = element_markdown())

ggsave(paste0("~/Buetow lab/eqtl_analysis/MatrixEQTL/violin/", "Gene1", "_", "rs80002225", ".png"), device="png")

但使用自定义函数RS_plots结合mapply批量绘制多个SNP对应同一基因的图时，出现错误：

Error in data.frame(..., check.names = FALSE) :  arguments imply differing number of rows: 0, 130

批量代码如下：

RS_plots = function(rs_id, gene_id){
  rs_id = data.frame(rs_info$rs_id)
  gene_id = data.frame(xp_info$gene_id)
  rs_id = cbind(rs_id, samp_id)
  gene_id = cbind(gene_id, samp_id)
  rs_id_gene_id = rs_id %>% inner_join( gene_id, by=c('samp_id'='samp_id'))
  colnames(rs_id_gene_id)[1] = "rs_id"
  colnames(rs_id_gene_id)[3] = "gene_id"
  rs_id_gene_id %>%
    mutate(rs_id = factor(rs_id,
                               levels=c("0",
                                        "1",
                                        "2")))
  zero_color = "#bf75d1"
  one_color =  "#778fe6"
  two_color = "#86b08e"
  
  rs_count = rs_id_gene_id %>%
    count(rs_id)
  zero_n = rs_count %>%
    filter(rs_id == "0") %>%
    pull(n)
  one_n = rs_count %>%
    filter(rs_id == "1") %>%
    pull(n)
  two_n = rs_count %>%
    filter(rs_id == "2") %>%
    pull(n)
  
  
  rs_id_gene_id %>%
  ggplot(aes(x=rs_id, y=gene_id, fill=factor(rs_id))) +
    geom_violin(trim = FALSE)+
    # geom_boxplot(show.legend=FALSE, outlier.shape=NA, alpha=0.25, width=0.6,
    #              coef=0)+
    stat_summary(fun.data = median_hilow, fun.args=0.50, show.legend=FALSE,
                 geom="crossbar", alpha=0.25, width=0.6) +
    geom_jitter(show.legend=FALSE, width=0.25, shape=21, color="black") +
    scale_fill_manual(values = c(zero_color, one_color, two_color)) +
    labs(x = "rs_id", y = "gene_id_Expression", fill = "")+
  theme_classic() +
  theme(axis.text.x = element_markdown())

ggsave(paste0("~/violin/", gene_id, "_", rs_id, ".png"), device="png")
}
  
hopefully = mapply(FUN = RS_plots, rs_id = rs_df$V1, gene_id = "gene_id")

问题原因及修复方案 #

核心错误点 #

1. 动态列引用失败：rs_info$rs_id是字面量引用，无法根据传入的rs_id参数动态提取对应SNP列，需改用rs_info[[rs_id]]；同理xp_info$gene_id要改为xp_info[[gene_id]]。
2. 参数名被覆盖：函数内用rs_id = data.frame(...)会把传入的SNP ID字符串覆盖为数据框，导致后续保存文件名出错，需给数据框用独立变量名（如rs_data、gene_data）。
3. 因子转换未生效：mutate操作未赋值回原数据框，导致rs_id仍为字符型而非因子，需添加赋值语句。
4. 调用参数错误：mapply中传入的gene_id = "gene_id"是无效字符串，需替换为实际基因名（如"Gene1"）。

修复后的完整代码 #

RS_plots = function(rs_id, gene_id){
  # 动态提取SNP和基因表达数据，避免覆盖参数名
  rs_data = data.frame(rs_info[[rs_id]])
  gene_data = data.frame(xp_info[[gene_id]])
  
  rs_data = cbind(rs_data, samp_id)
  gene_data = cbind(gene_data, samp_id)
  
  rs_gene_df = rs_data %>% inner_join(gene_data, by=c('samp_id'='samp_id'))
  colnames(rs_gene_df)[1] = rs_id  # 用实际SNP名命名列
  colnames(rs_gene_df)[3] = gene_id  # 用实际基因名命名列
  
  # 赋值回原数据框，确保因子转换生效
  rs_gene_df = rs_gene_df %>%
    mutate(!!rs_id := factor(!!sym(rs_id), levels=c("0", "1", "2")))
  
  zero_color = "#bf75d1"
  one_color =  "#778fe6"
  two_color = "#86b08e"
  
  rs_count = rs_gene_df %>%
    count(!!sym(rs_id))
  zero_n = rs_count %>% filter(!!sym(rs_id) == "0") %>% pull(n)
  one_n = rs_count %>% filter(!!sym(rs_id) == "1") %>% pull(n)
  two_n = rs_count %>% filter(!!sym(rs_id) == "2") %>% pull(n)
  
  # 动态映射变量名绘图
  p = rs_gene_df %>%
    ggplot(aes(x=!!sym(rs_id), y=!!sym(gene_id), fill=factor(!!sym(rs_id)))) +
    geom_violin(trim = FALSE)+
    stat_summary(fun.data = median_hilow, fun.args=0.50, show.legend=FALSE,
                 geom="crossbar", alpha=0.25, width=0.6) +
    geom_jitter(show.legend=FALSE, width=0.25, shape=21, color="black") +
    scale_fill_manual(values = c(zero_color, one_color, two_color)) +
    labs(x = rs_id, y = paste0(gene_id, "_Expression"), fill = "")+
    theme_classic() +
    theme(axis.text.x = element_markdown())
  
  # 指定绘图对象保存，避免空文件
  ggsave(paste0("~/violin/", gene_id, "_", rs_id, ".png"), plot = p, device="png")
}
  
# 传入实际基因名调用函数
hopefully = mapply(FUN = RS_plots, rs_id = rs_df$V1, gene_id = "Gene1")

额外注意事项 #

• 使用!!sym(var_name)实现dplyr和ggplot中的动态变量引用，避免硬编码。
• 提前检查rs_df$V1中的SNP ID是否全部存在于rs_info的列名中，若有不存在的SNP会生成空数据框，仍可能触发行数不匹配错误。

内容的提问来源于stack exchange，提问作者Daisy