Thanks Martin, I will try this.

On Thu, Nov 4, 2010 at 10:06 AM, Martin Morgan <mtmorgan@fhcrc.org> wrote:

> On 11/04/2010 09:45 AM, Changbin Du wrote:
> > Thanks, Jim!
> >
> > This is not what I want,  What I want is calculate the percentage of
> reads
> > bigger or equal to that reads in each position.MY output is like the
> > following:
>
> Hi Changbin -- I might be repeating myself, but the Bioconductor
> packages IRanges and GenomicRanges are designed to work with this sort
> of data, and include 'coverage' functions that do what you're interested
> in. Look into ?GRanges if interested.
>
>
>
> http://bioconductor.org/help/bioc-views/release/BiocViews.html#___HighThroughputSequencing
>
> Martin
>
> > for row 1, all the reads is >= 4, so the cover_per is 100,
> > for row 2, 99 % reads >=4, so the cover_per is 99.
> >> head(final)
> >   cover_per reads
> > 1       100     4
> > 2        99     8
> > 3        98    13
> > 4        97    14
> > 5        96    17
> > 6        95    20
> >
> > I attached the input file with this email. This file is only 100 rows,
> very
> > small. MY original data set is 3384766 rows.
> >
> >  >
> >
> matt<-read.table("/house/groupdirs/genetic_analysis/mjblow/ILLUMINA_ONLY_MICROBIAL_GENOME_ASSEMBLY/4083340/STANDARD_LIBRARY/GWZW.994.5.1129.trim_69.fastq.19621832.sub.sorted.bam.clone.depth",
> > sep="\t", skip=0, header=F,fill=T) #
> >> dim(matt)
> > [1] 3384766       2
> >
> > Thanks so much for your time!
> >
> >> matt<-read.table("/home/cdu/operon/dimer5_0623/matt_test.txt", sep="\t",
> > skip=0, header=F,fill=T) #
> >> names(matt)<-c("id","reads")
> >> dim(matt)
> > [1] 100   2
> >> cover<-matt$reads
> >>
> >>
> >> #calculate the cumulative coverage.
> >> cover_per<-function (data) {
> > + output<-numeric(0)
> > + for (i in data) {
> > +           x<-(100*sum(ifelse(data >= i, 1, 0))/length(data))
> > +           output<-c(output, x)
> > +                 }
> > + return(output)
> > + }
> >>
> >>
> >> result<-cover_per(cover)
> >> head(result)
> > [1] 100  99  98  97  96  95
> >>
> >> final<-data.frame(result, cover)
> >>
> >> names(final)<-c("cover_per", "reads")
> >> head(final)
> >   cover_per reads
> > 1       100     4
> > 2        99     8
> > 3        98    13
> > 4        97    14
> > 5        96    17
> > 6        95    20
> >
> >
> >
> >
> >
> > On Thu, Nov 4, 2010 at 9:18 AM, jim holtman <jholtman@gmail.com> wrote:
> >
> >> Is this what you want:
> >>
> >>> x
> >>            id reads
> >> 1   Contig79:1     4
> >> 2   Contig79:2     8
> >> 3   Contig79:3    13
> >> 4   Contig79:4    14
> >> 5   Contig79:5    17
> >> 6   Contig79:6    20
> >> 7   Contig79:7    25
> >> 8   Contig79:8    27
> >> 9   Contig79:9    32
> >> 10 Contig79:10    33
> >> 11 Contig79:11    34
> >>> x$percent <- x$reads / max(x$reads) * 100
> >>> x
> >>            id reads   percent
> >> 1   Contig79:1     4  11.76471
> >> 2   Contig79:2     8  23.52941
> >> 3   Contig79:3    13  38.23529
> >> 4   Contig79:4    14  41.17647
> >> 5   Contig79:5    17  50.00000
> >> 6   Contig79:6    20  58.82353
> >> 7   Contig79:7    25  73.52941
> >> 8   Contig79:8    27  79.41176
> >> 9   Contig79:9    32  94.11765
> >> 10 Contig79:10    33  97.05882
> >> 11 Contig79:11    34 100.00000
> >>
> >>
> >> On Thu, Nov 4, 2010 at 11:46 AM, Changbin Du <changbind@gmail.com>
> wrote:
> >>> HI, Dear R community,
> >>>
> >>> I have one data set like this,  What I want to do is to calculate the
> >>> cumulative coverage. The following codes works for small data set
> (#rows
> >> =
> >>> 100), but when feed the whole data set,  it still running after 24
> hours.
> >>> Can someone give some suggestions for long vector?
> >>>
> >>> id                reads
> >>> Contig79:1    4
> >>> Contig79:2    8
> >>> Contig79:3    13
> >>> Contig79:4    14
> >>> Contig79:5    17
> >>> Contig79:6    20
> >>> Contig79:7    25
> >>> Contig79:8    27
> >>> Contig79:9    32
> >>> Contig79:10    33
> >>> Contig79:11    34
> >>>
> >>>
> >>
> matt<-read.table("/house/groupdirs/genetic_analysis/mjblow/ILLUMINA_ONLY_MICROBIAL_GENOME_ASSEMBLY/4083340/STANDARD_LIBRARY/GWZW.994.5.1129.trim_69.fastq.19621832.sub.sorted.bam.clone.depth",
> >>> sep="\t", skip=0, header=F,fill=T) #
> >>> dim(matt)
> >>> [1] 3384766       2
> >>>
> >>> matt_plot<-function(matt, outputfile) {
> >>> names(matt)<-c("id","reads")
> >>>
> >>>  cover<-matt$reads
> >>>
> >>>
> >>> #calculate the cumulative coverage.
> >>> + cover_per<-function (data) {
> >>> + output<-numeric(0)
> >>> + for (i in data) {
> >>> +           x<-(100*sum(ifelse(data >= i, 1, 0))/length(data))
> >>> +           output<-c(output, x)
> >>> +                 }
> >>> + return(output)
> >>> + }
> >>>
> >>>
> >>>  result<-cover_per(cover)
> >>>
> >>>
> >>> Thanks so much!
> >>>
> >>>
> >>> --
> >>> Sincerely,
> >>> Changbin
> >>> --
> >>>
> >>>        [[alternative HTML version deleted]]
> >>>
> >>> ______________________________________________
> >>> R-help@r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible code.
> >>>
> >>
> >>
> >>
> >> --
> >> Jim Holtman
> >> Cincinnati, OH
> >> +1 513 646 9390
> >>
> >> What is the problem that you are trying to solve?
> >>
> >
> >
> >
> >
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
>
> --
> Computational Biology
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109
>
> Location: M1-B861
> Telephone: 206 667-2793
>



-- 
Sincerely,
Changbin
--

Changbin Du
DOE Joint Genome Institute
Bldg 400 Rm 457
2800 Mitchell Dr
Walnut Creet, CA 94598
Phone: 925-927-2856

	[[alternative HTML version deleted]]

