gmane.comp.lang.r.sequencing

bioc-sig-sequencing list will be removed

Dan Tenenbaum — 2011-09-30T19:51:11

Hi all,

We will be closing the bioc-sig-sequencing email list. The reasons are
because sequencing is now a mainstream part of the Bioconductor
project, posts to one list are often relevant to both, new users are
confused by the appropriate list for posting, and a large fraction of
bioc-sig-sequencing subscribers also subscribe to the main
Bioconductor mailing list.

The list archive will continue to be available at

 https://stat.ethz.ch/pipermail/bioc-sig-sequencing/

For continued support, you may use the Bioconductor mailing list

 http://bioconductor.org/help/mailing-list/

including a new feature allowing questions posted as a guest

 http://bioconductor.org/help/mailing-list/mailform/


The bioc-sig-sequencing list will be removed on Tuesday, October 4th.

Posting to the Bioconductor list without subscribing

Dan Tenenbaum — 2011-09-30T19:12:44

Hi all,

We've added a new feature to the Bioconductor web site: a form through which
you can post to the Bioconductor email list without first subscribing to it.
The form contains a captcha to protect against spammers.

The form is here:
http://bioconductor.org/help/mailing-list/mailform/

and it can also be reached through the Mailing Lists link on the main page
at http://bioconductor.org.

Enjoy,
Dan

[[alternative HTML version deleted]]

Re: Another ScanBamParam suggestion

Ivan Gregoretti — 2011-09-30T17:14:10

Hi Martin,

Thanks for flag=scanBamFlag(isValidVendorRead=TRUE). I didn't know it
existed.

Regarding parallelisation of I/O, I completely understand the challenge. I
can only add that in my experience with R and big sequencing files, the
bottleneck has been invariably the CPU and not the disk. I work in a 16 core
Intel Xeon X5570 2.93GHz with 144RAM. The disk is accessible through a
gigabit network.

I hope we get to hear some ideas from the community.

Thank you,

Ivan



Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1016 and 1-301-496-1592
Fax: 1-301-496-9878


On Fri, Sep 30, 2011 at 12:11 PM, Martin Morgan <mtmorgan-q9hIisBwmLrYtjvyW6yDsg< at >public.gmane.org> wrote:


[[alternative HTML version deleted]]

Re: Another ScanBamParam suggestion

Martin Morgan — 2011-09-30T16:11:12

Hi Ivan --

in principle, flag=scanBamFlag(isValidVendorRead=TRUE) will do this; it 
requires that the flag is set in the BAM file.


I'm actually revising the I/O a little at the moment; I'll implement a 
better strategy for reading in the data. When I look at 'top' on my 
system, I see the CPU running at say 50% which implies disk input is the 
bottleneck; probably this is on our system administration end, where the 
large storage required for BAM files doesn't have completely adequate 
performance. This I/O is tricky to guage, because the next time through 
the BAM file input _is_ CPU limited and much faster -- the disk system 
has done some clever buffering. But in real use cases I wouldn't see the 
benefit of that buffering since I wouldn't be revisiting the file.

In terms of parallel throughput it might often be appropriate to 
parallelize at a higher level, e.g., iterating over regions of interest 
(e.g., GRanges defining chromosomes, with the iteration via lapply) and 
a function FUN tasked with in

Another ScanBamParam suggestion

Ivan Gregoretti — 2011-09-30T14:48:32

Following Janet's example, I would also like to propose an upgrade to
ScanBamParam:

It would be great if we could tell ScanBamPram that we want to load
only the reads that passed the vendor's quality filter.

In other words, the functionality I am suggesting is analogous to the
filter in readAligned() from the ShortRead library.


With the new release of Illumina sequencing reagents (version 3) you
get 200 million reads per lane from the HiSeq 2000. In my view, with
samples that big becoming popular, any investment in "read in"
efficiency is a good investment. I would be happy to provide a sample
BAM for those interested in addressing this suggestion.

It is also my humble opinion that we should start considering
parallelisation for reading in. I hope that I am not just wishing too
much.

Thank you,

Ivan

Re: ScanBamParam suggestion

Martin Morgan — 2011-09-30T00:55:59

yep, should be in 1.5.68 when it becomes available. Thanks for the 
suggestion. Martin

ScanBamParam suggestion

Janet Young — 2011-09-30T00:26:24

Hi,  

I have a suggestion to make ScanBamParam easier to use for coding amateurs like myself (I'm still sometimes confused with the many ways to encode genomic regions):

Is it easy/possible to change bamWhich function to accept GRanges objects, rather than requiring RangesList?  See below...

thanks, as usual,

Janet



############
library(Rsamtools)
myGR <- GRanges(seqnames="chr1",ranges=IRanges(start=1,end=100))

# I can use GRanges as the "which" argument if I do it when I create the ScanBamParam object
myparams1 <- ScanBamParam(which=myGR)

#but not if I try to set later
myparams2 <- ScanBamParam()
bamWhich(myparams2) <- myGR
### Error in checkSlotAssignment(object, name, value) : 
###   assignment of an object of class "GRanges" is not valid for slot "which" 
### in an object of class "ScanBamParam"; is(value, "RangesList") is not TRUE

## it's OK, though coercion does work.  
bamWhich(myparams2) <- as(myGR,"RangesList")

sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: i386-apple-darwin9.8.0/i3

Re: as.data.frame on GRanges object with DNAStringSet in values

Michael Lawrence — 2011-09-29T21:17:30

I saw that all coercions to atomic vectors from AtomicList are now
deprecated. You had proposed deprecating as.vector(), because it should not
unlist, and I agreed. Really as.vector() should return an ordinary R list.
However, as.character(), as.numeric(), etc, in base R will unlist. I'd like
to keep consistency with base R. Do we really need to deprecate those, as
well?

Michael

2011/6/15 Michael Lawrence <michafla-RuTDbSqP/YI< at >public.gmane.org>


[[alternative HTML version deleted]]

_______________________________________________
Bioc-sig-sequencing mailing list
Bioc-sig-sequencing-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

Re: Reads in 3'utr

rohan bareja — 2011-09-29T02:46:36

Hi,
I have summed the counts now to the gene level for 3'UTR.
I want to assess the relative amount of each 3âUTR end usage such as what percentage of reads comes from each 3âUTR isoform?Â Â I want to identify the different 3âUTR ends for each gene to get alternative 3'UTR Â usage(disease vs control)?
Do you have any idea about how to proceed?
Thanks,Rohan


--- On Sat, 24/9/11, Valerie Obenchain <vobencha-q9hIisBwmLrYtjvyW6yDsg< at >public.gmane.org> wrote:

From: Valerie Obenchain <vobencha-q9hIisBwmLrYtjvyW6yDsg< at >public.gmane.org>
Subject: Re: [Bioc-sig-seq] Reads in 3'utr
To: "rohan bareja" <rohan_1925-/E1597aS9LQxFYw1CcD5bw< at >public.gmane.org>
Cc: bioc-sig-sequencing-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Date: Saturday, 24 September, 2011, 4:24 AM



  

    
  
  
    On 09/23/2011 02:57 PM, rohan bareja wrote:
    
      
        
          
            Hi,
              

              
              
                utr=threeUTRsByTranscript(txdb,use.names=FALSE)
                So,utr

Re: Reads in 3'utr

Valerie Obenchain — 2011-09-29T02:56:34

DESeq and edgeR vignettes.

Valerie


On 09/28/11 19:46, rohan bareja wrote:

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

Re: a wired problem

Harris A. Jaffee — 2011-09-28T21:45:19

I forgot to say that in this code,



you do NOT want to select from 'reads' using 'result' as is, since
its values are all 0 or 1:

 > str(result)
  int [1:25000] 0 0 0 0 0 0 0 0 0 0 ...

 > table(result)
result
     0     1
24959    41

I suppose you are getting 41 instances of your first read, namely

   40-letter "DNAString" instance
seq: GTTTTCTCATTNGAAATTTTGTTACCGCAAAANNNCCACC

when instead, you want something like

reads[result == 1]

On Sep 28, 2011, at 5:27 PM, Harris A. Jaffee wrote:

Re: a wired problem

Harris A. Jaffee — 2011-09-28T21:27:15

I don't understand where your 'subject1' sequence:

came from.

Modulo that, everything looks fine.

In my hands, your vcountPattern call finds 41 40-letter elements
of 'seqs', all equal to GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA,
that match your 41-letter pattern, subject to your

max.mismatch=1, with.indels=TRUE'

They are all equal to your pattern with its first letter deleted.

I haven't looked closely, but I don't see a reason to doubt this
result.

The following calls would seem to be the sanity checks you want.
They at least confirm the 41 fuzzy matches.

 > countPattern(PCR2rc, "GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA",
max.mismatch=1, with.indels=TRUE)
[1] 1

 > matchPattern(PCR2rc, "GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA",
max.mismatch=1, with.indels=TRUE)

Views on a 40-letter BString subject
subject: GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA
views:
     start end width
[1]     1  40    40 [GATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA]

On Sep 28, 2011, at 4:14 PM, wang peter wrote:

a wired problem

wang peter — 2011-09-28T20:14:24

dear all:
     Using vcountPattern, i found some matched sequences.
but those are not similar to the pattern.
     see such coding

rm(list=ls())
reads <- readFastq(fastqfile);#downloaded from
http://biocluster.ucr.edu/~tbackman/query.fastq
seqs <- sread(reads);
PCR2rc<-DNAString("AGATCGGAAGAGCTCGTATGCCGTCTTCTGCTTGAAACAAA")
result <- vcountPattern(PCR2rc, seqs, max.mismatch=1, min.mismatch=0,
with.indels=TRUE, algorithm="indels")
reads <- reads[result]
seqs <- sread(reads)
sum(result)
     then using countPattern, i found they are really not match

subject1 = "GTTGGTGCAAACATTAGTTCTTCTGTTGGTGCAACCTTTG"
result <- countPattern(PCR2rc, subject1, max.mismatch=1, min.mismatch=0,
with.indels=TRUE)
[1] 0

shan gao

[[alternative HTML version deleted]]

coverage vector longer than covered area (ScanBamand IRanges)?

Francesco Lescai — 2011-09-27T19:53:45

Hi there,
may be I am missing something obvious in my code.
I extracted information from a bam file on 11 regions.
when I calculate the coverage on one of these regions, I get a vector which is much longer than the number of bases of the region itself.

I tried several times to replicate this problem on a small IRanges object but - of course - everything was ok.
Anyone has an idea?

that's what I've done:

   chr     start       end
1    1 195915909 197867662
2    2   1199920   2234892 [...]

scanned my bam within those regions
which=GRanges(regions$chr,IRanges(regions$start,regions$end))
ggf<-scanBam("s_7_1_sequence.txt.novo.rmdup.bam_sorted.bam", 
             param=ScanBamParam(which=which,
                                what = c("pos", "qwidth"),
                                flag =scanBamFlag(isUnmappedQuery =FALSE)))

created a summary function, as suggested in the IRanges vignette
summaryFunction<-function(seqname,bamfile,...){
  x<-bamfile[[seqname]]
  coverage(IRanges(x[["pos"]],width=x[["qwidth"

Re: question about trimLRPatterns

Harris A. Jaffee — 2011-09-27T21:24:31

Let's hope this is the bug Herve mentioned yesterday, whose fix  
should appear soon.


 > showMethods(which.isMatchingStartingAt, includeDefs=TRUE)
then
 > Biostrings:::.matchPatternAt

rev is in base (see ?rev)  [different from *reverse*, from IRanges]

 > Biostrings:::normargPattern


Again, I can't tell what you're asking here.

Re: what is

Harris A. Jaffee — 2011-09-27T20:10:10

matchPattern does not support your non-integral setting  
max.mismatch=0.2.
The help page is very clear about that with the phrase "number of ...":

max.mismatch, min.mismatch: The maximum and minimum number of
           mismatching letters allowed

Your 0.2 is converted to 0L by as.integer(max.mismatch), as here:

 > Biostrings:::normargMaxMismatch(0.2)
[1] 0

I think the following call gives what you want:

 > matchPattern(pattern1, subject1, max.mismatch=1, with.indels=TRUE)
   Views on a 81-letter BString subject
subject:  
ATCGAGATCGGAAGAGCGGTTCAGCAGGAATGCC...TATGCCGTCTTCTGCTTGAAAAAAAAAAAATATT
views:
     start end width
[1]     1  65    65  
[ATCGAGATCGGAAGAGCGGTTCAGCAGG...CCGATCTCGTATGCCGTCTTCTGCTTG]

However, the trimLRPatterns call you wrote is suspicious because you  
don't
mention any 'subject' (only *subject1*):


But, if we use subject1 instead there, we get

 > trimLRPatterns(Lpattern = pattern1, subject = subject1,  
max.Lmismatch=0,
with.Lindels=TRUE)
[1] "AAAAAAAAAAAATATT"

And this happe

what is

wang peter — 2011-09-27T19:23:57

hello every one,

i have such coding

"ATCGAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAATATT"

"AATCGAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTG"
with.indels=TRUE)
  Views on a 81-letter BString subject
subject:
ATCGAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAATATT
views: NONE
max.Lmismatch=0,with.Lindels=TRUE)
[1] "AAAAAAAAAAAATATT"

i already allow the indels,but why matchPattern cannot find the pattern in
subject
what does with.indels mean?
i am confused
thx
shan gao

[[alternative HTML version deleted]]

Re: Reads in 3'utr

rohan bareja — 2011-09-26T16:16:33

Hi Valerie,
Thanks a lot..It worked finally.. So now I have a data frame for the geneIds ,TranscriptIds and the counts (3'utr) which is given below:
  GENE     TX     countsUTRctl[1,] "148398" "1121" "2"         [2,] "339451" "1118" "0"         [3,] "84069"  "1116" "0"         [4,] "84069"  "1119" "11"        [5,] "9636"   "1126" "11"        [6,] "375790" "1127" "0"     
Now I want to do differential expression of genes using DESeq,so do I have to merge the two same genes and its counts such as geneID 84069 (from above ) or i can proceed with the above dataframe?If I have to merge them how do I do that?
Thanks,Rohan
--- On Sat, 24/9/11, Valerie Obenchain <vobencha-q9hIisBwmLrYtjvyW6yDsg< at >public.gmane.org> wrote:

From: Valerie Obenchain <vobencha-q9hIisBwmLrYtjvyW6yDsg< at >public.gmane.org>
Subject: Re: [Bioc-sig-seq] Reads in 3'utr
To: "rohan bareja" <rohan_1925-/E1597aS9LQxFYw1CcD5bw< at >public.gmane.org>
Cc: bioc-sig-sequencing-0bNBQ1PAWB4BXFe83j6qeQ< at >public.gmane.org
Date: Saturda

Re: remove adaptor before remove barcode

Harris A. Jaffee — 2011-09-27T16:26:07

I'm sorry, but I can't guess what "my method" was, nor what you are  
asking now.
Also, whether you are saying that something works at one end of a  
'subject' but
not the other (so you had to reflect?), and what 'trimmedReads' is.

On Sep 26, 2011, at 10:39 PM, wang peter wrote:

Re: Read big fastq files in chunks

Anita Lerch — 2011-09-27T14:17:21

Great. I think that is the right way.
I know the reference classes, but I don't have any working experience
with them. I start with implementing and keep you informed.

Anita

Re: Read big fastq files in chunks

Martin Morgan — 2011-09-27T12:33:19

Hi Anita -- this needs to be added to ShortRead, and I can do so later 
this week. But if you'd like to help, I would like to create a reference 
class

   .FastqStreamer <- setRefClass("FastqStreamer", <...>)

a constructor

   FastqStreamer <- function(con, n=1e6..., verbose=FALSE)

and a method

   yield,FastqStreamer-method

so that streaming through a file looks like

   f <- FastqStreamer(con, n=1e6, verbose=FALSE)
   while (length(srq <- yield(f))) {
        ## work
   }

I think .FastqStreamer would share some infrastructure with 
.FastqSampler, and would take the same approach as in 
R/methods-Sampler.R -- treat the file at the R level as a connection 
that is input with readBin() so that R's connection interface can be 
used, do some R-level buffering, and translate the raw() to ShortReadQ 
objects in C (like src/sampler.c:sampler_rec_parser). Probably there is 
some code re-factoring and a little reference class hierarchy .FastqFile 
as parent of .FastqStreamer and .FastqSampler.

Martin