Skip to content
Snippets Groups Projects
Commit 27a85648 authored by Todor Kondić's avatar Todor Kondić
Browse files

Refactor MS1 chromatogram extraction.

parent 03780562
No related branches found
No related tags found
No related merge requests found
...@@ -335,45 +335,83 @@ extr_data <-function(m) { ...@@ -335,45 +335,83 @@ extr_data <-function(m) {
dpath = m$run$paths$data dpath = m$run$paths$data
## Big ms1 chromatogram table.
cg1 = new_ms1_cgm_table()
cg2 = new_ms2_cgm_table()
for (fn in fine[,unique(file)]) {
## Open all files.
fns = fine[,unique(file)]
lms = lapply(fns,function(fn) read_data_file(file.path(dpath,fn)))
names(lms) = fns
## Load all feature data.
lfdata = lapply(lms,get_fdata)
names(lfdata) = fns
## Extract MS1 chromatograms using "fine" tolerance.
cgram_ms1 = fine[,extr_cgrams_ms1(lms[[file]],.SD,lfdata[[file]]),
by="file",
.SDcols=c("iso_fine_min",
"iso_fine_max",
"rt_min",
"rt_max",
"precid")]
setkey(cgram_ms1,file,precid)
## Extract MS2 chromatograms.
cgram_ms2 = data.table(precid=integer(0),
an=integer(0),
ce=numeric(0),
rt=numeric(0),
intensity=numeric(0))
for (fn in names(lfdata)) {
browser()
x = lfdata[[file]]$ms2[cgram_ms1,.(an,ce,precid=i.precid),on="prec_idx==idx"]
}
1+1
## cg2 = new_ms2_cgm_table()
for (fn in fns) {
## Read data. ## Read data.
ms = read_data_file(file=file.path(dpath,fn)) ms = lms[[fn]]
## ## Get input entries for a particular file.
## ftab = fine[.(fn),on="file"]
## Get input entries for a particular file. ## Get feature data.
ftab = fine[.(fn),on="file"] fdata = lfdata[[fn]]
## Extract input data that is needed. ## Get tables with mz fine ranges and retention times (if
isotab = ftab[,.(iso_fine_min,iso_fine_max,rt_min,rt_max),by="precid"] ## any), used to filter the MS1 chromatogram.
isotab1 = isotab[is.na(rt_min),.(precid,iso_fine_min,iso_fine_max)] cginputs = get_inputs_4_cgram(ftab)
isotab2 = isotab[!is.na(rt_min),.(precid,iso_fine_min,iso_fine_max,rt_min,rt_max)]
## Extract MS1 chromatograms. ## Extract MS1 chromatograms.
fncg1 = new_ms1_cgm_table() cg1 = extr_ms1_cgm(ms=ms,
fncg1 = extr_ms1_cgm(ms=ms, isotab=cginputs$rtno,
isotab=isotab1, qrt=F,
qrt=F, cg1)
fncg1)
fncg1 = extr_ms1_cgm(ms=ms, cg1 = extr_ms1_cgm(ms=ms,
isotab=isotab2, isotab=cginputs$rtyes,
qrt=T, qrt=T,
fncg1) cg1)
cg1 = cg1[fncg1,.(precid,
rt,
intensity=i.intensity,
scan=i.scan)]
## ## Annotate fdata with precids.
## fdata$ms1[fncg1,precid:=i.precid,on="scan"]
## Extract MS2 chromatograms. ## Extract MS2 chromatograms.
fnfd = as.data.table(fData(ms),keep.rownames="scan",key="scan")
fncg2 = extr_ms2_cgm(fnfd,fncg1)
1+1
} }
......
...@@ -103,41 +103,3 @@ make_db_precursors <- function(m) { ...@@ -103,41 +103,3 @@ make_db_precursors <- function(m) {
m$db$precursors = masses m$db$precursors = masses
m m
} }
new_ms1_cgm_table <- function() {
## Creates a chromatogram.
data.table(precid=integer(0),
rt=numeric(0),
intensity=numeric(0),
scan=character(0),
key = c("precid","rt"))
}
new_ms2_cgm_cat <- function() {
## A kind of catalogue of ms2.
data.table(precid = integer(0),
an = integer(0))
}
new_ms2_cgm_table <- function() {
## Varios rt-related MS2 data.
data.table(precid = integer(0),
ce = numeric(0),
an = integer(0),
rt = numeric(0))
}
new_ms2_cgm_table <- function() {
## Creates a chromatogram.
data.table(precid = integer(0),
an = integer(0),
rt = numeric(0),
intensity = numeric(0),
scan = character(0),
key = c("precid","an"))
}
...@@ -611,6 +611,7 @@ gen_ms2_spec_blk <- function(spectra) { ...@@ -611,6 +611,7 @@ gen_ms2_spec_blk <- function(spectra) {
## NEW FUNCTIONS. ## NEW FUNCTIONS.
create_fine_table <- function(m) { create_fine_table <- function(m) {
## Select fine mz-ranges and split them into those with rt entries ## Select fine mz-ranges and split them into those with rt entries
## and those without. ## and those without.
...@@ -641,48 +642,56 @@ read_data_file <- function(file) { ...@@ -641,48 +642,56 @@ read_data_file <- function(file) {
MSnbase::readMSData(file=file,msLevel=c(1,2),mode="onDisk") MSnbase::readMSData(file=file,msLevel=c(1,2),mode="onDisk")
} }
extr_ms1_cgm <- function(ms,isotab,qrt,res) {
## Extract chromatograms.
## Get mz ranges in matrix format.
mzrng = as.matrix(isotab[,.(iso_fine_min,iso_fine_max)])
x = if (!qrt) {
## Call without rt argument.
MSnbase::chromatogram(ms,mz = mzrng)
} else {
## Call with rt argument (in seconds).
rtrng = as.matrix(isotab[,.(rt_min*60,rt_max*60)])
MSnbase::chromatogram(ms,mz = mzrng, rt = rtrng)
}
## If there were any input masses actually, extr_cgrams_ms1 <- function(ms,tab,fdata) {
if (dim(mzrng)[[1L]] > 0L) { ## Some helpers.
## fill the data table. new_restab <- function(intab,cgm) {
for (i in 1L:nrow(mzrng)) { base = intab[,.(precid=precid,cgmidx=.I)]
rt = rtime(x[i,1]) cgm = base[,{
precid = isotab[(i),precid] rt = rtime(cgm[cgmidx,1])
chunk = data.table(precid=precid, inte = intensity(cgm[cgmidx,1])
rt=rt, .(precid=precid,
intensity=intensity(x[i,1]), rt = rt,
scan = names(rt), intensity = inte,
key = c("precid","rt")) scan = names(rt))},
by="cgmidx"]
res = res[chunk,.(precid, setkey(cgm,scan)
rt, cgm[fdata$ms1,idx:=i.idx,on="scan"]
intensity=i.intensity, cgm
scan=i.scan)]
}
} }
trt = tab[!is.na(rt_min)]
tnort = tab[is.na(rt_min)]
resrt = if (nrow(trt)>0L) {
## Call with rt argument (in seconds).
mzrng = as.matrix(trt[,.(iso_fine_min,iso_fine_max)])
rtrng = as.matrix(trt[,.(rt_min*60,rt_max*60)])
new_restab(trt,MSnbase::chromatogram(ms,mz = mzrng, rt = rtrng))
} else data.table()
resnort = if (nrow(tnort)>0L) {
mzrng = as.matrix(tnort[,.(iso_fine_min,iso_fine_max)])
new_restab(tnort,MSnbase::chromatogram(ms,mz = mzrng))
} else data.table()
rbind(resnort,resrt,fill=T)
}
get_fdata <- function(ms) {
fdata = as.data.table(fData(ms),keep.rownames="scan")
setkey(fdata,scan)
res = list()
res$ms1 = fdata[msLevel==1L,.(scan,
idx=spIdx)]
res$ms2 = fdata[msLevel==2L,.(scan,
idx=spIdx,
an=acquisitionNum,
ce=collisionEnergy,
prec_mz=precursorMZ,
prec_idx=precursorScanNum)]
res res
}
extr_ms2_cgm <- function(feat_table,one_file_cg1) {
## Get scanIdx.
browser()
precs = feat_table[one_file_cg1,.(scan,scanIdx),on="scan"]
} }
...@@ -10,9 +10,7 @@ test_that("Extraction returns what is needed.",{ ...@@ -10,9 +10,7 @@ test_that("Extraction returns what is needed.",{
phase=c("setup","comptab","db")) phase=c("setup","comptab","db"))
cat = m$db$cat cat = m$db$cat
print(m$db$precursors[,.(mz,file)])
m = run(envopts=eo,m=m,phase="extract") m = run(envopts=eo,m=m,phase="extract")
print(m$db$extr$cgm$ms1)
expect_true(1==1) expect_true(1==1)
}) })
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment