蛋白组数据预处理可参考该文章：
(2 封私信) 生信分析系列干货 | 高分文章蛋白质组学数据预处理方法调研报告 - 知乎
DIA数据---中值归一化
setwd("E:\\OMV课题研究\\DIP-Fe-蛋白组\\R自行分析") # ① 加载所有包 library(readxl) library(dplyr) library(MSnbase) library(impute) # ② 读入数据 expr <- readxl::read_xlsx("protein_matrix.xlsx") expr <- as.data.frame(expr) rownames(expr) <- expr[[1]] expr[[1]] <- NULL expr <- expr %>% mutate_all(as.numeric) # ③ 将0替换为NA expr[expr == 0] <- NA cat("原始数据行数:", nrow(expr), "\n") cat("原始数据列数:", ncol(expr), "\n") # ④ log2转换（先） log2_raw <- log2(expr) # ⑤ 中位数归一化（再）——log2空间用减法 medians <- apply(log2_raw, 2, median, na.rm = TRUE) log2_norm <- sweep(log2_raw, 2, medians, "-") cat("\n各列归一化后中位数（应接近0）：\n") print(round(apply(log2_norm, 2, median, na.rm = TRUE), 4)) # ⑥ 检查归一化后分布 hist(as.matrix(log2_norm), main = "log2归一化后分布", xlab = "log2(normalized intensity)", breaks = 50, col = "lightblue") # ⑦ 缺失值过滤（后）——至少一组三个样本全部存在才保留 # 前3列为 control，后3列为 treat ctrl_idx <- 1:3 treat_idx <- 4:6 ctrl_complete <- apply(log2_norm[, ctrl_idx], 1, function(x) all(!is.na(x))) treat_complete <- apply(log2_norm[, treat_idx], 1, function(x) all(!is.na(x))) keep <- ctrl_complete | treat_complete cat("\n过滤前行数:", nrow(log2_norm), "\n") cat("过滤后行数:", sum(keep), "\n") cat(" 其中 control 组完整:", sum(ctrl_complete), "\n") cat(" 其中 treat 组完整:", sum(treat_complete), "\n") cat(" 两组均完整 :", sum(ctrl_complete & treat_complete), "\n") log2_filtered <- log2_norm[keep, ] # ⑧ 各列缺失率 col_missing_rate <- apply(log2_filtered, 2, function(x) sum(is.na(x)) / length(x)) cat("\n各列缺失率(%)：\n") print(round(col_missing_rate * 100, 2)) # ⑨ KNN填补 set.seed(42) knn_result <- impute.knn(as.matrix(log2_filtered), k = 10, rowmax = 0.5, colmax = 0.8) imputed_data <- knn_result$data # ⑩ 验证 cat("\n填充后是否还有NA：", any(is.na(imputed_data)), "\n") cat("数值范围：", range(imputed_data), "\n") # ⑪ 保存（log2空间） write.csv(imputed_data, "expr_knn.csv") cat("\n已保存至 expr_knn.csv\n")bacteria （DDA数据--FOT归一化）
setwd("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria") ###########计算每一组的完整率################# # 读取数据 annotation_reference3 <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\refenrence3_gokegg.csv") bacteria_row<- read.csv("filt_bacteria.csv",row.names = 1,check.names = F) bacteria_row <- bacteria_row[,-1] # 定义计算完整率的函数 calculate_completeness <- function(row, start_col, end_col) { # 提取当前组的列 group_data <- row[start_col:end_col] # 计算非零值的数量 non_zero_count <- sum(as.numeric(group_data) > 0, na.rm = TRUE) # 计算完整率 completeness <- non_zero_count / (end_col - start_col + 1) return(completeness) } # 计算每一组的完整率 bacteria_row$group1_completeness <- apply(bacteria_row, 1, function(row) { calculate_completeness(row, 1, 10) }) bacteria_row$group2_completeness <- apply(bacteria_row, 1, function(row) { calculate_completeness(row, 11, 20) }) bacteria_row$group3_completeness <- apply(bacteria_row, 1, function(row) { calculate_completeness(row, 21, 25) }) bacteria_row$group4_completeness <- apply(bacteria_row, 1, function(row) { calculate_completeness(row, 26, 30) }) ########筛选特异性蛋白################## # 条件：其中一组大于0.5，其他三组小于等于0.2 filtered_rows <- bacteria_row[ ( (bacteria_row$group1_completeness > 0.5 & bacteria_row$group2_completeness <= 0.2 & bacteria_row$group3_completeness <= 0.2 & bacteria_row$group4_completeness <= 0.2) | (bacteria_row$group2_completeness > 0.5 & bacteria_row$group1_completeness <= 0.2 & bacteria_row$group3_completeness <= 0.2 & bacteria_row$group4_completeness <= 0.2) | (bacteria_row$group3_completeness > 0.5 & bacteria_row$group1_completeness <= 0.2 & bacteria_row$group2_completeness <= 0.2 & bacteria_row$group4_completeness <= 0.2) | (bacteria_row$group4_completeness > 0.5 & bacteria_row$group1_completeness <= 0.2 & bacteria_row$group2_completeness <= 0.2 & bacteria_row$group3_completeness <= 0.2) ), ] filtered_rows$Entry <- rownames(filtered_rows) unique_annotation_bacteria <- left_join(filtered_rows,annotation_reference3,by="Entry") write.csv(unique_annotation_bacteria,"unique_annotation_bacteria.csv") ######筛选后续分析蛋白###### bacteria_filt <- bacteria_row[ ( (bacteria_row$group1_completeness > 0.5 | (bacteria_row$group2_completeness > 0.5| (bacteria_row$group3_completeness > 0.5| (bacteria_row$group4_completeness > 0.5))))), ] bacteria_filt <- bacteria_filt[,-c(31:34)] write.csv(bacteria_filt,"bacteria_filt.csv") ##################归一化########################## # 加载必要的库 library(dplyr) # 计算每个样本的总量，忽略NA值 data <- bacteria_filt sample_sums <- colSums(data, na.rm = TRUE) # 选择Sample1作为参考样本 reference_sample <- sample_sums[1] # 计算归一化系数 normalization_factors <- sample_sums / reference_sample # 归一化数据，忽略NA值 normalized_data <- sweep(data, 2, normalization_factors, "/") # 计算归一化后的平均值，忽略NA值 averages <- rowMeans(normalized_data, na.rm = TRUE) # 均一化数据，忽略NA值 scaled_data <- normalized_data / averages * 100 ###################log2转换#################### log_data <- log2(scaled_data+1) boxplot(as.data.frame(log_data),main="Original") #据此是否归一化###################log2转换#################### ####################缺失值填充################### min_value <- min(log_data, na.rm = TRUE) # 将所有缺失值填充为最小值 data_filled <- log_data data_filled[is.na(data_filled)] <- min_value write.csv(data_filled,"bacteria_filled_scaled.csv") #####计算每一组的平均值######## data1 <- data_filled data1$'CR-hv'<-apply(data1[,1:10],1,mean) data1$'CS-hv'<-apply(data1[,11:20],1,mean) data1$'CR-c'<-apply(data1[,21:25],1,mean) data1$'CS-c'<-apply(data1[,26:30],1,mean) #提取平均值列到新的数据集 groupmean<-data1[,31:34] ##############limma差异分析################# setwd("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\差异分析") bacteria_card <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\VF_CARD\\bacteria_card.csv") bacteria_filled <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\bacteria_filled_scaled.csv",row.names = 1,check.names = F) bacteriaexpr_card <- bacteria_filled[bacteria_card$id,] rownames(bacteriaexpr_card) <- bacteria_card$name library(limma) #构建分组矩阵 data_filled <- bacteriaexpr_card[,c(1:20,21:30)] group <- c(rep("treat",20),rep("con",10)) group <- factor(group,levels = c("con","treat")) design <- model.matrix(~0+group) colnames(design) <- levels(group) design #构建比较矩阵 contrast.matrix <- makeContrasts(treat - con,levels=design) #线性拟合模型构建 fit <- lmFit(data_filled,design)#线性模型拟合 fit <- contrasts.fit(fit, contrast.matrix) fit <- eBayes(fit)#贝叶斯检验 #最终得到差异分析结果 allDiff=topTable(fit,number=Inf) data <- allDiff select.log2FC <- abs(data$logFC)>0.585 select.Pvalue <- data$adj.P.Val<0.05 select.vec <- select.log2FC & select.Pvalue table(select.vec) data$change <- "Normal" data$change[data$logFC>=0.585 & data$adj.P.Val<0.05]="Up" data$change[data$logFC<=-0.585 & data$adj.P.Val<0.05]="Down" write.csv(data,"./毒力非毒力之间的耐药比较/hv vs c_allDiff.CSV") ###################火山图############################ library(ggplot2) library(ggrepel) data$logP <- -log10(data$adj.P.Val) # 创建火山图 ggplot(data, aes(x = logFC, y = logP, color = change)) + xlim(-7, 12) + ylim(0, 10) + # 设置x轴和y轴的范围 geom_point(alpha=0.4,size = 3.5) + # 画散点图，调整点的大小 theme_bw() + # 使用经典主题 scale_color_manual(values = c("blue4", "grey", "red3")) + # 点的颜色 geom_hline(yintercept = -log10(0.05), linetype = 4, size = 0.5,lwd=0.5) + # 添加水平线 geom_vline(xintercept = c(-0.585, 0.585), linetype = 4, size = 0.5,lwd=0.5) + # 添加垂直线 theme( title = element_text(size = 14), # 设置标题字体大小 text = element_text(size = 15) # 设置文本字体大小 ) + labs(x = "log2(fold change)", y = "-log10(adj.P-value)") # 设置坐标轴标签 ###################################热图############################## library(pheatmap) # 加载pheatmap这个R包 #读取热图数据文件 #df = read.delim("https://www.bioladder.cn/shiny/zyp/demoData/heatmap/data.heatmap.txt", #文件名称 注意文件路径，格式 header = T, # 是否有标题 sep = "\t", # 分隔符是Tab键 row.names = 1, # 指定第一列是行名 fill=T) # 是否自动填充，一般选择是 #读取分组数据文件 #dfSample = read.delim("https://www.bioladder.cn/shiny/zyp/demoData/heatmap/sample.class.txt",header = T,row.names = 1,fill = T,sep = "\t") #dfGene = read.delim("https://www.bioladder.cn/shiny/zyp/demoData/heatmap/gene.class.txt",header = T,row.names = 1,fill = T,sep = "\t") #绘图 #pheatmap(df, annotation_row=dfGene, # （可选）指定行分组文件 annotation_col=dfSample, # （可选）指定列分组文件 show_colnames = TRUE, # 是否显示列名 show_rownames=TRUE, # 是否显示行名 fontsize=10, # 字体大小 color = colorRampPalette(c('#0000ff','#ffffff','#ff0000'))(50), # 指定热图的颜色 annotation_legend=TRUE, # 是否显示图例 border_color=NA, # 边框颜色 NA表示没有 scale="row", # 指定归一化的方式。"row"按行归一化，"column"按列归一化，"none"不处理 cluster_rows = TRUE, # 是否对行聚类 cluster_cols = TRUE # 是否对列聚类 ) bacteria_card <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\VF_CARD\\bacteria_card.csv") bacteria_vf <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\VF_CARD\\bacteria_VF.csv") bacteria_filled <- read.csv("C:\\Users\\HUAWEI\\Desktop\\更正2蛋白组\\bacteria\\bacteria_filled_scaled.csv",row.names = 1,check.names = F) bacteriaexpr_vf <- bacteria_filled[bacteria_vf$id,] rownames(bacteriaexpr_vf) <- bacteria_vf$p_names dfProtein <- data.frame(bacteria_vf[,4]) rownames(dfProtein) <- bacteria_vf$p_names colnames(dfProtein) <- "Mechanism" dfSample <- data.frame(Group =c(rep("CR-hv",10),rep("CS-hv",10),rep("CR-c",5),rep("CS-c",5))) rownames(dfSample) <-colnames(bacteriaexpr_vf) df <- bacteriaexpr_vf pdf("bacteriaheatmap_vf.pdf",height = 6,width = 10) pheatmap(df, annotation_col=dfSample, annotation_row=dfProtein,# （可选）指定列分组文件 fontsize=10, # 字体大小 show_colnames = F, # 是否显示列名 show_rownames=T, color = colorRampPalette(c("navy","white","firebrick3"))(50), # 指定热图的颜色 annotation_legend=T , # 是否显示图例 border_color=NA, # 边框颜色 NA表示没有 scale="row", # 指定归一化的方式。"row"按行归一化，"column"按列归一化，"none"不处理 cluster_rows = T, # 是否对行聚类 cluster_cols = F,# 是否对列聚类 annotation_names_col = TRUE, # 显示列注释名称 annotation_names_row = F, annotation_colors = list( Group = c("CR-hv" = "#FFD47F","CS-hv" = "#F7C1CF","CR-c" = "#7B92C7","CS-c"="#ADD9EE"))) dev.off()
