[Bioc-sig-seq] Bioc short read directions
mtmorgan at fhcrc.org
mtmorgan at fhcrc.org
Wed Apr 2 15:49:31 CEST 2008
Hi Loyal --
Quoting Loyal Goff <lgoff at rci.rutgers.edu>:
> This is a great start...thanks to both Martin and Herve. The speed is
> indeed impressive! I do have one question. Would it be advantageous
> to reduce the data to a unique list of read sequences, and in doing so
> both retain counts in a separate slot and reduce the matrix size? It
> seems to me this would speed everything along as well. (ie. only
> attempt to align a unique sequence once). Does anyone have a need to
> retain independent reads after a quality score cutoff?
In terms of matching, PDict does the right thing and it makes (little)
difference in terms of performance whether the reads have been made unique,
while introducing a bookkeeping nightmare for the user. Hopefully other
algorithms will be implemented to cope with duplicates effectively, if
appropriate.
I think the bookkeeping argument also weighs quite heavily in favor of dealing
with the data 'as-is'. I suspect also that quality scores will play an
increasingly important part in algorithms, so a simple cut-off will not be a
good solution. And probably the time / space 'savings' is only a fraction (10%
duplicate reads?? Maybe someone on the list has experience with this?) of the
requirements, so doesn't really change things that much.
My first inclination is to leave the data accessible in all its glory.
Martin
>
> Loyal
>
> Loyal A. Goff
>
> Rutgers Stem Cell Research Center
> Rutgers: The State University of New Jersey
> Nelson Biology Labs D-251
> 604 Allison Rd,
> Piscataway, NJ 08854
> lgoff at rci.rutgers.edu
>
> On Apr 1, 2008, at 11:20 PM, Martin Morgan wrote:
>
> > Short-readers!
> >
> > I want to take this opportunity to provide an update on the directions
> > we've initiated for short reads. It would be great to get your
> > feedback,
> > and to hear of your own efforts.
> >
> > WARNING: This software is in development, and is far from final. In
> > particular, functions in the BiostringsCinterfaceDemo package are NOT
> > meant to be final or for use in other packages; they're here to
> > demonstrate 'proof-of-concept' and to illustrate how users can access
> > the Biostrings package from C code. I'll indicate which functions are
> > from BiostringsCinterfaceDemo below. Expect a short-read 'base'
> > package
> > to materialize in the devel branch of Bioconductor in the not too
> > distant future.
> >
> > WARNING: The data used to illustrate functionality are not meant to be
> > indiciative of the quality of data produced by Solexa; they are
> > generally 'first runs' that present a number of interesting challenges
> > for interpretation.
> >
> > HARDWARE AND SOFTWARE: The following include timing and object size
> > measurements. The machine being used is fast, but we're not doing
> > anything fancy to, e.g., exploit multiple processors. The machine
> > has a
> > very large amount of memory; we used about 10 GB below, looking at
> > three
> > different data sets. The following uses the R-2-7-branch (this is
> > different from R-devel). The Biostrings and BiostringsCinterfaceDemo
> > packages are updated very regularly, so be prepared for broken or
> > outdated functions.
> >
> > Herve Pages is responsible for the clever code; I am just a scribe.
> >
> > Ok, first a convenience function to print out 'size' in megabytes,
> > 'cause objects are large!
> >
> >> mb <- function(sz) noquote(sprintf("%.1f MB", sz / (1024^2)))
> >
> >
> > * Starters...
> >
> > We load the BiostringsCinterfaceDemo, which requires Biostrings.
> > Both of
> > these need to be from the 'development' branch of Bioconductor. Both
> > are
> > changing rapidly, and should be obtained from svn and updated
> > regularly
> > (http://wiki.fhcrc.org/bioc/DeveloperPage,
> > http://wiki.fhcrc.org/bioc/SvnHowTo).
> >
> >> library(BiostringsCinterfaceDemo)
> >
> >
> > * I/O, DNAStringSet, and alphabetFrequency
> >
> > We next read in a fasta file derived from a lane of solexa reads.
> > Here's
> > what the data looks like:
> >
> >> readLines(fastaFile, 10)
> > [1] ">5_1_102_368"
> > [2] "TAAGAGGTTTAAATTTTCTTCAGGTCAGTATTCTTT"
> > [3] ">5_1_120_254"
> > [4] "TTAATTCGTAAACAAGCAGTAGTAATTCCTGCTTTT"
> > [5] ">5_1_110_385"
> > [6] "GCTAATTTGCCTACTAACCAAGAACTTGATTTCTTC"
> > [7] ">5_1_118_88"
> > [8] "GTTTGGAGTGATACTGACCGCTCTCGTGGTCGTCGC"
> > [9] ">5_1_113_327"
> > [10] "GCTTGCGTTTATGGTACGCTGGACTTTGTAGGATAC"
> >
> > This is a single lane from a Solexa training run; the data are not
> > filtered, and the run is not meant to be representative in terms of
> > quality or other characteristics. The DNA used for the reads is from
> > phage phiX-174. Here's how we read it in (countLines and
> > readSolexaFastA
> > are in BiostringsCinterfaceDemo)
> >
> >> countLines(fastaFile)
> > s_5.fasta
> > 18955056
> >> system.time({
> > + seqa <- readSolexaFastA(fastaFile)
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 67.48 2.08 69.68
> >> mb(object.size(seqa))
> > [1] 397.7 MB
> >> seqa
> > A DNAStringSet instance of length 9477528
> > width seq
> > [1] 36 TAAGAGGTTTAAATTTTCTTCAGGTCAGTATTCTTT
> > [2] 36 TTAATTCGTAAACAAGCAGTAGTAATTCCTGCTTTT
> > [3] 36 GCTAATTTGCCTACTAACCAAGAACTTGATTTCTTC
> > [4] 36 GTTTGGAGTGATACTGACCGCTCTCGTGGTCGTCGC
> > [5] 36 GCTTGCGTTTATGGTACGCTGGACTTTGTAGGATAC
> > [6] 36 TGACCCTCAGCAATCTTAAACTTCTTAGACGAATCA
> > [7] 36 GCTGGTTCTCACTTCTGTTACTCCAGCTTCTTCGGC
> > [8] 36 TTTAGGTGTCTGTAAAACAGGTGCCGAAGAAGCTGG
> > [9] 36 GGTCTGTTGAACACGACCAGAAAACTGGCCTAACGA
> > ... ... ...
> > [9477520] 36 TACGCAGTTTTGCCGTATACTCGTTGTTCTGACTCT
> > [9477521] 36 TATACCCCCCCTCCTACTTGTGCTGTTTCTCATGTT
> > [9477522] 36 CAGGTTGTTTCTGTTGGTGCTGATATTTCTTTTTTT
> > [9477523] 36 GTCTTCCTTGCTTGTCAGATTGGTCGTCTTATTACC
> > [9477524] 36 ATACGAAAGACCAGGTATATGCACAAAATGAGTTGC
> > [9477525] 36 ACCACAAACGCGCTCGTTTATGCTTGCCTCTATTAC
> > [9477526] 36 ------------------------------------
> > [9477527] 36 CCAGCAAGGAAGCCAAGATGGGAAAGGTCATGCGGC
> > [9477528] 36 CATTGTTGACCACCTACATACCACAGACGAGCACCT
> >
> > It takes just over a minute to read in the nearly 9.5 million reads.
> > The
> > reads are stored efficiently, without the overhead of R character
> > strings. The data structure (a DNAStringSet, from Biostrings and
> > therefore stable) will not copy the large data, but instead contains
> > 'views' into it.
> >
> > A basic question is about the nucleotides present in the reads.
> > alphabetFrequency (from Biostrings) scans all the sequences and
> > tallies
> > nucleotides
> >
> >> system.time({
> > + alf1 <- alphabetFrequency(seqa, collapse=TRUE, freq=TRUE)
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 0.612 0.000 0.612
> >> alf1
> > A C G T M R W S Y K
> > 0.2449 0.2119 0.2201 0.3030 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
> > V H D B N -
> > 0.0000 0.0000 0.0000 0.0000 0.0000 0.0201
> >
> > (bases are recored with IUPAC symbols, see ?IUPAC_CODE_MAP). This
> > executes very efficiently. A variant produces a matrix with rows
> > corresponding to reads and columns to bases.
> >
> >> alf <- alphabetFrequency(seqa, baseOnly=TRUE)
> >> dim(alf)
> > [1] 9477528 5
> >> head(alf)
> > A C G T other
> > [1,] 9 4 6 17 0
> > [2,] 11 6 5 14 0
> > [3,] 10 9 4 13 0
> > [4,] 4 9 12 11 0
> > [5,] 6 6 11 13 0
> > [6,] 12 10 4 10 0
> >
> > This can be remarkably useful. For instance, to select just the
> > 'clean'
> > sequences (those without ambiguous base calls), one can
> >
> >> cleanSeqs <- seqa[alf[,"other"]==0]
> >> length(seqa)
> > [1] 9477528
> >> length(cleanSeqs)
> > [1] 9207292
> >
> > This creates a new DNAStringSet with just the clean sequences. It
> > executes very quickly, because the DNAStringSet is a view into the
> > original. The memory associated with the reads themselves is not
> > copied.
> > here is the alphabetFrequency of the 'clean' reads.
> >
> >> cleanAlf <- alphabetFrequency(cleanSeqs, baseOnly=TRUE)
> >
> > Again this is very useful, for instance the distribution of GC content
> > among clean reads is
> >
> >> plot(density(rowSums(cleanAlf[,c("G", "C")]) / rowSums(cleanAlf)))
> >
> >
> > * PDict, countPDict, matchPDict
> >
> > A 'PDict' (defined in Biostrings) is a dictionary-like structure that
> > can be used for very efficient exact- and partially-exact matching
> > algorithms. To illustrate, we'll use data from about a million reads
> > of
> > the Solexa BAC cloning vector. These reads again come from an early
> > run
> > on the Solexa instrumentation here, and results should not be taken to
> > be representative of performance.
> >
> > We read and clean the sequences as above, resulting in
> >
> >> length(cleanSeqs)
> > [1] 923680
> >
> > We then create a PDict from our DNAStringSet with
> >
> >> system.time({
> > + pDict <- PDict(cleanSeqs)
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 1.09 0.00 1.10
> >> pDict
> > 923680-pattern constant width PDict object of width 25 (patterns
> > have no
> > names)
> >> mb(object.size(pDict))
> > [1] 160.4 MB
> >
> > This is created quickly. It is a larger object, but the size allows
> > fast
> > searches. Here we'll use Biostrings readFASTA to read in the
> > sequence to
> > which the data are to be aligned.
> >
> >> bac <- read.DNAStringSet(bacFile, "fasta")[[1]]
> > Read 2479 items
> >> length(bac)
> > [1] 173427
> >
> > This is a BAC clone. We'll match our pDict to the BAC subject, finding
> > all EXACT matches;
> >
> >> system.time({
> > + counts <- countPDict(pDict, bac)
> > + }, gcFirst=FALSE)
> > user system elapsed
> > 0.200 0.048 0.268
> >> length(counts)
> > [1] 923680
> >> table(counts)[1:5]/sum(table(counts))
> > counts
> > 0 1 2 3 4
> > 0.53954 0.42738 0.01136 0.00528 0.00334
> >
> > This is very fast, partly because the subject against which the
> > PDict is
> > being matched is short. A more realistic use case is to match
> > against a
> > genome. The integration with BSgenome packages is very smooth. To
> > match
> > our pDict against human chromosome 6 (where the BAC cloning vector
> > comes
> > from), we can load the appropriate BSgenome package and chromosome
> >
> >> library(BSgenome.Hsapiens.UCSC.hg18)
> >
> >> Hsapiens[["chr6"]]
> > 170899992-letter "DNAString" instance
> > seq:
> > NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
> > ...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
> >
> > (The N's represent the chromsome telomeres, whose seuqences is not
> > available). We look for exact matches
> >
> >> system.time({
> > + hcount <- countPDict(pDict, Hsapiens[["chr6"]])
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 24.2 0.0 24.2
> >> table(hcount)[1:5] / sum(table(hcount))
> > hcount
> > 0 1 2 3 4
> > 0.50466 0.32996 0.02043 0.00959 0.00761
> >
> > About 1/2 the sequences exactly match one or more locations on
> > chromosome 6. Some sequences match many times, though the reason (in
> > this case, anyway) is not too surprising:
> >
> >> max(hcount)
> > [1] 8286
> >> maxIdx = which(hcount==max(hcount))
> >> cleanSeqs[maxIdx]
> > A DNAStringSet instance of length 70
> > width seq
> > [1] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [2] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [3] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [4] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [5] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [6] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [7] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [8] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [9] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > ... ... ...
> > [62] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [63] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [64] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [65] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [66] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [67] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [68] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [69] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> > [70] 25 TTTTTTTTTTTTTTTTTTTTTTTTT
> >
> > (that these are identical can be checked with
> > unique(as.character(cleanSeqs[maxIdx])))
> >
> > The current implementation of PDict allows for a 'trusted band' of
> > nucleotides that need to match exactly, allowing for approximate
> > matches
> > in the remaining nucleotides. Here we trust the first 12 bases, and
> > allow up to 3 mismatches. Also, we use the more informative
> > matchPDict,
> > which allows us to find the positions of all matches
> >
> >> trusted <- PDict(cleanSeqs, tb.end=12)
> >> system.time({
> > + mmMatch <- matchPDict(trusted, Hsapiens[[6]], max.mismatch=3)
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 195.58 3.85 199.44
> >> table(cut(countIndex(mmMatch), c(0, 10^(0:5)), right=FALSE))
> >
> > [0,1) [1,10) [10,100) [100,1e+03) [1e+03,1e+04)
> > 377233 337806 72983 66322 60818
> > [1e+04,1e+05)
> > 8518
> >
> > Execution time increases as the stringency of the match decreases.
> > PDict
> > facilities do not yet incorporate quality scores, but filtering
> > results
> > based on quality (qualities are discussed further below) represents a
> > natural direction for development.
> >
> >
> >
> > * alphabetByCycle
> >
> > alphabetByCycle uses a small C function in BiostringsCinterfaceDemo,
> > alphabet_by_cycle. It and the read* functions in this package
> > illustrate
> > how to access DNAStringSet objects at the C level. alphabetByCycle
> > is a
> > matrix that tallies nucleotide use per cycle
> >
> >> system.time({
> > + abc <- alphabetByCycle(seqa)
> > + })
> > user system elapsed
> > 2.68 0.00 2.68
> >
> > Again this can be quite useful. For instance, we can find out the
> > number
> > of bases that were not called, as a function of cycle
> >
> >> abc["-",]
> > [1] 22181 108829 123173 180382 225091 225055 216787 208538 208881
> > [10] 213104 148936 142966 141030 148163 178747 204304 211538 211303
> > [19] 213722 211167 208021 208715 165441 158359 147110 151462 204781
> > [28] 221008 223922 221171 227622 226936 232232 236606 242387 241718
> >
> > and the number of 'T' nucleotides as a function of cycle
> >
> >> abc["T",] / colSums(abc[1:4,])
> > [1] 0.286 0.292 0.292 0.284 0.292 0.280 0.293 0.297 0.299 0.287 0.290
> > [12] 0.294 0.294 0.299 0.300 0.301 0.304 0.305 0.306 0.309 0.310 0.311
> > [23] 0.312 0.316 0.315 0.319 0.320 0.322 0.326 0.328 0.331 0.335 0.339
> > [34] 0.342 0.346 0.358
> >
> > That's quite a striking increase after cycle 25!
> >
> >
> > * Qualities
> >
> > Solexa reads have quality scores associated with each base call. These
> > are summarized in files formatted like:
> >
> >> readLines(fastqFile, n=8)
> > [1] "@HWI-EAS88_1_1_1_1001_499"
> > [2] "GGACTTTGTAGGATACCCTCGCTTTCCTTCTCCTGT"
> > [3] "+HWI-EAS88_1_1_1_1001_499"
> > [4] "]]]]]]]]]]]]Y]Y]]]]]]]]]]]]VCHVMPLAS"
> > [5] "@HWI-EAS88_1_1_1_898_392"
> > [6] "GATTTCTTACCTATTAGTGGTTGAACAGCATCGGAC"
> > [7] "+HWI-EAS88_1_1_1_898_392"
> > [8] "]]]]]]]]]]]]Y]]]]]]]]]YPV]T][PZPICCK"
> >
> > A record consists of four lines. The first and third lines are
> > identifiers (repeated), the second line is the sequence, the fourth
> > line
> > an ASCII character representing the score (']' is good, Z is better
> > than
> > A). Here we read a quality file into a data structure defined in
> > BiostringsCinterfaceDemo designed to coordinate the sequence, name,
> > and
> > quality information.
> >
> >> system.time({
> > + seqq <- readSolexaFastQ(fastqFile)
> > + }, gcFirst=TRUE)
> > user system elapsed
> > 17.533 0.417 17.969
> >> mb(object.size(seqq))
> > [1] 254.5 MB
> >> seqq
> > class: SolexaSequenceQ
> > length: 2218237
> >
> >> sequences(seqq)[1:5]
> > A DNAStringSet instance of length 5
> > width seq
> > [1] 36 GGACTTTGTAGGATACCCTCGCTTTCCTTCTCCTGT
> > [2] 36 GATTTCTTACCTATTAGTGGTTGAACAGCATCGGAC
> > [3] 36 GCGGTGGTCTATAGTGTTATTAATATCAATTTGGGT
> > [4] 36 GTTACCATGATGTTATTTCTTCATTTGGAGGTAAAA
> > [5] 36 GTATGTTTCTCCTGCTTATCACCTTCTTGAAGGCTT
> >> names(seqq)[1:5]
> > A BStringSet instance of length 5
> > width seq
> > [1] 24 HWI-EAS88_1_1_1_1001_499
> > [2] 23 HWI-EAS88_1_1_1_898_392
> > [3] 23 HWI-EAS88_1_1_1_922_465
> > [4] 23 HWI-EAS88_1_1_1_895_493
> > [5] 23 HWI-EAS88_1_1_1_953_493
> >> scores(seqq)[1:5]
> > A BStringSet instance of length 5
> > width seq
> > [1] 36 ]]]]]]]]]]]]Y]Y]]]]]]]]]]]]VCHVMPLAS
> > [2] 36 ]]]]]]]]]]]]Y]]]]]]]]]YPV]T][PZPICCK
> > [3] 36 ]]]]Y]]]]]V]T]]]]]T]]]]]V]TMJEUXEFLA
> > [4] 36 ]]]]]]]]]]]]]]]]]]]]]]T]]]]RJRZTQLOA
> > [5] 36 ]]]]]]]]]]]]]]]]]T]]]]]]]]]]MJUJVLSS
> >
> > The object returned by readSolexaFastQ is an S4 object with three
> > slots.
> > Each slot contains a XStringSet, where 'X' is DNA for sequences, and
> > 'BString' for names and scores. This represents one way of structuring
> > quality data; the S4 class coordinates subsetting, and provides
> > (read-only) accessors to the underlying objects. A likely addition to
> > this class as it matures is the inclusion of lane-specific phenotype
> > (sample) information, much as an ExpressionSet coordinates sample and
> > expression values.
> >
> > We can gain some basic insight into the sequences as before, e.g.,
> >
> >> abc <- alphabetByCycle(sequences(seqq))
> >> abc["N",]
> > [1] 0 0 0 0 0 0 0 0 0 0 0
> > [12] 0 1213 1631 1155 1240 721 418 8503 526 6493 703
> > [23] 14999 718 1623 737 243 40704 811 590 1964 961 809
> > [34] 910 477 208
> >
> > Solexa provides files with quality information after filtering reads
> > based on 'purity', a measure that precludes uncertain bases (IUPAC
> > code
> > 'N') from a user-specified region (the first 12 cycles, by default).
> >
> >> abc["T",] / colSums(abc[1:4,])
> > [1] 0.298 0.290 0.292 0.290 0.295 0.279 0.292 0.298 0.301 0.293 0.295
> > [12] 0.299 0.296 0.293 0.294 0.294 0.293 0.295 0.294 0.297 0.297 0.299
> > [23] 0.298 0.303 0.302 0.307 0.312 0.312 0.321 0.331 0.334 0.351 0.361
> > [34] 0.374 0.380 0.423
> >
> > We can also summarize quality information by cycle, using an alphabet
> > that reflects the encoded scores:
> >
> >> alphabet <- sapply(as.raw(32:93), rawToChar)
> >> abcScore <- alphabetByCycle(scores(seqq), alphabet=alphabet)
> >>
> >> rowSums(abcScore)
> > ! " # $ %
> > & '
> > 0 0 0 0 0 0
> > 0 0
> > ( ) * + ,
> > - . /
> > 0 0 0 0 0 0
> > 0 0
> > 0 1 2 3 4 5
> > 6 7
> > 0 0 0 0 0 0
> > 0 0
> > 8 9 : ; < =
> > > ?
> > 0 0 0 0 0 0
> > 0 0
> > @ A B C D E
> > F G
> > 0 1367337 0 2035064 0 1006372
> > 819775 0
> > H I J K L M
> > N O
> > 1926112 116395 1614109 356496 517731 1602298 578172
> > 2016009
> > P Q R S T U
> > V W
> > 3043393 401196 2152160 1553511 2149231 254108 3896558
> > 232409
> > X Y Z [ \\ ]
> > 706452 4353706 1292309 732210 0 45133419
> >> abcScore[34:62,c(1:4, 33:36)]
> > [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
> > A 1 1788 1798 8 260944 299631 332500 354717
> > B 0 0 0 0 0 0 0 0
> > C 59 36426 12705 227 39 131268 140461 154643
> > D 0 0 0 0 0 0 0 0
> > E 133 13391 4033 180 120135 0 0 0
> > F 0 0 0 0 0 259549 272238 287988
> > G 0 0 0 0 0 0 0 0
> > H 272 8815 2933 357 242741 234088 239917 239157
> > I 0 0 0 0 116395 0 0 0
> > J 472 4678 1656 586 0 195991 198114 187535
> > K 5 59 25 11 110052 84345 83744 77353
> > L 0 0 0 0 102438 144459 142049 128785
> > M 75 171 121 92 93945 119209 116492 103765
> > N 653 2992 1232 769 85940 91924 89084 79246
> > O 1227 2719 1754 1576 175403 187339 176997 160100
> > P 65834 57830 60855 65863 61910 0 0 0
> > Q 0 0 0 0 236290 0 0 0
> > R 3211 4482 4321 4215 0 0 0 0
> > S 0 0 0 0 68378 470434 426641 444948
> > T 4444 5373 5868 5637 0 0 0 0
> > U 0 0 0 0 0 0 0 0
> > V 9146 10270 11922 11547 543627 0 0 0
> > W 0 0 0 0 0 0 0 0
> > X 0 0 0 0 0 0 0 0
> > Y 28363 29080 34093 32873 0 0 0 0
> > Z 0 0 0 0 0 0 0 0
> > [ 0 0 0 0 0 0 0 0
> > \\ 0 0 0 0 0 0 0 0
> > ] 2104342 2040163 2074921 2094296 0 0 0 0
> >
> > Output from the last line shows how scores decrease from the first
> > four
> > cycles to the last four. Standard R and Biostrings commands can be
> > used
> > to ask many interesting questions, such as an overall quality score of
> > reads (e.g., summing the scores of individual nucleotides) and the
> > relationship between sequence characteristics (e.g., frequency of 'T')
> > and read quality.
> >
> >> atgc <- alphabetFrequency(sequences(seqq), baseOnly=TRUE)
> >> qscore <- alphabetFrequency(scores(seqq))
> >> dim(qscore)
> > [1] 2218237 256
> >> mb(object.size(qscore))
> > [1] 2166.2 MB
> >
> >> quality <- colSums(t(qscore) * 1:ncol(qscore))
> >> plot(density(quality)) # small secondary peak at low quality
> >> scores(seqq)[which(quality<2925)]
> > A BStringSet instance of length 87696
> > width seq
> > [1] 36 PPPPPPPPPPPPPEPPPPPOPPMOOPPOMMMPOOOJ
> > [2] 36 PPPPPPPPPPPPPPPPPPPPHPPPOPPMOPPPNKMA
> > [3] 36 PPPPPPPPPPPPPPPPPOPPPPPPEPMPPPMPOFAF
> > [4] 36 PPPPPPPPPPPPPPPPPPPPOPPPPPPPPPOPOOHK
> > [5] 36 PPPPPPPPPPPPPPPPPPPPPPPOPPOPOPPMKMLF
> > [6] 36 PPPPPPPPPPPPPPPPOPPPPPCPPPPHOOPPOOOO
> > [7] 36 PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPOOKO
> > [8] 36 PPPPPPPPPPPPPPPPPPPOPPPPPPPPJPOPOMAO
> > [9] 36 PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPOOOL
> > ... ... ...
> > [87688] 36 PPPPPPPPMPOPCOOCOOJEOMMCCOEEHCCMAAAA
> > [87689] 36 PTYOR]NTVVVJOHJCCPMEHCMHOJCECCCCAACA
> > [87690] 36 Y]]]NTVYYT]TPNRPNRYHCTECOHCHHCECHAAC
> > [87691] 36 YPRVNVYVPYHORCOCOPEPECJCCHCCJECHAAAA
> > [87692] 36 PPPPPPPPPPPPPHPPOCPOPPOPCHMMPPPEKHOO
> > [87693] 36 JJHPTTJYV]]]]JJCJCJJTCCOVCMHMCCOAAAA
> > [87694] 36 PPPPPPPPPPPPPPPPPOPPPJPPPCPOHHCMOCAF
> > [87695] 36 TYR]R]TN]YOEPTNERPPTVTCTVCCCRHOMIFAC
> > [87696] 36 YRRYJVVTPHVPPECHPNCMCCJCEECOCCCCAAAA
> >>
> >> t <- atgc[,"T"] / rowSums(atgc[,1:4])
> >> cor(t, quality)
> > [1] 0.154
> >
> > All of these operations are quick enough to perform in an interactive
> > session; the qscore is a large matrix (it can be made smaller by
> > choosing bounds that reflect allowable scores, e.g., 32:127), and its
> > transposition is relatively expensive.
> >
> > A final point to remember is that R stores a matrix m as a vector of
> > length nrow(m) * ncol(m). R has an internal limit on the size of a
> > vector equal to 2^32-1, so the maximum number of reads whose scores
> > can
> > be represented by alphabetFrequency is less than 2^32 / 256, or
> > about 16
> > million reads; this number of reads might be approached in a single
> > Solexa lane; a simple solution is to divide the DNAStringSet into
> > pieces
> > that are processed separately.
> >
> >
> >
> > I hope that the forgoing provides some indication of where we stand at
> > the moment. Again, it would be great to have feedback, and to hear of
> > other efforts. And again, the programming credit goes to Herve Pages.
> >
> > Martin
> >
> > The obligatory sessionInfo, plus some stats on processing this
> > document
> > (referencing, in the end, three different data sets).
> >
> >> sessionInfo()
> > R version 2.7.0 alpha (2008-03-28 r44972)
> > x86_64-unknown-linux-gnu
> >
> > locale:
> > LC_CTYPE
> > =
> > en_US
> > .UTF
> > -8
> > ;LC_NUMERIC
> > =
> > C
> > ;LC_TIME
> > =
> > en_US
> > .UTF
> > -8
> > ;LC_COLLATE
> > =
> > en_US
> > .UTF
> > -8
> > ;LC_MONETARY
> > =
> > C
> > ;LC_MESSAGES
> > =
> > en_US
> > .UTF
> > -8
> > ;LC_PAPER
> > =
> > en_US
> > .UTF
> > -8
> > ;LC_NAME
> > =
> > C
> > ;LC_ADDRESS
> > =C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] tools stats graphics grDevices utils datasets
> > [7] methods base
> >
> > other attached packages:
> > [1] BiostringsCinterfaceDemo_0.1.2
> > [2] BSgenome.Hsapiens.UCSC.hg18_1.3.2
> > [3] BSgenome_1.7.4
> > [4] Biostrings_2.7.41
> > [5] Biobase_1.99.4
> >> gc()
> > used (Mb) gc trigger (Mb) max used (Mb)
> > Ncells 9.37e+05 50.1 4.87e+06 260 1.93e+07 1033
> > Vcells 6.88e+08 5252.6 1.31e+09 10031 1.31e+09 9998
> >> proc.time()
> > user system elapsed
> > 356.7 16.6 378.5
> >
> > _______________________________________________
> > Bioc-sig-sequencing mailing list
> > Bioc-sig-sequencing at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
>
> _______________________________________________
> Bioc-sig-sequencing mailing list
> Bioc-sig-sequencing at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
>
More information about the Bioc-sig-sequencing
mailing list