蛋白组数据预处理可参考该文章(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,byEntry) 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_data1) boxplot(as.data.frame(log_data),mainOriginal) #据此是否归一化###################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(~0group) colnames(design) - levels(group) design #构建比较矩阵 contrast.matrix - makeContrasts(treat - con,levelsdesign) #线性拟合模型构建 fit - lmFit(data_filled,design)#线性模型拟合 fit - contrasts.fit(fit, contrast.matrix) fit - eBayes(fit)#贝叶斯检验 #最终得到差异分析结果 allDifftopTable(fit,numberInf) data - allDiff select.log2FC - abs(data$logFC)0.585 select.Pvalue - data$adj.P.Val0.05 select.vec - select.log2FC select.Pvalue table(select.vec) data$change - Normal data$change[data$logFC0.585 data$adj.P.Val0.05]Up data$change[data$logFC-0.585 data$adj.P.Val0.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(alpha0.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,lwd0.5) # 添加水平线 geom_vline(xintercept c(-0.585, 0.585), linetype 4, size 0.5,lwd0.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, # 指定第一列是行名 fillT) # 是否自动填充一般选择是 #读取分组数据文件 #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_rowdfGene, # 可选指定行分组文件 annotation_coldfSample, # 可选指定列分组文件 show_colnames TRUE, # 是否显示列名 show_rownamesTRUE, # 是否显示行名 fontsize10, # 字体大小 color colorRampPalette(c(#0000ff,#ffffff,#ff0000))(50), # 指定热图的颜色 annotation_legendTRUE, # 是否显示图例 border_colorNA, # 边框颜色 NA表示没有 scalerow, # 指定归一化的方式。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_coldfSample, annotation_rowdfProtein,# 可选指定列分组文件 fontsize10, # 字体大小 show_colnames F, # 是否显示列名 show_rownamesT, color colorRampPalette(c(navy,white,firebrick3))(50), # 指定热图的颜色 annotation_legendT , # 是否显示图例 border_colorNA, # 边框颜色 NA表示没有 scalerow, # 指定归一化的方式。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()