There is no such function in Rsubread currently, although we might add it in the future. But what I do for looking at duplicated reads was writing a simple C program to retrieve read sequences from the FASTQ file (you can do this using R as well and the running time will not be much different because the major computational cost here is input and output), and then run unix commands 'sort' and 'uniq -c' to get the number of occurrences of each read sequence. Similarly, you can write a simple program to retrieve the mapping locations (chr id + coordinate) of each read from the alignment results (SAM file), save them to a file and run 'sort' and 'uniq -c' commands to get the number of mapped reads for each chromosomal location. But you should be aware that not all duplicate alignments could come from PCR contamination, because RNA-seq data have much higher coverage than gDNA-seq data (with the same library sizes) since the transcriptome size is only about 2% of the genome size. An RNA-seq dataset which includes 10M 100 bp SE reads has a ~17X coverage. So you will expect to see quite a few duplicate alignments in your RNA-seq data. Hope this help. Wei On Oct 13, 2011, at 7:37 PM, Tim Triche, Jr. wrote: > Cool! I started using featureCounts() by 'gene' and by 'exon' almost immediately. However, when we run the same RNA-seq data through our local Bowtie/Tophat pipeline, we find that one of the samples has extensive PCR contamination (lots and lots of duplicate alignments). Is there a simple way to detect the same in Rsubread? The latter is MUCH faster... > > Thanks again for a great package and a very fast aligner. > > --t > > On Wed, Oct 12, 2011 at 10:39 PM, Wei Shi <shi@wehi.edu.au> wrote: > Hi Tim, > > The featureCounts() function in Rsubread counts the number of reads mapped to genes or exons and returns a list which includes a read count matrix and annotation information. The returned list object can be easily transformed to a DGEList object which can then be fed into downstream expression analysis packages like edgeR or limma. The Exact tests or GLM likelihood ratio tests in edgeR, or voom() function in limma, can then be used to perform differential expression analysis for the RNA-seq data. > > We have established a pipeline to perform complete analysis for RNA- seq data using R packages including Rsubread, edgeR and limma. This pipeline starts with the FASTQ/FASTA files, which is the output of next-generation sequencers, through to read alignment, read count summarization for genes or exons, and to biological variation estimation and statistical testing of differential expression, and finally yields lists of differentially expressed genes. The repertoire of functions included in limma and edgeR can be readily used for further downstream analysis such as gene set tests and so on. > > This is a quick and easy-to-use pipeline and it is now being used in our routine RNA-seq data analysis. Hope you will find it useful. > > Cheers, > Wei > > On Oct 11, 2011, at 9:32 AM, Tim Triche, Jr. wrote: >> but as you pointed out, a very nice feature of Rsubread for me is featureCounts. It's also really fast, so if I can get away with using it, I am thinking that I will do so from now on... >> >> Thanks for a fast and easy to use alignment package, >> >> --Tim >> >> On Mon, Oct 10, 2011 at 3:26 PM, Wei Shi <shi@wehi.edu.au> wrote: >> Dear Tim, >> >> Yes, Subread performs read mapping by mapping a set of 16 mers extracted >> from each read to the genome, counting the number of mapped 16 mers at >> each candidate location and then choosing the one which is mapped to by >> the majority of the 16 mers as the mapping location of the read. The >> fundamental difference between Subread and other aligners is that it uses >> a voting method rather than a extension method to determine the mapping >> locations of the reads, which makes it a lot faster and more sensitive. >> >> Subread has a both a C version and an R version. The C version is freely >> available from sourceforge and the R version is included in the Rsubread >> package in Bioc. >> >> The Rsubread package also includes a function called featureCounts, which >> can be used to count the number of reads for each exon or gene. So this >> function will be useful for you to look at the differential expression at >> both gene level and exon level. >> >> Another function which might be useful for your data analysis is the >> subjunc funtion, which is designed to discover exon junctions. Subjunc >> uses an idea similar to that of Subread. Our preliminary results showed >> that subjunc outperformed competing junction detectors in terms of speed, >> sensitivity and accuracy. >> >> The devel version of Rsubread package includes a lot of our recent >> development for both Subread aligner and Subjunc junction detector, so I >> would recommend using the devel version if you want to try the Rsubread >> package. >> >> >> Hope this helps. >> >> Cheers, >> Wei >> >> >> >> >> >> > A professor sent me a bunch of raw RNA-seq reads (as FASTQ files) and I >> > want >> > to align them, and I couldn't really make heads or tails of the options, >> > so >> > I listened to what Phil Green told me at a conference and looked around >> > for >> > a sensible word-nucleated aligner like he described. It seems that >> > Rsubread >> > works this way? >> > >> > http://sourceforge.net/projects/subread/ >> > >> > I would like BAM files as intermediate output, but my real interest is >> > differential exon usage in differentiating cells. Given that the reads I >> > have to align are relatively short (36bp, SE), is there an advantage or >> > disadvantage in using subread compared to other options? And when I'm >> > done >> > trimming and aligning, I could choose raw counts, conditional quantile >> > normalized counts, or something like RPKM to summarize how often a given >> > exon seems to have been transcribed. I read this: >> > >> > http://seqanswers.com/forums/showthread.php?t=586 >> > >> > and I see that packages using a Gamma prior for the dispersion of a >> > Poisson >> > count model benefit from having raw counts. If I am after correlated >> > changes in exon usage depending on other sequence features, is it >> > reasonable >> > to use (say) 'cqn' on the raw counts, then log-transform and work with >> > those >> > normalized counts? >> > >> > Thanks for any suggestions, >> > >> > -- >> > Tim Triche, Jr. >> > USC Biostatistics >> > >> > [[alternative HTML version deleted]] >> > >> > _______________________________________________ >> > Bioconductor mailing list >> > Bioconductor@r-project.org >> > https://stat.ethz.ch/mailman/listinfo/bioconductor >> > Search the archives: >> > http://news.gmane.org/gmane.science.biology.informatics.conductor >> > >> >> >> >> ______________________________________________________________________ >> The information in this email is confidential and intended solely for the addressee. >> You must not disclose, forward, print or use it without the permission of the sender. >> ______________________________________________________________________ >> >> >> >> -- >> If people do not believe that mathematics is simple, >> it is only because they do not realize how complicated life is. >> John von Neumann >> >> > > > ______________________________________________________________________ > The information in this email is confidential and inte...{{dropped:25}}