百度360必应搜狗淘宝本站头条
HTML5 标签标签HTML5 标签 (已废弃)HTML 简介键盘快捷键
当前位置：网站首页 > 技术文章 > 正文

分析代谢组数据，大名鼎鼎的xcms来了

zhezhongyun 2025-01-27 01:14 25 浏览

xcms是基于R语言设计的程序包(R package)，可以用分析代谢组数据等。下面我们就来介绍一下使用方法。


操作步骤

示例数据是fatty acid amide hydrolase (FAAH) 基因敲除鼠的脊柱LC-MS，六个基因敲除个体，六个野生型，使用的是centroid mode ，正离子模式，200-600 m/z ，2500-4500 seconds。


source("https://bioconductor.org/biocLite.R")

yum install netcdf

yum install netcdf-devel.x86_64

biocLite("xcms")

biocLite("ncdf4")

biocLite("faahKO")

biocLite("MSnbase")

library(xcms)

library(faahKO)

library(RColorBrewer)

library(pander)



#读取数据


cdfs <- dir(system.file("cdf", package = "faahKO"), full.names = TRUE,recursive = TRUE)

> cdfs

[1] "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko15.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko16.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko18.CDF"

[4] "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko19.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko21.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/KO/ko22.CDF"

[7] "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt15.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt16.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt18.CDF"

[10] "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt19.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt21.CDF" "C:/Program Files/R/R-3.4.2/library/faahKO/cdf/WT/wt22.CDF"


#构建样式矩阵


pd <- data.frame(sample_name = sub(basename(cdfs), pattern = ".CDF",replacement = "", fixed = TRUE),sample_group = c(rep("KO", 6), rep("WT", 6)),stringsAsFactors = FALSE)

> pd

sample_name sample_group

1 ko15 KO

2 ko16 KO

3 ko18 KO

4 ko19 KO

5 ko21 KO

6 ko22 KO

7 wt15 WT

8 wt16 WT

9 wt18 WT

10 wt19 WT

11 wt21 WT

12 wt22 WT


#读取数据


raw_data <- readMSData(files = cdfs, pdata = new("NAnnotatedDataFrame", pd),mode = "onDisk")



#查看保留时间


> head(rtime(raw_data))

F01.S0001 F01.S0002 F01.S0003 F01.S0004 F01.S0005 F01.S0006

2501.378 2502.943 2504.508 2506.073 2507.638 2509.203


#查看质荷比


> head(mz(raw_data))

$F01.S0001

[1] 200.1 201.0 201.9 202.9 203.8 204.2 205.1 206.0 207.0 208.0 209.1 210.0 211.0 212.0 213.0 214.0 215.1 216.1 217.1 218.0 219.0 220.0 220.9 222.0 223.1 224.1 225.0 226.0 227.1 228.0


#查看强度


> head(intensity(raw_data))

$F01.S0001

[1] 1716 1723 2814 1961 667 676 1765 747 2044 757 1810 926 3381 1442 1688 1223 1465 1624 2446 1309 2167 900 5471 873 2285 1355 2610 1797 6494 2314


#按文件分割数据


mzs <- mz(raw_data)

mzs_by_file <- split(mzs, f = fromFile(raw_data))

length(mzs_by_file)


#总体评价图


bpis <- chromatogram(raw_data, aggregationFun = "max")

group_colors <- brewer.pal(3, "Set1")[1:2]

names(group_colors) <- c("KO", "WT")

plot(bpis, col = group_colors[raw_data$sample_group])



#查看某个个体


bpi_1 <- bpis[1, 1]

plot(bpi_1, col = group_colors[raw_data$sample_group])



#查看样品离子流


tc <- split(tic(raw_data), f = fromFile(raw_data))

boxplot(tc, col = group_colors[raw_data$sample_group],ylab = "intensity", main = "Total ion current")


#定义保留时间与质荷比，提取特定的峰


rtr <- c(2700, 2900)

mzr <- c(334.9, 335.1)

chr_raw <- chromatogram(raw_data, mz = mzr, rt = rtr)

plot(chr_raw, col = group_colors[chr_raw$sample_group])


#提取质谱数据


msd_raw <- extractMsData(raw_data, mz = mzr, rt = rtr)

plotMsData(msd_raw[[1]])



#定义峰宽与噪音，自动找出所有的峰


cwp <- CentWaveParam(peakwidth = c(30, 80), noise = 1000)

xdata <- findChromPeaks(raw_data, param = cwp)

head(chromPeaks(xdata))


#统计检测到的峰


summary_fun <- function(z) {

c(peak_count = nrow(z), rt = quantile(z[, "rtmax"] - z[, "rtmin"]))

}

T <- lapply(split.data.frame(chromPeaks(xdata),

f = chromPeaks(xdata)[, "sample"]),

FUN = summary_fun)

T <- do.call(rbind, T)

rownames(T) <- basename(fileNames(xdata))

pandoc.table(T,

caption = paste0("Summary statistics on identified chromatographic",

" peaks. Shown are number of identified peaks per",

" sample and widths/duration of chromatographic ",

"peaks."))



#画某个样品的峰图


plotChromPeaks(xdata, file = 3)


#对所有样品画热图


plotChromPeakImage(xdata)



#标注某个峰的差异


plot(chr_raw, col = group_colors[chr_raw$sample_group], lwd = 2)

highlightChromPeaks(xdata, border = group_colors[chr_raw$sample_group],lty = 3, rt = rtr, mz= mzr)


#提取某个峰所有样品的数据


pander(chromPeaks(xdata, mz = mzr, rt = rtr),caption = paste("Identified chromatographic peaks in a selected ","m/z and retention time range."))



#画峰强的箱线图


ints <- split(log2(chromPeaks(xdata)[, "into"]),f = chromPeaks(xdata)[, "sample"])

boxplot(ints, varwidth = TRUE, col = group_colors[xdata$sample_group],ylab =expression(log[2]~intensity), main = "Peak intensities")

grid(nx = NA, ny = NULL)



# 设定binSize


xdata <- adjustRtime(xdata, param = ObiwarpParam(binSize = 0.6))



#查看校正过的保留时间


head(adjustedRtime(xdata))



#查看校正前的保留时间


> head(rtime(xdata, adjusted = FALSE))



#画校正保留时间后的峰图及校正后与校正前的差异


bpis_adj <- chromatogram(xdata, aggregationFun = "max")

par(mfrow = c(2, 1), mar = c(4.5, 4.2, 1, 0.5))

plot(bpis_adj, col = group_colors[bpis_adj$sample_group])


#查看数据是否校正过时间


hasAdjustedRtime(xdata)

[1] TRUE


#恢复到没校正的状态


xdata <- dropAdjustedRtime(xdata)

hasAdjustedRtime(xdata)

[1] FALSE


#根据样品组设定参数


pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,minFraction = 0.8)

#根据组来提取数据

xdata <- groupChromPeaks(xdata, param = pdp)

#根据组来校正时间

pgp <- PeakGroupsParam(minFraction = 0.85)

xdata <- adjustRtime(xdata, param = pgp)

#对校正前和校正后的结果作图

plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group],peakGroupsCol = "grey", peakGroupsPch = 1)


#对校正前和校正后的某个峰作图


par(mfrow = c(2, 1))

plot(chr_raw, col = group_colors[chr_raw$sample_group])

chr_adj <- chromatogram(xdata, rt = rtr, mz = mzr)

plot(chr_adj, col = group_colors[chr_raw$sample_group])



#选择一种质荷比，提取峰图


mzr <- c(305.05, 305.15)

chr_mzr <- chromatogram(xdata, mz = mzr, rt = c(2500, 4000))

par(mfrow = c(3, 1), mar = c(1, 4, 1, 0.5))

cols <- group_colors[chr_mzr$sample_group]

plot(chr_mzr, col = cols, xaxt = "n", xlab = "")

highlightChromPeaks(xdata, mz = mzr, col = cols, type = "point", pch = 16)



#画一级质谱和二级质谱的峰图


mzr <- c(305.05, 305.15)

chr_mzr_ms1 <- chromatogram(filterMsLevel(xdata, 1), mz = mzr, rt = c(2500, 4000))

plot(chr_mzr_ms1)

chr_mzr_ms2 <- chromatogram(filterMsLevel(xdata, 2), mz = mzr, rt = c(2500, 4000))

plot(chr_mzr_ms2)


#定义峰提取参数


pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,minFraction = 0.4, bw = 30)

par(mar = c(4, 4, 1, 0.5))

plotChromPeakDensity(xdata, mz = mzr, col = cols, param = pdp,pch = 16, xlim = c(2500, 4000))


#用不同的峰提取参数，bw定义了距离多少的两个峰合并为一个峰


pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,minFraction = 0.4, bw = 20)

plotChromPeakDensity(xdata, mz = mzr, col = cols, param = pdp,pch = 16, xlim = c(2500, 4000))


#提取数据，minFraction是占所有样本百分之多少以上的峰视为正确的数据


pdp <- PeakDensityParam(sampleGroups = xdata$sample_group,minFraction = 0.4, bw = 20)

xdata <- groupChromPeaks(xdata, param = pdp)

featureDefinitions(xdata)



#对结果分组

head(featureValues(xdata, value = "into"))

#利用原始数据对NA进行回填

xdata <- fillChromPeaks(xdata)

head(featureValues(xdata))

#查看回填前的NA

apply(featureValues(xdata, filled = FALSE), MARGIN = 2,FUN = function(z) sum(is.na(z)))

#查看回填后的NA

apply(featureValues(xdata), MARGIN = 2,FUN = function(z) sum(is.na(z)))

#PCA


ft_ints <- log2(featureValues(xdata, value = "into"))

pc <- prcomp(t(na.omit(ft_ints)), center = TRUE)

cols <- group_colors[xdata$sample_group]

pcSummary <- summary(pc)

plot(pc$x[, 1], pc$x[,2], pch = 21, main = "", xlab = paste0("PC1: ", format(pcSummary$importance[2, 1] * 100,digits = 3), " % variance"),ylab = paste0("PC2: ", format(pcSummary$importance[2, 2] * 100,digits = 3), " % variance"),col = "darkgrey", bg = cols, cex = 2)

grid()

text(pc$x[, 1], pc$x[,2], labels = xdata$sample_name, col = "darkgrey",pos = 3, cex = 2)


#查看数据处理历史，经过了Peak detection、Peak grouping、Retention time correction、Peak grouping、Missing peak filling

processHistory(xdata)


#提取某一步的数据

ph <- processHistory(xdata, type = "Retention time correction")

ph

#提取数据的参数

processParam(ph[[1]])

#提取某个文件的数据

subs <- filterFile(xdata, file = c(2, 4))

#提取数据并留取保留时间

subs <- filterFile(xdata, keepAdjustedRtime = TRUE)

#按保留时间提取数据

subs <- filterRt(xdata, rt = c(3000, 3500))

range(rtime(subs))

#提取某个文件的所有数据

one_file <- filterFile(xdata, file = 3)

one_file_2 <- xdata[fromFile(xdata) == 3]

#查看峰文件

head(chromPeaks(xdata))


#导出数据

result<-cbind(as.data.frame(featureDefinitions(xdata)),featureValues(xdata, value = "into"))

write.table(result,file="xcms_result.txt",sep="\t",quote=F)


#PCA


values <- groupval(xdata)

data <- t(values)

pca.result <- pca(data)

library(pcaMethods)

pca.result <- pca(data)

loadings <- pca.result@loadings

scores <- pca.result@scores

plotPcs(pca.result, type = "scores",col=as.integer(sampclass(xdata)) + 1)


#MDS


library(MASS)

for (r in 1:ncol(data))

data[,r] <- data[,r] / max(data[,r])

data.dist <- dist(data)

mds <- isoMDS(data.dist)

plot(mds$points, type = "n")

text(mds$points, labels = rownames(data),col=as.integer(sampclass(xdata))+ 1 )


好了，以上就是使用xcms分析代谢组数据的操作。如果文章对你有所帮助，请转发给你身边需要的人噢！

关注我们Get更多科研小工具


相关推荐

JPA实体类注解，看这篇就全会了

基本注解@Entity标注于实体类声明语句之前，指出该Java类为实体类，将映射到指定的数据库表。name(可选)：实体名称。缺省为实体类的非限定名称。该名称用于引用查询中的实体。不与@Tab...

Dify教程02 - Dify+Deepseek零代码赋能，普通人也能开发AI应用

开始今天的教程之前，先解决昨天遇到的一个问题，docker安装Dify的时候有个报错，进入Dify面板的时候会出现“InternalServerError”的提示，log日志报错：S3_USE_A...

用离散标记重塑人体姿态:VQ-VAE实现关键点组合关系编码

在人体姿态估计领域，传统方法通常将关键点作为基本处理单元，这些关键点在人体骨架结构上代表关节位置（如肘部、膝盖和头部）的空间坐标。现有模型对这些关键点的预测主要采用两种范式：直接通过坐标回归或间接通过...

B 客户端流RPC (clientstream Client Stream)

客户端编写一系列消息并将其发送到服务器，同样使用提供的流。一旦客户端写完消息，它就等待服务器读取消息并返回响应gRPC再次保证了单个RPC调用中的消息排序在客户端流RPC模式中，客户端会发送多个请...

我的模型我做主02——训练自己的大模型:简易入门指南

模型训练往往需要较高的配置，为了满足友友们的好奇心，这里我们不要内存，不要gpu，用最简单的方式，让大家感受一下什么是模型训练。基于你的硬件配置，我们可以设计一个完全在CPU上运行的简易模型训练方案。...

开源项目MessageNest打造个性化消息推送平台多种通知方式

今天介绍一个开源项目，MessageNest-可以打造个性化消息推送平台，整合邮件、钉钉、企业微信等多种通知方式。定制你的消息，让通知方式更灵活多样。开源地址：https://github.c...

使用投机规则API加快页面加载速度

当今的网络用户要求快速导航，从一个页面移动到另一个页面时应尽量减少延迟。投机规则应用程序接口（SpeculationRulesAPI）的出现改变了网络应用程序接口（WebAPI）领域的游戏规则。...

JSONP安全攻防技术

关于JSONPJSONP全称是JSONwithPadding，是基于JSON格式的为解决跨域请求资源而产生的解决方案。它的基本原理是利用HTML的元素标签，远程调用JSON文件来实现数据传递。如果...

大数据Doris(六):编译 Doris遇到的问题

编译Doris遇到的问题一、js_generator.cc:(.text+0xfc3c):undefinedreferenceto`well_known_types_js’查找Doris...

网页内嵌PDF获取的办法

最近女王大人为了通过某认证考试，交了2000RMB，官方居然没有给线下教材资料，直接给的是在线教材，教材是PDF的但是是内嵌在网页内，可惜却没有给具体的PDF地址，无法下载，看到女王大人一点点的截图保...

印度女孩被邻居家客人性骚扰，父亲上门警告，反被围殴致死

微信的规则进行了调整希望大家看完故事多点“在看”，喜欢的话也点个分享和赞这样事儿君的推送才能继续出现在你的订阅列表里才能继续跟大家分享每个开怀大笑或拍案惊奇的好故事啦~话说只要稍微关注新闻的人，应该...

下周重要财经数据日程一览 (1229-0103)

下周焦点全球制造业PMI美国消费者信心指数美国首申失业救济人数值得注意的是，下周一希腊还将举行第三轮总统选举需要谷歌日历同步及部分智能手机（安卓，iPhone）同步日历功能的朋友请点击此链接，数据公布...

PyTorch 深度学习实战(38):注意力机制全面解析

在上一篇文章中，我们探讨了分布式训练实战。本文将深入解析注意力机制的完整发展历程，从最初的Seq2Seq模型到革命性的Transformer架构。我们将使用PyTorch实现2个关键阶段的注意力机制变...

聊聊Spring AI的EmbeddingModel

序本文主要研究一下SpringAI的EmbeddingModelEmbeddingModelspring-ai-core/src/main/java/org/springframework/ai/e...

前端分享-少年了解过iframe么

iframe就像是HTML的「内嵌画布」，允许在页面中加载独立网页，如同在画布上叠加另一幅动态画卷。核心特性包括：独立上下文：每个iframe都拥有独立的DOM/CSS/JS环境（类似浏...

一周热门
最近发表
标签列表