[Bioc-sig-seq] Bioc short read directions
Loyal Goff
lgoff at rci.rutgers.edu
Wed Apr 2 15:27:18 CEST 2008
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?
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
More information about the Bioc-sig-sequencing
mailing list