[Hparker]

Sorry the error message is actually this: (I changed the number of samples >=2 to 1 briefly to check whether this was the problem)

Error in which(rowSums(BS.cov[, bsseq.data$Type == "SSA"] >= 2) >= 2 & :
error in evaluating the argument 'x' in selecting a method for function 'which': Error in rowSums(BS.cov[, bsseq.data$Type == "SSA"] >= 2) :
'x' must be an array of at least two dimensions

In addition when I look at the number of CpGs covered by at least 1 read in all samples I get only 151

sum(rowSums(getCoverage(bsseq.data) >= 1) == 14
[151]

[Hparker]

I also get an error:

Error in which(rowSums(BS.cov[, bsseq.data$Type == "SSA", drop = F] >= :
error in evaluating the argument 'x' in selecting a method for function 'which': Error in rowSums(BS.cov[, bsseq.data$Type == "SSA", drop = F] >= 2) :
'x' must be an array of at least two dimensions

[dpryan]

Oh, that's the problem then. You don't seem to have enough coverage of enough sites to do anything really. Perhaps it'd help if you described the underlying biological experiment and how the sequencing data was processed. It's likely that something went amiss earlier in the process.

[Hparker]

The approach used Zymo's Fluidigm's array to simultaneously bisulfite sequence 48 samples targeting 48 genetic loci. I used a heterogenous mix of samples, mostly DNA from fresh tissue but also some FFPE samples.

How many reads should I be getting to make the analysis? It could be that only one sample which amplified poorly is affecting this? Or is it more likely to be accross samples?

Is it possible for example to remove a genetic loci region from the analysis and if so how could I then do this in R?

[dpryan]

That method was designed with something like 5x average coverage in mind (I think that's what they used in the paper, but you'd have to double check that). With targeted BSseq this would be pretty simple to achieve (5x is probably an absolute minimum too).

Yes, you can subset the objects however you like. Just "obj <- obj[idx,]", where "idx" contains the rows you want to keep. See help(SummarizedExperiment) for further information.

[Hparker]

So what do you think the problem is when I try to use the: keepLoci (below) ?

[dpryan]

No clue, I'd have to play with your dataset. Perhaps post something on the bioconductor support site and see if the BSseq authors reply.

[Hparker]

Ok will post on bioconductor site, thank you.

Just a couple more questions:

- when i use getCoverage and then View this information I have the samples (V1, V2...) as columns and CpG sites as rows. How can I rename the samples as the actual samples not V1 etc. and also label the CpG sites as specific regions? Because at the moment I can't tell which samples correspond to areas of poor coverage.
The numbers that are shown in this table for the coverage are not the actual number of reads, but is an algorithm used to determine these numbers based on the number of reads? Is there a cutoff for what might be considered low read coverage?

Many thanks

[dpryan]

No clue and I don't have the time to check at the moment.