SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics

Similar Threads
Thread Thread Starter Forum Replies Last Post
Perl script for annotating DEGseq results umnklang Bioinformatics 12 06-27-2012 07:34 AM
problem for velvet parameters and results hequn Bioinformatics 2 11-17-2011 05:56 AM
ERANGE mRNA counts problem yh253 Bioinformatics 0 08-06-2010 08:39 AM
some problem of solexa results biocc Illumina/Solexa 1 06-08-2010 02:22 PM

Reply
 
Thread Tools
Old 09-19-2012, 02:30 AM   #1
Jluis
Member
 
Location: Bilbao

Join Date: Apr 2012
Posts: 44
Default Problem annotating results derived from an table of counts

Dear all,

I am having problems annotating results from a table of counts obtained using the following R-script code:

library(Rsamtools)
library(GenomicFeatures)

bamlist=list()
src_files=list.files(pattern="*.bam")

for(filename in src_files)
{
tmp=readBamGappedAlignments(filename)
bamlist[[length(bamlist)+1]]=GRanges(seqnames=rname(tmp),
ranges=IRanges(start=start(tmp),end=end(tmp)),
strand=rep("*",length(tmp)))
}

names(bamlist)=src_files

#--------------------Fetching gene information
txdb=makeTranscriptDbFromUCSC(genome='mm9',tablename='knownGene')

#--------------------GRangesList object
tx_by_gene=transcriptsBy(txdb,'gene')

#--------------------include reads that overlap exons
ex_by_gene=exonsBy(txdb,'gene')

#------------------Finally, we get the number of reads that overlap with each gene

#Initialize table of counts
toc=data.frame(rep(NA,length(tx_by_gene)))
for(i in 1:length(bamlist))
{
toc[,i]=countOverlaps(tx_by_gene,bamlist[[i]])
}
#Fix up labels
rownames(toc)=names(tx_by_gene)
print(ncol(toc))
print(length(bamlist))
colnames(toc)=names(bamlist)

dim(toc)


head(toc)

write.table(toc, file="mm9_ToC.tsv",quote=FALSE,sep="\t")



The thing is that using this reference genome:

genome='mm9',tablename='knownGene'

The resulting table of counts shows up this IDs format:

IDs s1.bam s2.bam

100009600 121 83
100009609 0 3
100009614 0 0
100012 0 0
100017 121 79
100019 37 24


I can't explain myself which kind of identifier is this and how to convert it to other format that may allow me to infer the Gene Function related to the Gene Id...

May there be any problem with the code so that the Table of Counts isn't
generated correctly?

I'll appreciate any help and/or comments, because it's been a real headache trying to correlate such Ids to any G.O.
Jluis is offline   Reply With Quote
Old 09-19-2012, 02:39 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

The manual for GenomicFeatures says that those should be Entrez gene IDs, which is what I would have guessed from looking at them. That's usually more convenient for downstream processing than using actual gene names, which can overlap.
dpryan is offline   Reply With Quote
Old 09-19-2012, 03:00 AM   #3
Jluis
Member
 
Location: Bilbao

Join Date: Apr 2012
Posts: 44
Default

Dear dpryan,

Thank you very much for your explanation, I'll try to convert those Ids to ensembl and ncbi formats to help bench scientist to correlate them with GO functions.

Best wishes
Jluis is offline   Reply With Quote
Reply

Tags
annotation, gene onthology, table of counts

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:36 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO