你给的建议很好,但是还是画不了图。我不用VIP以及他的filter,请检查我的代码:> # 诊断 expr 的类型
> class(expr)
[1] "data.frame"
> str(expr)
'data.frame': 161 obs. of 89 variables:
$ P16_3 : num 67961 73167 320130 4810 6134 ...
$ P36_3 : num 52925 44644 212557 22673 2472 ...
$ P43_3 : num 101281 63640 257677 19905 1526 ...
$ P46_3 : num 45285 55774 318102 19930 5156 ...
$ P47_3 : num 36459 59168 144029 16062 5380 ...
$ P48_3 : num 50039 43018 217704 17692 3336 ...
$ P50_3 : num 142051 135219 312985 36347 13803 ...
$ P51_3 : num 102816 70851 381088 31105 7889 ...
$ P52_3 : num 181595 86421 363186 17640 4404 ...
$ P54_3 : num 60192 66310 269933 17393 3591 ...
$ P55_3 : num 112513 76074 377750 10409 5777 ...
$ P56_3 : num 63428 91943 272754 19098 7613 ...
$ P57_3 : num 53328 53201 435370 14830 2016 ...
$ P59_3 : num 92725 62218 379374 12953 3689 ...
$ P60_3 : num 112171 80718 264273 32179 9370 ...
$ P61_3 : num 60565 100499 272527 23977 4612 ...
$ P62_3 : num 48670 52608 98222 5578 3504 ...
$ P63_3 : num 114774 77114 271784 18584 5601 ...
$ P65_3 : num 49578 61393 194017 15295 3373 ...
$ P66_3 : num 41415 41540 151165 9391 15594 ...
$ P67_3 : num 96594 56126 202153 6491 2760 ...
$ P68_3 : num 81789 58478 229933 10618 6264 ...
$ P69_3 : num 105061 68379 360112 8881 4169 ...
$ P70_3 : num 38773 47576 236634 10249 6530 ...
$ P71_3 : num 51274 43792 241896 9435 7655 ...
$ P72_3 : num 59814 42736 197908 9421 19136 ...
$ P73_3 : num 33669 69854 198013 14996 4121 ...
$ P74_3 : num 86124 38867 198205 8631 8794 ...
$ P76_3 : num 97585 62923 279572 18384 3032 ...
$ P77_3 : num 66949 54318 208424 10820 2663 ...
$ P78_3 : num 77367 53656 276004 5962 4517 ...
$ P79_3 : num 72756 63482 233224 9619 5108 ...
$ P80_3 : num 39551 46655 302934 16321 4596 ...
$ P81_3 : num 31690 39816 208057 9416 12160 ...
$ P84_3 : num 119611 75287 284189 9701 15640 ...
$ P85_3 : num 95084 91458 279023 21939 12734 ...
$ P86_3 : num 112292 82161 266722 8601 12648 ...
$ P88_3 : num 46441 66239 272629 8037 7225 ...
$ P90_3 : num 168933 85360 606717 23284 5520 ...
$ P16_1 : num 21.5 40681.4 46612.5 122701.8 7664.7 ...
$ P36_1 : num 22 139753 58252 338929 14448 ...
$ P41_1 : num 20.7 69095.4 46032.4 280569.2 40188.4 ...
$ P43_1 : num 22.7 79791.3 55343.3 255376.4 25606.3 ...
$ P45_1 : num 21.8 152495 70270.4 273617.2 11787.2 ...
$ P46_1 : num 21.8 136249.9 49970.2 405493.1 10059.3 ...
$ P47_1 : num 22.6 62073.4 64774.4 139282.5 8248.2 ...
$ P48_1 : num 22.9 372832.5 135796.7 661319.7 29087.9 ...
$ P49_1 : num 21.4 45771.9 52365.5 216571.3 7617.1 ...
$ P49_1b: num 21.4 46999.4 51212.7 214301.1 8107.5 ...
$ P50_1 : num 28.4 197789.5 77415.7 304874.1 18558.8 ...
$ P51_1 : num 26.6 119592.3 100841.1 408581.7 22670.8 ...
$ P52_1 : num 16.9 148078.5 76749.4 302118.3 16679 ...
$ P54_1 : num 20.6 117074.4 68253.2 461336.5 8106.3 ...
$ P55_1 : num 19.4 109283.2 99106.1 340218.4 8625.9 ...
$ P56_1 : num 22.8 47497.5 53221.2 164683.8 6228.9 ...
$ P57_1 : num 21.6 155447.8 74952.2 454655.6 10733.2 ...
$ P59_1 : num 24.6 89654.5 58800.7 308751.7 7407 ...
$ P60_1 : num 22.9 185536.4 96177.3 284125.3 15737.8 ...
$ P61_1 : num 19.4 103549.7 98127 420091.7 8466.2 ...
$ P62_1 : num 22.5 74573.2 63210 205234.6 10337.9 ...
$ P63_1 : num 22.6 117431.6 88632.6 494701.1 17176.8 ...
$ P65_1 : num 22 95796 62277 231937 5411 ...
$ P66_1 : num 20.1 64070.2 46722.3 255246.2 9304.1 ...
$ P67_1 : num 27.6 104951.4 64030.3 246236 7704.3 ...
$ P68_1 : num 26.1 106696.8 75855.5 265469.7 11721.8 ...
$ P69_1 : num 25.9 167477.1 89639.3 492901.9 8993.1 ...
$ P70_1 : num 25.6 88506.5 58846.3 222547.4 11396.8 ...
$ P71_1 : num 19.8 58958 70824.2 213633.9 7770.6 ...
$ P72_1 : num 22.6 63663.2 50911.5 139844.2 6106.1 ...
$ P73_1 : num 21.1 59101.7 62963.3 306584.4 7313.9 ...
$ P74_1 : num 19.8 71750.2 39777.2 211361.8 5414.4 ...
$ P75_1 : num 25.8 95451.2 78107 230287.5 10401.4 ...
$ P76_1 : num 21.1 103836.4 83110 220113.7 12030.3 ...
$ P77_1 : num 23.5 229180.4 69815.6 394408.4 7417.1 ...
$ P78_1 : num 29.2 157328.1 71213.7 272158.1 5756.9 ...
$ P79_1 : num 20.1 65004.3 71737.8 243271.2 6913.3 ...
$ P80_1 : num 37.7 145372.1 77984.5 212674 16661.1 ...
$ P81_1 : num 20.6 186422.6 101018.8 246182.6 17902.6 ...
$ P81_1b: num 20.6 116079.1 63878.5 191847 9013.7 ...
$ P84_1 : num 23.9 73586.5 62464.7 186149.7 5505.6 ...
$ P84_1b: num 23.9 40100.4 29762.6 156280.5 4522.2 ...
$ P85_1 : num 26.2 186562.7 70271.4 373234.8 18496.7 ...
$ P85_1b: num 26.2 112463.6 46783.4 304528.3 12101.9 ...
$ P86_1 : num 19.6 111811.8 62070.5 257121.8 16918.5 ...
$ P86_1b: num 19.6 85547.3 63011.4 217995.3 12866.6 ...
$ P88_1 : num 15.6 113491.2 78512.4 217897.4 7291 ...
$ P88_1b: num 15.6 93704.7 64928.8 206878.2 4876.7 ...
$ P90_1 : num 20 215425 92211 341472 28393 ...
$ P90_1b: num 20 174678 77047 343126 18825 ...
>
> # 诊断 Type 的类型
> class(Type)
[1] "factor"
> str(Type)
Factor w/ 2 levels "3m","Pre": 1 1 1 1 1 1 1 1 1 1 ...
- attr(*, "names")= chr [1:89] "P16_3" "P36_3" "P43_3" "P46_3" ...
>
> # 因子转字符
> Type_char <- as.character(Type)
>
> # 字符转数值(无法转的会是 NA)
> Type_num <- as.numeric(Type_char)
警告信息:
强制改变过程中产生了NA
>
> # 结果查看
> Type_num
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[31] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[61] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
> x <- colnames(expr) # 你真正要检查的对象
>
> class(x)
[1] "character"
> str(x)
chr [1:89] "P16_3" "P36_3" "P43_3" "P46_3" "P47_3" "P48_3" "P50_3" "P51_3" "P52_3" "P54_3" ...
>
> # 字符向量转数值(全是字符,所以会变成 NA)
> x_numeric <- as.numeric(x)
警告信息:
强制改变过程中产生了NA
>
> # 如果你强行当成因子处理并提取数值码:
> x_factor <- factor(x)
> x_numeric2 <- as.numeric(as.character(x_factor)) # 会全部 NA
警告信息:
强制改变过程中产生了NA
>
> # 若有 NA 结果,跳过该代谢物(避免进入判断)
> if(any(is.na(c(conMeans, treatMeans, conMed, treatMed, diffMed, logFC, pvalue)))){
+ next
+ }
错误: 找不到对象'conMed'
> outTab <- data.frame()
> for(i in rownames(expr)){
+ metData <- data.frame(
+ expression = as.numeric(expr[i, ]),
+ Type = Type
+ )
+
+ # t 检验自动忽略 NA
+ tTest <- t.test(expression ~ Type, data = metData)
+ pvalue <- tTest$p.value
+
+ # --- 关键修改:全部加上 na.rm = TRUE ---
+
+ conMeans <- mean(expr[i, Type == "Pre"], na.rm = TRUE)
+ treatMeans <- mean(expr[i, Type == "3m"], na.rm = TRUE)
+ logFC <- log2(treatMeans + 1) - log2(conMeans + 1)
+
+ conMed <- median(expr[i, Type == "Pre"], na.rm = TRUE)
+ treatMed <- median(expr[i, Type == "3m"], na.rm = TRUE)
+ diffMed <- treatMed - conMed
+
+ # 若有 NA 结果,跳过该代谢物(避免进入判断)
+ if(any(is.na(c(conMeans, treatMeans, conMed, treatMed, diffMed, logFC, pvalue)))){
+ next
+ }
+
+ # 一致性方向筛选
+ if(((logFC > 0) & (diffMed > 0)) | ((logFC < 0) & (diffMed < 0))){
+ outTab <- rbind(outTab,
+ cbind(
+ metabolite = i,
+ conMean = conMeans,
+ treatMean = treatMeans,
+ logFC = logFC,
+ pvalue = pvalue
+ ))
+ }
+ }
错误于median.default(expr[i, Type == "Pre"], na.rm = TRUE):
需要数值数据
此外: 警告信息:
1: In mean.default(expr[i, Type == "Pre"], na.rm = TRUE) :
参数不是数值也不是逻辑值:返回NA
2: In mean.default(expr[i, Type == "3m"], na.rm = TRUE) :
参数不是数值也不是逻辑值:返回NA
> row.names(outTab) <- outTab$metabolite
> outTab <- outTab[order(as.numeric(outTab$pvalue)),]
> # 输出差异结果 ----------------------------------------------------------
> write.table(outTab, file="met.all.txt", sep="\t", row.names=FALSE, quote=FALSE)
> # 显著差异
> outDiff <- outTab[
+ abs(as.numeric(outTab$logFC)) > logFCfilter &
+ as.numeric(outTab$pvalue) < pvalueFilter,
+ ]
> write.table(outDiff, file="met.diff.txt", sep="\t", row.names=FALSE, quote=FALSE)
> # 热图 -----------------------------------------------------------------
> metaboliteNum = 20
> outDiff <- outDiff[order(outDiff$logFC),]
错误于order(outDiff$logFC): 参数1不是向量
> diffNames <- outDiff$metabolite
> diffLen <- length(diffNames)
> if(diffLen > 2 * metaboliteNum){
+ hmNames <- diffNames[c(1:metaboliteNum, (diffLen-metaboliteNum+1):diffLen)]
+ } else {
+ hmNames <- diffNames
+ }
> hmExpr <- log2(expr[hmNames,] + 1)
> Type_df <- data.frame(Type = Type)
> rownames(Type_df) <- colnames(expr)
> pdf("met.heatmap.pdf", width=9, height=6)
> pheatmap(hmExpr,
+ annotation_col = Type_df,
+ color = colorRampPalette(c("blue", "white", "red"))(50),
+ cluster_cols = FALSE,
+ show_colnames = FALSE,
+ scale = "row")
错误于seq.default(-m, m, length.out = n + 1):
'from' must be a finite number
此外: 警告信息:
1: In min(x, na.rm = T) : min里所有的参数都不存在; 返回Inf
2: In max(x, na.rm = T) : max里所有的参数都不存在;返回-Inf
> dev.off()
null device
1
> # 火山图 + VIP ----------------------------------------------------------
> # 读取 VIP 文件
> vipTab <- read.table(vipFile, sep="\t", header=TRUE, check.names=FALSE)
错误: 找不到对象'vipFile'
> colnames(vipTab) <- c("metabolite", "VIP")
错误: 找不到对象'vipTab'
> # 合并 VIP
> outTab$VIP <- vipTab$VIP[match(outTab$metabolite, vipTab$metabolite)]
错误: 找不到对象'vipTab'
> # 火山图显著性定义
> outTab$pvalue <- as.numeric(outTab$pvalue)
> outTab$logFC <- as.numeric(outTab$logFC)
> Significant = ifelse(
+ pvalueFilter & abs(outTab$logFC) > logFCfilter,
+ ifelse(outTab$logFC > 0, "Up", "Down"),
+ "Not"
+ )
> # 火山图绘制 --------------------------------------------------------------
> p <- ggplot(outTab, aes(logFC, -log10(pvalue))) +
+ geom_point(aes(color = Significant, size = VIP), alpha = 0.8) +
+ scale_color_manual(values = c("Down"="#6BAED6", "Not"="grey", "Up"="#FD7446")) +
+ theme_bw() +
+ labs(title = "Volcano Plot: 3m vs Pre")
> pdf("met.vol.pdf", width=6, height=5)
> print(p)
Error in `geom_point()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 1st layer.
Caused by error:
! 找不到对象'VIP'
Run `rlang::last_trace()` to see where the error occurred.
> # 火山图绘制 --------------------------------------------------------------
> p <- ggplot(outTab, aes(logFC, -log10(pvalue))) +
+ geom_point(aes(color = Significant, size = 2.4), alpha = 0.8) +
+ scale_color_manual(values = c("Down"="#6BAED6", "Not"="grey", "Up"="#FD7446")) +
+ theme_bw() +
+ labs(title = "Volcano Plot: 3m vs Pre")
> pdf("met.vol.pdf", width=6, height=5)
> print(p)
警告信息:
No shared levels found between `names(values)` of the manual scale and the data's colour
values.