Ricca authored on 26/10/2021 13:08:30
| ... | ... |
@@ -176,10 +176,11 @@ betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 176 | 176 |
#' @export |
| 177 | 177 |
#' @noRd |
| 178 | 178 |
setMethod("show","betweenAlignment",
|
| 179 |
- function(object){
|
|
| 179 |
+ function(object) {
|
|
| 180 | 180 |
cat("An object of class \"", class(object), "\"\n", sep = "")
|
| 181 | 181 |
cat(length(object@mergedPeaksDataset@peaksrt), "groups:", |
| 182 |
- sapply(object@mergedPeaksDataset@peaksrt, length), "merged peaks\n" |
|
| 182 |
+ sapply(object@mergedPeaksDataset@peaksrt, length), |
|
| 183 |
+ "merged peaks\n" |
|
| 183 | 184 |
) |
| 184 | 185 |
} |
| 185 | 186 |
) |
Ricca authored on 10/04/2019 13:24:37
| ... | ... |
@@ -1,7 +1,53 @@ |
| 1 |
- |
|
| 1 |
+#' Data Structure for "between" alignment of many GCMS samples |
|
| 2 |
+#' |
|
| 3 |
+#' This function creates a "between" alignment (i.e. comparing merged peaks) |
|
| 4 |
+#' |
|
| 5 |
+#' \code{betweenAlignment} objects gives the data structure which stores the
|
|
| 6 |
+#' result of an alignment across several "pseudo" datasets. These pseudo |
|
| 7 |
+#' datasets are constructed by merging the "within" alignments. |
|
| 8 |
+#' |
|
| 9 |
+#' @aliases betweenAlignment betweenAlignment-class betweenAlignment-show show, |
|
| 10 |
+#' betweenAlignment-method |
|
| 11 |
+#' @param pD a \code{peaksDataset} object
|
|
| 12 |
+#' @param cAList \code{list} of \code{clusterAlignment} objects, one for each
|
|
| 13 |
+#' experimental group |
|
| 14 |
+#' @param pAList \code{list} of \code{progressiveAlignment} objects, one for
|
|
| 15 |
+#' each experimental group |
|
| 16 |
+#' @param impList \code{list} of imputation lists
|
|
| 17 |
+#' @param filterMin minimum number of peaks within a merged peak to be kept in |
|
| 18 |
+#' the analysis |
|
| 19 |
+#' @param gap gap parameter |
|
| 20 |
+#' @param D retention time penalty parameter |
|
| 21 |
+#' @param usePeaks logical, whether to use peaks (if \code{TRUE}) or the full
|
|
| 22 |
+#' 2D profile alignment (if \code{FALSE})
|
|
| 23 |
+#' @param df distance from diagonal to calculate similarity |
|
| 24 |
+#' @param verbose logical, whether to print information |
|
| 25 |
+#' @param metric numeric, different algorithm to calculate the similarity |
|
| 26 |
+#' matrix between two mass spectrum. \code{metric=1} call
|
|
| 27 |
+#' \code{normDotProduct()}; \code{metric=2} call \code{ndpRT()};
|
|
| 28 |
+#' \code{metric=3} call \code{corPrt()}
|
|
| 29 |
+#' @param type numeric, two different type of alignment function |
|
| 30 |
+#' @param penality penalization applied to the matching between two mass |
|
| 31 |
+#' spectra if \code{(t1-t2)>D}
|
|
| 32 |
+#' @param compress logical whether to compress the similarity matrix into a |
|
| 33 |
+#' sparse format. |
|
| 34 |
+#' @return \code{betweenAlignment} object
|
|
| 35 |
+#' @author Mark Robinson |
|
| 36 |
+#' @seealso \code{\link{multipleAlignment}}
|
|
| 37 |
+#' @references Mark D Robinson (2008). Methods for the analysis of gas |
|
| 38 |
+#' chromatography - mass spectrometry data \emph{PhD dissertation} University
|
|
| 39 |
+#' of Melbourne. |
|
| 40 |
+#' @keywords classes |
|
| 41 |
+#' @examples |
|
| 42 |
+#' |
|
| 43 |
+#' require(gcspikelite) |
|
| 44 |
+#' ## see 'multipleAlignment' |
|
| 45 |
+#' @importFrom stats median |
|
| 46 |
+#' @export betweenAlignment |
|
| 2 | 47 |
betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 3 | 48 |
gap = 0.7, D = 10, usePeaks = TRUE, df = 30, |
| 4 |
- verbose = TRUE, metric = 2, type = 2, penality = 0.2){
|
|
| 49 |
+ verbose = TRUE, metric = 2, type = 2, |
|
| 50 |
+ penality = 0.2, compress = FALSE){
|
|
| 5 | 51 |
n <- length(pAList) |
| 6 | 52 |
if(length(filterMin) == 1) |
| 7 | 53 |
{
|
| ... | ... |
@@ -61,9 +107,11 @@ betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 61 | 107 |
## return(wRA) |
| 62 | 108 |
## filtind <- |
| 63 | 109 |
## RR |
| 64 |
- cA <- clusterAlignment(wRA, runs = 1:n, gap = gap, D = D, df = df, metric = metric, type = type, |
|
| 65 |
- compress = TRUE, penality = penality) ## bug here with filterMin > 1 |
|
| 66 |
- pA <- progressiveAlignment(wRA, cA, gap = gap, D = D, df = df, compress = TRUE, type = type) |
|
| 110 |
+ cA <- clusterAlignment(wRA, runs = 1:n, gap = gap, D = D, df = df, |
|
| 111 |
+ metric = metric, type = type, |
|
| 112 |
+ compress = compress, penality = penality) ## bug here with filterMin > 1 |
|
| 113 |
+ pA <- progressiveAlignment(wRA, cA, gap = gap, D = D, df = df, |
|
| 114 |
+ compress = compress, type = type) |
|
| 67 | 115 |
## |
| 68 | 116 |
## cA <- clusterAlignment(wRA, 1:n, gap = gap, D = D, df = df, verbose = verbose) |
| 69 | 117 |
## pA <- progressiveAlignment(wRA, cA, gap = gap, D = D, df = df, verbose = verbose) |
| ... | ... |
@@ -79,7 +127,7 @@ betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 79 | 127 |
full.newind[[i]] <- matrix(NA, nrow(ind), length(pD@rawdata)) |
| 80 | 128 |
} |
| 81 | 129 |
names(full.newind) <- c("apex", "start", "end")
|
| 82 |
- |
|
| 130 |
+ |
|
| 83 | 131 |
col <- 1 |
| 84 | 132 |
for(i in 1:length(groups)){
|
| 85 | 133 |
g <- groups[i] |
| ... | ... |
@@ -99,7 +147,7 @@ betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 99 | 147 |
} |
| 100 | 148 |
col <- col + ncol(this.ind) |
| 101 | 149 |
} |
| 102 |
- |
|
| 150 |
+ |
|
| 103 | 151 |
new("betweenAlignment",
|
| 104 | 152 |
mergedPeaksDataset = wRA, |
| 105 | 153 |
ind = full.ind, |
| ... | ... |
@@ -113,10 +161,25 @@ betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
| 113 | 161 |
) |
| 114 | 162 |
} |
| 115 | 163 |
|
| 164 |
+ |
|
| 165 |
+#' Show method for "between" alignment of many GCMS samples |
|
| 166 |
+#' |
|
| 167 |
+#' This function show the results of a "between" alignment |
|
| 168 |
+#' |
|
| 169 |
+#' \code{betweenAlignment} objects gives the data structure which stores the
|
|
| 170 |
+#' result of an alignment across several "pseudo" datasets. These pseudo |
|
| 171 |
+#' datasets are constructed by merging the "within" alignments. |
|
| 172 |
+#' @title betweenAlignment-show |
|
| 173 |
+#' @param object |
|
| 174 |
+#' @return \code{betweenAlignment} object
|
|
| 175 |
+#' @author Mark Robinson |
|
| 176 |
+#' @export |
|
| 177 |
+#' @noRd |
|
| 116 | 178 |
setMethod("show","betweenAlignment",
|
| 117 | 179 |
function(object){
|
| 118 |
- cat("An object of class \"", class(object), "\"\n", sep = "")
|
|
| 119 |
- cat(length(object@mergedPeaksDataset@peaksrt), "groups:", |
|
| 120 |
- sapply(object@mergedPeaksDataset@peaksrt, length), "merged peaks\n" |
|
| 121 |
- ) |
|
| 122 |
-}) |
|
| 180 |
+ cat("An object of class \"", class(object), "\"\n", sep = "")
|
|
| 181 |
+ cat(length(object@mergedPeaksDataset@peaksrt), "groups:", |
|
| 182 |
+ sapply(object@mergedPeaksDataset@peaksrt, length), "merged peaks\n" |
|
| 183 |
+ ) |
|
| 184 |
+ } |
|
| 185 |
+ ) |
Ricca authored on 10/10/2017 08:30:53
| ... | ... |
@@ -1,85 +1,122 @@ |
| 1 | 1 |
|
| 2 |
+betweenAlignment <- function(pD, cAList, pAList, impList, filterMin = 1, |
|
| 3 |
+ gap = 0.7, D = 10, usePeaks = TRUE, df = 30, |
|
| 4 |
+ verbose = TRUE, metric = 2, type = 2, penality = 0.2){
|
|
| 5 |
+ n <- length(pAList) |
|
| 6 |
+ if(length(filterMin) == 1) |
|
| 7 |
+ {
|
|
| 8 |
+ filterMin <- rep(filterMin, n) |
|
| 9 |
+ } |
|
| 10 |
+ pkd <- vector("list", n)
|
|
| 11 |
+ pkr <- vector("list", n)
|
|
| 12 |
+ filtind <- vector("list", n) # list of filtered indices
|
|
| 13 |
+ newind <- vector("list", n)
|
|
| 14 |
+ for(g in 1:n){
|
|
| 15 |
+ pa <- pAList[[g]] |
|
| 16 |
+ nm <- length(pa@merges) |
|
| 17 |
+ runs <- pa@merges[[nm]]$runs |
|
| 18 |
+ ind <- pa@merges[[nm]]$ind |
|
| 19 |
+ keep <- rowSums(!is.na(ind)) >= filterMin[g] ## RR |
|
| 20 |
+ if (verbose) |
|
| 21 |
+ {
|
|
| 22 |
+ cat("[betweenAlignment]", names(pAList)[g], ": Removing", nrow(ind) - sum(keep),"peaks.\n")
|
|
| 23 |
+ ind <- ind[keep,] |
|
| 24 |
+ } |
|
| 25 |
+ newind[[g]] <- lapply(impList[[g]], FUN = function(u,k) u[k,], k = keep) |
|
| 26 |
+ filtind[[g]] <- ind |
|
| 27 |
+ rt <- matrix(NA, nrow(ind), ncol(ind)) |
|
| 28 |
+ mz <- pD@mz |
|
| 29 |
+ peaksdata <- matrix(0, nrow = length(pD@mz), ncol = nrow(ind)) |
|
| 30 |
+ for(j in 1:ncol(ind)){
|
|
| 31 |
+ if(usePeaks) |
|
| 32 |
+ {
|
|
| 33 |
+ cur.ds <- pD@peaksdata[[runs[j]]] |
|
| 34 |
+ cur.rt <- pD@peaksrt |
|
| 35 |
+ } |
|
| 36 |
+ else |
|
| 37 |
+ {
|
|
| 38 |
+ cur.ds <- pD@rawdata[[runs[j]]] |
|
| 39 |
+ cur.rt <- pD@rawrt |
|
| 40 |
+ } |
|
| 41 |
+ for(i in 1:nrow(ind)){
|
|
| 42 |
+ if(!is.na(ind[i, j])) |
|
| 43 |
+ {
|
|
| 44 |
+ peaksdata[,i] <- peaksdata[,i] + cur.ds[,ind[i,j]] |
|
| 45 |
+ rt[i,j] <- cur.rt[[runs[j]]][ind[i,j]] |
|
| 46 |
+ } |
|
| 47 |
+ } |
|
| 48 |
+ } |
|
| 49 |
+ pkd[[g]] <- peaksdata |
|
| 50 |
+ if(!usePeaks) |
|
| 51 |
+ {
|
|
| 52 |
+ pkd[[g]] <- pkd[[g]] / outer(rep(1, nrow(peaksdata)), rowSums(!is.na(ind))) |
|
| 53 |
+ } # average over the number of samples |
|
| 54 |
+ pkr[[g]] <- apply(rt, 1, median, na.rm = TRUE) |
|
| 55 |
+ } |
|
| 56 |
+ # wRA <- new("peaksDataset", mz = mz, rawdata = NULL, files = names(cAList),
|
|
| 57 |
+ # rawrt = NULL, peaksdata = pkd, peaksrt = pkr, filtind = filtind) |
|
| 58 |
+ wRA <- new("peaksDataset",
|
|
| 59 |
+ mz = mz, rawdata = list(NULL), files = names(cAList), |
|
| 60 |
+ rawrt = list(NULL), peaksdata = pkd, peaksrt = pkr) |
|
| 61 |
+ ## return(wRA) |
|
| 62 |
+ ## filtind <- |
|
| 63 |
+ ## RR |
|
| 64 |
+ cA <- clusterAlignment(wRA, runs = 1:n, gap = gap, D = D, df = df, metric = metric, type = type, |
|
| 65 |
+ compress = TRUE, penality = penality) ## bug here with filterMin > 1 |
|
| 66 |
+ pA <- progressiveAlignment(wRA, cA, gap = gap, D = D, df = df, compress = TRUE, type = type) |
|
| 67 |
+ ## |
|
| 68 |
+ ## cA <- clusterAlignment(wRA, 1:n, gap = gap, D = D, df = df, verbose = verbose) |
|
| 69 |
+ ## pA <- progressiveAlignment(wRA, cA, gap = gap, D = D, df = df, verbose = verbose) |
|
| 2 | 70 |
|
| 3 |
-betweenAlignment<-function(pD,cAList,pAList,impList,filterMin=3,gap=0.7,D=10,usePeaks=TRUE,df=30,verbose=TRUE) {
|
|
| 4 |
- n<-length(pAList) |
|
| 5 |
- if(length(filterMin)==1) filterMin<-rep(filterMin,n) |
|
| 6 |
- pkd<-vector("list",n)
|
|
| 7 |
- pkr<-vector("list",n)
|
|
| 8 |
- filtind<-vector("list",n) # list of filtered indices
|
|
| 9 |
- newind<-vector("list",n)
|
|
| 10 |
- for(g in 1:n) {
|
|
| 11 |
- pa<-pAList[[g]] |
|
| 12 |
- nm<-length(pa@merges) |
|
| 13 |
- runs<-pa@merges[[nm]]$runs |
|
| 14 |
- ind<-pa@merges[[nm]]$ind |
|
| 15 |
- keep<-rowSums(!is.na(ind))>filterMin[g] |
|
| 16 |
- if (verbose) |
|
| 17 |
- cat("[betweenAlignment]", names(pAList)[g], ": Removing", nrow(ind)-sum(keep),"peaks.\n")
|
|
| 18 |
- ind<-ind[keep,] |
|
| 19 |
- newind[[g]]<-lapply(impList[[g]],FUN=function(u,k) u[k,],k=keep) |
|
| 20 |
- filtind[[g]]<-ind |
|
| 21 |
- rt<-matrix(NA,nrow(ind),ncol(ind)) |
|
| 22 |
- mz<-pD@mz |
|
| 23 |
- peaksdata<-matrix(0,nrow=length(pD@mz),ncol=nrow(ind)) |
|
| 24 |
- for(j in 1:ncol(ind)) {
|
|
| 25 |
- if (usePeaks) {
|
|
| 26 |
- cur.ds<-pD@peaksdata[[runs[j]]] |
|
| 27 |
- cur.rt<-pD@peaksrt |
|
| 28 |
- } else {
|
|
| 29 |
- cur.ds<-pD@rawdata[[runs[j]]] |
|
| 30 |
- cur.rt<-pD@rawrt |
|
| 31 |
- } |
|
| 32 |
- for(i in 1:nrow(ind)) {
|
|
| 33 |
- if(!is.na(ind[i,j])) {
|
|
| 34 |
- peaksdata[,i]<-peaksdata[,i]+cur.ds[,ind[i,j]] |
|
| 35 |
- rt[i,j]<-cur.rt[[runs[j]]][ind[i,j]] |
|
| 36 |
- } |
|
| 37 |
- } |
|
| 71 |
+ nm <- length(pA@merges) |
|
| 72 |
+ ind <- pA@merges[[nm]]$ind |
|
| 73 |
+ groups <- pA@merges[[nm]]$runs |
|
| 74 |
+ full.runs <- NULL |
|
| 75 |
+ full.groups <- NULL |
|
| 76 |
+ full.ind <- matrix(NA, nrow(ind), length(pD@rawdata)) |
|
| 77 |
+ full.newind <- vector("list", 3)
|
|
| 78 |
+ for(i in 1:length(full.newind)){
|
|
| 79 |
+ full.newind[[i]] <- matrix(NA, nrow(ind), length(pD@rawdata)) |
|
| 80 |
+ } |
|
| 81 |
+ names(full.newind) <- c("apex", "start", "end")
|
|
| 82 |
+ |
|
| 83 |
+ col <- 1 |
|
| 84 |
+ for(i in 1:length(groups)){
|
|
| 85 |
+ g <- groups[i] |
|
| 86 |
+ n <- length(pAList[[g]]@merges) |
|
| 87 |
+ full.runs <- c(full.runs, pAList[[g]]@merges[[n]]$runs) |
|
| 88 |
+ nr <- length(pAList[[g]]@merges[[n]]$runs) |
|
| 89 |
+ full.groups <- c(full.groups, rep(wRA@files[g],nr)) |
|
| 90 |
+ this.ind <- filtind[[g]] |
|
| 91 |
+ rw <- which(!is.na(ind[,i])) |
|
| 92 |
+ cl <- col:(col + ncol(this.ind) - 1) |
|
| 93 |
+ full.ind[rw, cl] <- this.ind |
|
| 94 |
+ if(!is.null(newind[[g]]$apex)) |
|
| 95 |
+ {
|
|
| 96 |
+ for(j in 1:length(full.newind)){
|
|
| 97 |
+ full.newind[[j]] <- newind[[g]][[j]] |
|
| 98 |
+ } |
|
| 99 |
+ } |
|
| 100 |
+ col <- col + ncol(this.ind) |
|
| 38 | 101 |
} |
| 39 |
- pkd[[g]]<-peaksdata |
|
| 40 |
- if (!usePeaks) |
|
| 41 |
- pkd[[g]]<-pkd[[g]] / outer(rep(1,nrow(peaksdata)),rowSums(!is.na(ind))) # average over the number of samples |
|
| 42 |
- pkr[[g]]<-apply(rt,1,median,na.rm=TRUE) |
|
| 43 |
- } |
|
| 44 |
- #wRA<-new("peaksDataset",mz=mz,rawdata=NULL,files=names(cAList),rawrt=NULL,peaksdata=pkd,peaksrt=pkr,filtind=filtind)
|
|
| 45 |
- wRA<-new("peaksDataset",mz=mz,rawdata=list(NULL),files=names(cAList),rawrt=list(NULL),peaksdata=pkd,peaksrt=pkr)
|
|
| 46 |
- #return(wRA) |
|
| 47 |
- #filtind<- |
|
| 48 |
- cA<-clusterAlignment(wRA,1:n,gap=gap,D=D,df=df,verbose=verbose) |
|
| 49 |
- pA<-progressiveAlignment(wRA,cA,gap=gap,D=D,df=df,verbose=verbose) |
|
| 50 |
- nm<-length(pA@merges) |
|
| 51 |
- ind<-pA@merges[[nm]]$ind |
|
| 52 |
- groups<-pA@merges[[nm]]$runs |
|
| 53 |
- full.runs<-NULL |
|
| 54 |
- full.groups<-NULL |
|
| 55 |
- full.ind<-matrix(NA,nrow(ind),length(pD@rawdata)) |
|
| 56 |
- full.newind<-vector("list",3)
|
|
| 57 |
- for(i in 1:length(full.newind)) |
|
| 58 |
- full.newind[[i]]<-matrix(NA,nrow(ind),length(pD@rawdata)) |
|
| 59 |
- names(full.newind)<-c("apex","start","end")
|
|
| 60 |
- |
|
| 61 |
- col<-1 |
|
| 62 |
- for(i in 1:length(groups)) {
|
|
| 63 |
- g<-groups[i] |
|
| 64 |
- n<-length(pAList[[g]]@merges) |
|
| 65 |
- full.runs<-c(full.runs,pAList[[g]]@merges[[n]]$runs) |
|
| 66 |
- nr<-length(pAList[[g]]@merges[[n]]$runs) |
|
| 67 |
- full.groups<-c(full.groups,rep(wRA@files[g],nr)) |
|
| 68 |
- this.ind<-filtind[[g]] |
|
| 69 |
- rw<-which(!is.na(ind[,i])) |
|
| 70 |
- cl<-col:(col+ncol(this.ind)-1) |
|
| 71 |
- full.ind[rw,cl]<-this.ind |
|
| 72 |
- if(!is.null(newind[[g]]$apex)) {
|
|
| 73 |
- for(j in 1:length(full.newind)) |
|
| 74 |
- full.newind[[j]]<-newind[[g]][[j]] |
|
| 75 |
- } |
|
| 76 |
- col<-col+ncol(this.ind) |
|
| 77 |
- } |
|
| 78 |
- new("betweenAlignment",mergedPeaksDataset=wRA,ind=full.ind,imputeind=full.newind,runs=full.runs,groups=full.groups,cA=cA,pA=pA,filtind=filtind,newind=newind)
|
|
| 102 |
+ |
|
| 103 |
+ new("betweenAlignment",
|
|
| 104 |
+ mergedPeaksDataset = wRA, |
|
| 105 |
+ ind = full.ind, |
|
| 106 |
+ imputeind = full.newind, |
|
| 107 |
+ runs = full.runs, |
|
| 108 |
+ groups = full.groups, |
|
| 109 |
+ cA = cA, |
|
| 110 |
+ pA = pA, |
|
| 111 |
+ filtind = filtind, |
|
| 112 |
+ newind = newind |
|
| 113 |
+ ) |
|
| 79 | 114 |
} |
| 80 | 115 |
|
| 81 | 116 |
setMethod("show","betweenAlignment",
|
| 82 |
-function(object) {
|
|
| 83 |
- cat("An object of class \"",class(object),"\"\n",sep="")
|
|
| 84 |
- cat(length(object@mergedPeaksDataset@peaksrt),"groups:",sapply(object@mergedPeaksDataset@peaksrt,length),"merged peaks\n") |
|
| 117 |
+ function(object){
|
|
| 118 |
+ cat("An object of class \"", class(object), "\"\n", sep = "")
|
|
| 119 |
+ cat(length(object@mergedPeaksDataset@peaksrt), "groups:", |
|
| 120 |
+ sapply(object@mergedPeaksDataset@peaksrt, length), "merged peaks\n" |
|
| 121 |
+ ) |
|
| 85 | 122 |
}) |
git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/flagme@114594 bc3139a8-67e5-0310-9ffc-ced21a209358
Riccardo Romoli authored on 10/03/2016 11:16:22
git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/flagme@64110 bc3139a8-67e5-0310-9ffc-ced21a209358
Mark Robinson authored on 19/03/2012 18:39:46
| ... | ... |
@@ -39,7 +39,7 @@ betweenAlignment<-function(pD,cAList,pAList,impList,filterMin=3,gap=0.7,D=10,use |
| 39 | 39 |
pkd[[g]]<-peaksdata |
| 40 | 40 |
if (!usePeaks) |
| 41 | 41 |
pkd[[g]]<-pkd[[g]] / outer(rep(1,nrow(peaksdata)),rowSums(!is.na(ind))) # average over the number of samples |
| 42 |
- pkr[[g]]<-apply(rt,1,median,na.rm=T) |
|
| 42 |
+ pkr[[g]]<-apply(rt,1,median,na.rm=TRUE) |
|
| 43 | 43 |
} |
| 44 | 44 |
#wRA<-new("peaksDataset",mz=mz,rawdata=NULL,files=names(cAList),rawrt=NULL,peaksdata=pkd,peaksrt=pkr,filtind=filtind)
|
| 45 | 45 |
wRA<-new("peaksDataset",mz=mz,rawdata=list(NULL),files=names(cAList),rawrt=list(NULL),peaksdata=pkd,peaksrt=pkr)
|
| ... | ... |
@@ -82,4 +82,4 @@ setMethod("show","betweenAlignment",
|
| 82 | 82 |
function(object) {
|
| 83 | 83 |
cat("An object of class \"",class(object),"\"\n",sep="")
|
| 84 | 84 |
cat(length(object@mergedPeaksDataset@peaksrt),"groups:",sapply(object@mergedPeaksDataset@peaksrt,length),"merged peaks\n") |
| 85 |
-}) |
|
| 86 | 85 |
\ No newline at end of file |
| 86 |
+}) |
git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/flagme@34957 bc3139a8-67e5-0310-9ffc-ced21a209358
Marc Carlson authored on 28/10/2008 20:36:46
| 1 | 1 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,85 @@ |
| 1 |
+ |
|
| 2 |
+ |
|
| 3 |
+betweenAlignment<-function(pD,cAList,pAList,impList,filterMin=3,gap=0.7,D=10,usePeaks=TRUE,df=30,verbose=TRUE) {
|
|
| 4 |
+ n<-length(pAList) |
|
| 5 |
+ if(length(filterMin)==1) filterMin<-rep(filterMin,n) |
|
| 6 |
+ pkd<-vector("list",n)
|
|
| 7 |
+ pkr<-vector("list",n)
|
|
| 8 |
+ filtind<-vector("list",n) # list of filtered indices
|
|
| 9 |
+ newind<-vector("list",n)
|
|
| 10 |
+ for(g in 1:n) {
|
|
| 11 |
+ pa<-pAList[[g]] |
|
| 12 |
+ nm<-length(pa@merges) |
|
| 13 |
+ runs<-pa@merges[[nm]]$runs |
|
| 14 |
+ ind<-pa@merges[[nm]]$ind |
|
| 15 |
+ keep<-rowSums(!is.na(ind))>filterMin[g] |
|
| 16 |
+ if (verbose) |
|
| 17 |
+ cat("[betweenAlignment]", names(pAList)[g], ": Removing", nrow(ind)-sum(keep),"peaks.\n")
|
|
| 18 |
+ ind<-ind[keep,] |
|
| 19 |
+ newind[[g]]<-lapply(impList[[g]],FUN=function(u,k) u[k,],k=keep) |
|
| 20 |
+ filtind[[g]]<-ind |
|
| 21 |
+ rt<-matrix(NA,nrow(ind),ncol(ind)) |
|
| 22 |
+ mz<-pD@mz |
|
| 23 |
+ peaksdata<-matrix(0,nrow=length(pD@mz),ncol=nrow(ind)) |
|
| 24 |
+ for(j in 1:ncol(ind)) {
|
|
| 25 |
+ if (usePeaks) {
|
|
| 26 |
+ cur.ds<-pD@peaksdata[[runs[j]]] |
|
| 27 |
+ cur.rt<-pD@peaksrt |
|
| 28 |
+ } else {
|
|
| 29 |
+ cur.ds<-pD@rawdata[[runs[j]]] |
|
| 30 |
+ cur.rt<-pD@rawrt |
|
| 31 |
+ } |
|
| 32 |
+ for(i in 1:nrow(ind)) {
|
|
| 33 |
+ if(!is.na(ind[i,j])) {
|
|
| 34 |
+ peaksdata[,i]<-peaksdata[,i]+cur.ds[,ind[i,j]] |
|
| 35 |
+ rt[i,j]<-cur.rt[[runs[j]]][ind[i,j]] |
|
| 36 |
+ } |
|
| 37 |
+ } |
|
| 38 |
+ } |
|
| 39 |
+ pkd[[g]]<-peaksdata |
|
| 40 |
+ if (!usePeaks) |
|
| 41 |
+ pkd[[g]]<-pkd[[g]] / outer(rep(1,nrow(peaksdata)),rowSums(!is.na(ind))) # average over the number of samples |
|
| 42 |
+ pkr[[g]]<-apply(rt,1,median,na.rm=T) |
|
| 43 |
+ } |
|
| 44 |
+ #wRA<-new("peaksDataset",mz=mz,rawdata=NULL,files=names(cAList),rawrt=NULL,peaksdata=pkd,peaksrt=pkr,filtind=filtind)
|
|
| 45 |
+ wRA<-new("peaksDataset",mz=mz,rawdata=list(NULL),files=names(cAList),rawrt=list(NULL),peaksdata=pkd,peaksrt=pkr)
|
|
| 46 |
+ #return(wRA) |
|
| 47 |
+ #filtind<- |
|
| 48 |
+ cA<-clusterAlignment(wRA,1:n,gap=gap,D=D,df=df,verbose=verbose) |
|
| 49 |
+ pA<-progressiveAlignment(wRA,cA,gap=gap,D=D,df=df,verbose=verbose) |
|
| 50 |
+ nm<-length(pA@merges) |
|
| 51 |
+ ind<-pA@merges[[nm]]$ind |
|
| 52 |
+ groups<-pA@merges[[nm]]$runs |
|
| 53 |
+ full.runs<-NULL |
|
| 54 |
+ full.groups<-NULL |
|
| 55 |
+ full.ind<-matrix(NA,nrow(ind),length(pD@rawdata)) |
|
| 56 |
+ full.newind<-vector("list",3)
|
|
| 57 |
+ for(i in 1:length(full.newind)) |
|
| 58 |
+ full.newind[[i]]<-matrix(NA,nrow(ind),length(pD@rawdata)) |
|
| 59 |
+ names(full.newind)<-c("apex","start","end")
|
|
| 60 |
+ |
|
| 61 |
+ col<-1 |
|
| 62 |
+ for(i in 1:length(groups)) {
|
|
| 63 |
+ g<-groups[i] |
|
| 64 |
+ n<-length(pAList[[g]]@merges) |
|
| 65 |
+ full.runs<-c(full.runs,pAList[[g]]@merges[[n]]$runs) |
|
| 66 |
+ nr<-length(pAList[[g]]@merges[[n]]$runs) |
|
| 67 |
+ full.groups<-c(full.groups,rep(wRA@files[g],nr)) |
|
| 68 |
+ this.ind<-filtind[[g]] |
|
| 69 |
+ rw<-which(!is.na(ind[,i])) |
|
| 70 |
+ cl<-col:(col+ncol(this.ind)-1) |
|
| 71 |
+ full.ind[rw,cl]<-this.ind |
|
| 72 |
+ if(!is.null(newind[[g]]$apex)) {
|
|
| 73 |
+ for(j in 1:length(full.newind)) |
|
| 74 |
+ full.newind[[j]]<-newind[[g]][[j]] |
|
| 75 |
+ } |
|
| 76 |
+ col<-col+ncol(this.ind) |
|
| 77 |
+ } |
|
| 78 |
+ new("betweenAlignment",mergedPeaksDataset=wRA,ind=full.ind,imputeind=full.newind,runs=full.runs,groups=full.groups,cA=cA,pA=pA,filtind=filtind,newind=newind)
|
|
| 79 |
+} |
|
| 80 |
+ |
|
| 81 |
+setMethod("show","betweenAlignment",
|
|
| 82 |
+function(object) {
|
|
| 83 |
+ cat("An object of class \"",class(object),"\"\n",sep="")
|
|
| 84 |
+ cat(length(object@mergedPeaksDataset@peaksrt),"groups:",sapply(object@mergedPeaksDataset@peaksrt,length),"merged peaks\n") |
|
| 85 |
+}) |
|
| 0 | 86 |
\ No newline at end of file |