Dear Bioconductor List, I am working through the easyRNASeq human sequence use case with a difference. I am using Bowtie/Tophat/Cufflinks aligned reads and transcripts. It works fine for getting counts of exons, transcripts, and genes for the genes.gtf file from the UCSD assembly. The difficulty comes in when I try to get counts for novel exons, genes and transcripts. To get this, I am uses the merged.gtf file which is generated by cuffmerge rather than the genes.gtf file from UCSD. I get Exon counts just fine (although I lose the old Ensembl Exon symbols). What I can't get is transcript counts. Here are excerpts from my session: > library(easyRNASeq) Loading required package: parallel Loading required package: genomeIntervals Loading required package: intervals Loading required package: BiocGenerics Attaching package: ‘BiocGenerics’ The following object(s) are masked from ‘package:stats’: xtabs The following object(s) are masked from ‘package:base’: anyDuplicated, cbind, colnames, duplicated, eval, Filter, Find, get, intersect, lapply, Map, mapply, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rbind, Reduce, rep.int, rownames, sapply, setdiff, table, tapply, union, unique Loading required package: Biobase Welcome to Bioconductor Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'. Loading required package: biomaRt Loading required package: edgeR Loading required package: limma Loading required package: Biostrings Loading required package: IRanges Attaching package: ‘IRanges’ The following object(s) are masked from ‘package:intervals’: expand, reduce Attaching package: ‘Biostrings’ The following object(s) are masked from ‘package:intervals’: type Loading required package: BSgenome Loading required package: GenomicRanges Loading required package: DESeq Loading required package: locfit locfit 1.5-8 2012-04-25 Attaching package: ‘locfit’ The following object(s) are masked from ‘package:GenomicRanges’: left, right Loading required package: lattice Attaching package: ‘DESeq’ The following object(s) are masked from ‘package:limma’: plotMA Loading required package: Rsamtools Loading required package: ShortRead Loading required package: latticeExtra Loading required package: RColorBrewer Warning messages: 1: replacing previous import ‘coerce’ when loading ‘intervals’ 2: replacing previous import ‘initialize’ when loading ‘intervals’ > library(BSgenome.Hsapiens.UCSC.hg19) > library(Rsamtools) > > chr.sizes=seqlengths(Hsapiens) > chr.sizes chr1 chr2 chr3 chr4 249250621 243199373 198022430 191154276 chr5 chr6 chr7 chr8 180915260 171115067 159138663 146364022 chr9 chr10 chr11 chr12 141213431 135534747 135006516 133851895 chr13 chr14 chr15 chr16 115169878 107349540 102531392 90354753 chr17 chr18 chr19 chr20 81195210 78077248 59128983 63025520 chr21 chr22 chrX chrY 48129895 51304566 155270560 59373566 chrM chr1_gl000191_random chr1_gl000192_random chr4_ctg9_hap1 16571 106433 547496 590426 chr4_gl000193_random chr4_gl000194_random chr6_apd_hap1 chr6_cox_hap2 189789 191469 4622290 4795371 chr6_dbb_hap3 chr6_mann_hap4 chr6_mcf_hap5 chr6_qbl_hap6 4610396 4683263 4833398 4611984 chr6_ssto_hap7 chr7_gl000195_random chr8_gl000196_random chr8_gl000197_random 4928567 182896 38914 37175 chr9_gl000198_random chr9_gl000199_random chr9_gl000200_random chr9_gl000201_random 90085 169874 187035 36148 chr11_gl000202_random chr17_ctg5_hap1 chr17_gl000203_random chr17_gl000204_random 40103 1680828 37498 81310 chr17_gl000205_random chr17_gl000206_random chr18_gl000207_random chr19_gl000208_random 174588 41001 4262 92689 chr19_gl000209_random chr21_gl000210_random chrUn_gl000211 chrUn_gl000212 159169 27682 166566 186858 chrUn_gl000213 chrUn_gl000214 chrUn_gl000215 chrUn_gl000216 164239 137718 172545 172294 chrUn_gl000217 chrUn_gl000218 chrUn_gl000219 chrUn_gl000220 172149 161147 179198 161802 chrUn_gl000221 chrUn_gl000222 chrUn_gl000223 chrUn_gl000224 155397 186861 180455 179693 chrUn_gl000225 chrUn_gl000226 chrUn_gl000227 chrUn_gl000228 211173 15008 128374 129120 chrUn_gl000229 chrUn_gl000230 chrUn_gl000231 chrUn_gl000232 19913 43691 27386 40652 chrUn_gl000233 chrUn_gl000234 chrUn_gl000235 chrUn_gl000236 45941 40531 34474 41934 chrUn_gl000237 chrUn_gl000238 chrUn_gl000239 chrUn_gl000240 45867 39939 33824 41933 chrUn_gl000241 chrUn_gl000242 chrUn_gl000243 chrUn_gl000244 42152 43523 43341 39929 chrUn_gl000245 chrUn_gl000246 chrUn_gl000247 chrUn_gl000248 36651 38154 36422 39786 chrUn_gl000249 38502 countrabtex <- easyRNASeq(filesDirectory=getwd(), + organism="Hsapiens", chr.sizes=chr.sizes, + readLength=58L, + annotationMethod="gtf", + annotationFile="merged.gtf", + count="exons", + gapped=TRUE, + filenames=bamfiles[1]) Checking arguments... Fetching annotations... Read 454815 records Summarizing counts... Processing SRR490224.bam Preparing output Warning messages: 1: In .Method(..., deparse.level = deparse.level) : number of columns of result is not a multiple of vector length (arg 1) 2: In easyRNASeq(filesDirectory = getwd(), organism = "Hsapiens", chr.sizes = chr.sizes, : There are 229182 features/exons defined in your annotation that overlap! This implies that some reads will be counted more than once! Is that really what you want? 3: In easyRNASeq(filesDirectory = getwd(), organism = "Hsapiens", chr.sizes = chr.sizes, : You enforce UCSC chromosome conventions, however the provided annotation is not compliant. Correcting it. 4: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : You enforce UCSC chromosome conventions, however the provided alignments are not compliant. Correcting it. 5: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : The read length stored in the object (probably provided as argument): 58 is not the same as the ones: 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 determined from the file: /Documents/clients/Phyllis/learning_easyRNAseq/SRR490224.bam Updating the readLength slot. 6: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : Not all the chromosome names in your chromosome size list 'chr.sizes' are present in your read file(s) (aln or bam). 7: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : The available chromosomes in both your read file(s) (aln or bam) and 'chr.sizes' list were restricted to their common term. These are: chr1, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr2, chr20, chr21, chr22, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chrM, chrX, chrY. > write.csv(countrabtex,"countrabtexmerged.csv" ) THE OUTPUT FILE CONTAINED AN EXON COUNT LIST. EXONS WERE LABELLED "XLOC". QUESTION: IS THERE A WAY TO GET THE ENSEMBL NUMBER FOR KNOWN EXONS AND RESERVE "XLOC" FOR NEW ONES? NOW FOR TRANSCRIPTS: > RnaSeq <- easyRNASeq(filesDirectory=getwd(), + organism="Hsapiens", chr.sizes=chr.sizes, + readLength=58L, + annotationMethod="gtf", + annotationFile="merged.gtf", + count="exons", + gapped=TRUE, + filenames=bamfiles[1], + outputFormat="RNAseq") Checking arguments... Fetching annotations... Read 454815 records Summarizing counts... Processing SRR490224.bam Preparing output Warning messages: 1: In .Method(..., deparse.level = deparse.level) : number of columns of result is not a multiple of vector length (arg 1) 2: In easyRNASeq(filesDirectory = getwd(), organism = "Hsapiens", chr.sizes = chr.sizes, : There are 229182 features/exons defined in your annotation that overlap! This implies that some reads will be counted more than once! Is that really what you want? 3: In easyRNASeq(filesDirectory = getwd(), organism = "Hsapiens", chr.sizes = chr.sizes, : You enforce UCSC chromosome conventions, however the provided annotation is not compliant. Correcting it. 4: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : You enforce UCSC chromosome conventions, however the provided alignments are not compliant. Correcting it. 5: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : The read length stored in the object (probably provided as argument): 58 is not the same as the ones: 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 determined from the file: /Documents/clients/Phyllis/learning_easyRNAseq/SRR490224.bam Updating the readLength slot. 6: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : Not all the chromosome names in your chromosome size list 'chr.sizes' are present in your read file(s) (aln or bam). 7: In fetchCoverage(rnaSeq, format = format, filename = filename, filter = filter, : The available chromosomes in both your read file(s) (aln or bam) and 'chr.sizes' list were restricted to their common term. These are: chr1, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr2, chr20, chr21, chr22, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chrM, chrX, chrY. > class(RnaSeq) [1] "RNAseq" attr(,"package") [1] "easyRNASeq" > > RnaSeqTranscripts<-transcriptCounts(RnaSeq) Error in aggregate.data.frame(as.data.frame(x), ...) : arguments must have same length > sessionInfo() R version 2.15.2 (2012-10-26) Platform: i386-apple-darwin9.8.0/i386 (32-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 easyRNASeq_1.4.2 [3] ShortRead_1.16.3 latticeExtra_0.6-24 [5] RColorBrewer_1.0-5 Rsamtools_1.10.2 [7] DESeq_1.10.1 lattice_0.20-10 [9] locfit_1.5-8 BSgenome_1.26.1 [11] GenomicRanges_1.10.5 Biostrings_2.26.2 [13] IRanges_1.16.4 edgeR_3.0.7 [15] limma_3.14.3 biomaRt_2.14.0 [17] Biobase_2.18.0 genomeIntervals_1.14.0 [19] BiocGenerics_0.4.0 intervals_0.13.3 loaded via a namespace (and not attached): [1] annotate_1.36.0 AnnotationDbi_1.20.3 bitops_1.0-4.2 DBI_0.2-5 [5] genefilter_1.40.0 geneplotter_1.36.0 grid_2.15.2 hwriter_1.3 [9] RCurl_1.95-3 RSQLite_0.11.2 splines_2.15.2 stats4_2.15.2 [13] survival_2.37-2 XML_3.95-0.1 xtable_1.7-0 zlibbioc_1.4.0 > QUESTION: IS THERE ANY WAY TO GET TRANSCRIPT COUNTS INCLUDING NEW TRANSCRIPTS ASSEMBLED BY CUFFMERGE? QUESTION: DO YOU HAVE ANY SUGGESTIONS FOR ANOTHER WAY OF GETTING COUNTS OF NOVEL TRANSCRIPTS ASSEMBLED BY CUFFMERGE? CUFFDIFF GIVES NON-INTEGRAL TRANSCRIPT NUMBERS SO I THINK THAT IT DOES SOME SORT OF NORMALIZATION WHICH I DON'T WANT. Thanks and best wishes, Rich Richard A. Friedman, PhD Associate Research Scientist, Biomedical Informatics Shared Resource Herbert Irving Comprehensive Cancer Center (HICCC) Lecturer, Department of Biomedical Informatics (DBMI) Educational Coordinator, Center for Computational Biology and Bioinformatics (C2B2)/ National Center for Multiscale Analysis of Genomic Networks (MAGNet)/ Columbia Initiative in Systems Biology Room 824 Irving Cancer Research Center Columbia University 1130 St. Nicholas Ave New York, NY 10032 (212)851-4765 (voice) friedman@cancercenter.columbia.edu http://cancercenter.columbia.edu/~friedman/ "Complex numbers! Ha! Ha! There is nothing weirder than imaginary numbers. Architects don't need to know complex numbers. Whenever I get a negative root for an area, I throw it out. And don't talk to me about quaternions. I am not going into computer animation." -Rose Friedman, age 16 [[alternative HTML version deleted]]