求助：Roary输出.Rtab基因数据清洗及韦恩图绘制方案

一、用Awk预处理数据（清理+统计） #

通过Awk过滤掉所有菌株均无该基因（全0）的行，同时统计各基因分类的数量：

awk '
NR==1 {next}                # 跳过表头行
$2+$3+$4 == 0 {next}        # 删除三个菌株均为0的行
{
    pattern = $2""$3""$4    # 拼接三个菌株的存在状态为标识
    count[pattern]++
}
END {
    # 按需求输出各分类的基因数
    print "核心基因(111):", count["111"]
    print "Accessory a+b(110):", count["110"]
    print "Accessory b+c(011):", count["011"]
    print "Accessory a+c(101):", count["101"]
    print "unique a(100):", count["100"]
    print "unique b(010):", count["010"]
    print "unique c(001):", count["001"]
}' input.Rtab > gene_counts.txt

执行后会生成gene_counts.txt文件，包含各分类的基因数量统计结果。

二、用RStudio绘制韦恩图 #

提供两种实现方式，可根据习惯选择：

方式1：基于Awk统计结果绘制 #

# 安装依赖包（首次运行需执行）
install.packages("VennDiagram")

# 加载包
library(VennDiagram)

# 读取统计数据
counts <- read.table("gene_counts.txt", sep=":", stringsAsFactors=FALSE)
counts$V2 <- as.integer(trimws(counts$V2))  # 清理空格并转为整数

# 提取各分类数值
core <- counts[counts$V1 == "核心基因(111)", "V2"]
a_b <- counts[counts$V1 == "Accessory a+b(110)", "V2"]
b_c <- counts[counts$V1 == "Accessory b+c(011)", "V2"]
a_c <- counts[counts$V1 == "Accessory a+c(101)", "V2"]
unique_a <- counts[counts$V1 == "unique a(100)", "V2"]
unique_b <- counts[counts$V1 == "unique b(010)", "V2"]
unique_c <- counts[counts$V1 == "unique c(001)", "V2"]

# 生成韦恩图
venn.diagram(
    x = list(
        StrainA = rep(1, unique_a + a_b + a_c + core),
        StrainB = rep(1, unique_b + a_b + b_c + core),
        StrainC = rep(1, unique_c + a_c + b_c + core)
    ),
    filename = "gene_venn.png",
    fill = c("#FF6B6B", "#4ECDC4", "#45B7D1"),
    alpha = 0.5,
    cat.cex = 1.1,
    main = "基因存在/缺失韦恩图",
    main.cex = 1.3
)

运行后会在工作目录生成gene_venn.png韦恩图文件。

方式2：直接处理原始.Rtab文件 #

无需Awk预处理，直接在R中完成数据清理、统计和绘图：

# 安装依赖包（首次运行需执行）
install.packages("VennDiagram")

# 加载包
library(VennDiagram)

# 读取原始.Rtab文件
gene_data <- read.table("input.Rtab", header = TRUE, stringsAsFactors = FALSE)

# 过滤全0行
filtered_data <- gene_data[rowSums(gene_data[, 2:4]) != 0, ]

# 绘制韦恩图
venn.diagram(
    x = list(
        StrainA = which(filtered_data$StrainA == 1),
        StrainB = which(filtered_data$StrainB == 1),
        StrainC = which(filtered_data$StrainC == 1)
    ),
    filename = "gene_venn_direct.png",
    fill = c("#E74C3C", "#3498DB", "#2ECC71"),
    alpha = 0.6,
    cat.pos = c(-30, 30, 180),
    cat.dist = 0.12,
    main = "Roary 基因存在/缺失分析"
)

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