Hi Martin,
there appears to be a minor bug in one of the qa() plots in the code;

ShortRead:::.plotReadOccurrences

If the code chunk which generates the dataset for plotting the read
occurences is run manually on an example dataset, the cumulative
proportion of reads can be slightly greater than 1. It looks like a slight
error in the denominator of the calculation.

library(ShortRead)

exptPath <- system.file("extdata", package = "ShortRead")
sp       <- SolexaPath(exptPath)
qa       <- qa(sp)

df       <- qa[["sequenceDistribution"]]

tmp <- with(df, {
        nOccur <- tapply(nOccurrences, lane, c)
        cumulative <- tapply(nOccurrences * nReads, lane, function(elt) {
            cs <- cumsum(elt)
            (cs - cs[1] + 1)/max(1, diff(range(cs))) ## denominator under
estimates by 1
        })
        lane <- tapply(lane, lane, c)
        data.frame(nOccurrences = unlist(nOccur), cumulative =
unlist(cumulative),
            lane = unlist(lane))
    })

print(tmp)
                 nOccurrences   cumulative lane
s_2_export.txt1             1 0.0008347245    1
s_2_export.txt2             2 0.0091819699    1
s_2_export.txt3             3 0.0116861436    1
s_2_export.txt4             5 0.0158597663    1
s_2_export.txt5            10 0.0242070117    1
s_2_export.txt6             1 0.6435726210    1
s_2_export.txt7             2 0.6519198664    1
s_2_export.txt8             3 0.6544240401    1
s_2_export.txt9             9 0.6619365609    1
s_2_export.txt10            1 0.9991652755    1
s_2_export.txt11            2 1.0008347245    1

tmp[which(tmp$cumulative>1),]
                 nOccurrences cumulative lane
s_2_export.txt11            2   1.000835    1

tmp <- with(df, {
        nOccur <- tapply(nOccurrences, lane, c)
        cumulative <- tapply(nOccurrences * nReads, lane, function(elt) {
            cs <- cumsum(elt)
            (cs - cs[1] + 1)/max(1, (diff(range(cs))+1))  ## Add 1 to
difference
        })
        lane <- tapply(as.character(lane), lane, c) ## as.character to stop
factor coercing
       data.frame(nOccurrences = unlist(nOccur), cumulative =
unlist(cumulative),
            lane = unlist(lane), row.names=NULL)  ## Removing row names
    })

print(tmp)
   nOccurrences   cumulative           lane
1             1 0.0008340284 s_2_export.txt
2             2 0.0091743119 s_2_export.txt
3             3 0.0116763970 s_2_export.txt
4             5 0.0158465388 s_2_export.txt
5            10 0.0241868224 s_2_export.txt
6             1 0.6430358632 s_2_export.txt
7             2 0.6513761468 s_2_export.txt
8             3 0.6538782319 s_2_export.txt
9             9 0.6613844871 s_2_export.txt
10            1 0.9983319433 s_2_export.txt
11            2 1.0000000000 s_2_export.txt



One suggestion is that the first part of this code (or similar) would be
useful  in the
help so example(qa) could be run and testing with reporoducible examples can
be quicker.

cheers,

Marcus


sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=English_New Zealand.1252  LC_CTYPE=English_New Zealand.1252

[3] LC_MONETARY=English_New Zealand.1252 LC_NUMERIC=C

[5] LC_TIME=English_New Zealand.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ShortRead_1.10.3    Rsamtools_1.4.1     lattice_0.19-23
[4] Biostrings_2.20.1   GenomicRanges_1.4.3 IRanges_1.10.3

loaded via a namespace (and not attached):
[1] Biobase_2.12.1 grid_2.13.0    hwriter_1.3

        [[alternative HTML version deleted]]

_______________________________________________
Bioc-sig-sequencing mailing list
Bioc-sig-sequencing@r-project.org
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

Reply via email to