[aroma.affymetrix] Gene ST array summarisation by probeset?

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

[aroma.affymetrix] Gene ST array summarisation by probeset?

Sophie Marion de Procé
Dear all,

I'm analysing a Rat Gene ST 2.1 array. I would like to filter the dataset using thresholds of expression for a minimum proportion of samples in each group.
I've been following the paper by Rodrigo-Domingo et al. (2014) Reproducible probe-level analysis of the Affymetrix Exon 1.0 ST array with R/Bioconductor.

They have two steps to the filtering step, first they filter probesets and then the transcripts. I'm getting stuck at the filtering step (code chunk number 7: intensityFiltering in the R code provided in this paper), as I can't find a way to summarise the Gene St array at the probeset level in order to do the first step. I have used plm <- RmaPlm(can) to summarise my data at the transcript level, but the plmPs  <- RmaPlm(csN, mergeGroups = FALSE) seems to summarise by transcripts as well.

So my questions are:
1) Is it possible to summarise a Gene ST array at the probeset level? If yes, how?
2) Less specific to the aroma-affymetrix package, but is it necessary to have the probeset-level dataset in order to filter present/absent probes/transcripts? What would be an appropriate workflow for this?

Thank you very much for your help,
Best wishes,
Sophie.

Here is the code chunk from that paper:

###################################################
### code chunk number 7: intensityFiltering
###################################################

# ** user-defined groups; default: groups defined by the treatment column in SampleInformation.txt
sample
.info$number <- seq(1, nrow(sample.info))
groups
<- split(sample.info$number, sample.info$treatment)

# check whether the filtering is already performed, perform it otherwise
# at probeset level
if(file.exists(file = paste(output.folder, "/PresentProbesets", ds, ".Rdata", sep = ""))){
 load
(file = paste(output.folder, "/PresentProbesets", ds, ".Rdata", sep = ""))
} else {
 
# remove cross-hybridising probesets
 non
.crosshyb.probesets <- probesets.NetAffx.32$probeset_id[probesets.NetAffx.32$crosshyb_type == 1]
 
if(!exists("exFit")){
  load
(file= paste(output.folder, "CoreExonIntensities.Rdata", sep = "/"))
 
}
 exon
.intensities <- exFit[exFit$groupName %in% non.crosshyb.probesets, -c(3:5)]
 rm
("exFit")
 gc
()

 
# ** user-defined criteria for absent/present probesets
 
# ** criterion 1: probeset absent/present in a group of samples: present in at least half of the samples
 present
.exons <- lapply(groups, FUN = function(group){
   
if(length(group) > 1){
     apply
(log2(exon.intensities[, group + 2]) < 3, 1, sum)/length(group) < 0.5
   
} else {
     log2
(exon.intensities[, group + 2]) > 3
   
}
 
})
 present
.exons <- t(do.call(rbind, present.exons))  # convert to dataframe
 rownames
(present.exons) <- exon.intensities$groupName  # use probeset identities

 
# ** criterion 2: probeset absent/present in the dataset
 
# remove probesets not present in any of the groups
 absent
.exons <- apply(present.exons, 1, AllFalse)
 probesets
.to.keep <- absent.exons[absent.exons == FALSE]
 probesets
.to.keep <- as.factor(names(probesets.to.keep))
 n
.present.exons <- length(probesets.to.keep)
 rm
(list = c("absent.exons"))
 save
(probesets.to.keep, file = paste(output.folder, "/PresentProbesets", ds, ".Rdata", sep = ""))
}
# check whether transcript filtering has been performed; perform it if not
if(file.exists(file = paste(output.folder, "/PresentTranscripts", ds, ".Rdata", sep = ""))){
 load
(file = paste(output.folder, "/PresentTranscripts", ds, ".Rdata", sep = ""))
} else {
 
if(!exists("exon.intensities")){
  load
(file= paste(output.folder, "CoreExonIntensities.Rdata", sep = "/"))
  exon
.intensities <- exFit[exFit$groupName %in% non.crosshyb.probesets, -c(3:5)]
  rm
("exFit")
  gc
()
 
}

 
# ** user-defined criteria for absent transcripts:
 
# criterion 1: half or more of probesets of transcript present in sample
 
# --> transcript present in sample
 core
.transcripts <- unique(exon.intensities$unitName)
 
# create a list of transcript clusters present/absent in each sample
 present
.genes    <- lapply(core.transcripts, FUN = function(gene){
   apply
(log2(exon.intensities[exon.intensities$unitName == gene, -c(1:2)]) < 3, 2, sum)/
   length
(exon.intensities[exon.intensities$unitName == gene,]$groupName) < 0.5
 
})# FALSE: gene not present in sample
 names
(present.genes) <- core.transcripts
 present
.genes <- do.call(rbind, present.genes) # convert to logical matrix

 
# criterion 2: transcript present in at least half of the samples of a group
 
# --> transcript present in the group
 present
.genes.in.group <- lapply(groups, FUN = function(group){
   
if(length(group) > 1){
     apply
(present.genes[ , group], 1, sum)/length(group) >= 0.5
   
} else {
     present
.genes[ , group]
   
}
 
})
 present
.genes.in.group <- do.call(rbind, present.genes.in.group) # logical matrix
 present
.genes.in.group <- t(present.genes.in.group)
 
# keep genes only present in at least two groups
 transcripts
.to.keep <- apply(present.genes.in.group, 1, TwoOrMoreTrue)
 transcripts
.to.keep <- names(transcripts.to.keep[transcripts.to.keep == TRUE])
 save
(transcripts.to.keep, file = paste(output.folder, "/PresentTranscripts", ds, ".Rdata", sep = ""))
 rm
("exon.intensities")
 gc
()
}
 n
.present.transcript.clusters <- length(transcripts.to.keep)


--
--
When reporting problems on aroma.affymetrix, make sure 1) to run the latest version of the package, 2) to report the output of sessionInfo() and traceback(), and 3) to post a complete code example.
 
 
You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group with website http://www.aroma-project.org/.
To post to this group, send email to [hidden email]
To unsubscribe and other options, go to http://www.aroma-project.org/forum/

---
You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group.
To unsubscribe from this group and stop receiving emails from it, send an email to [hidden email].
For more options, visit https://groups.google.com/d/optout.