SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
What to include in my count table(s) for DESeq rndouglas Bioinformatics 5 11-08-2013 05:54 AM
Count table from BAM file with custom gtf tamari Bioinformatics 3 09-05-2013 02:42 PM
Alternative to cufflinks? - Bam to count table? hbt Bioinformatics 3 11-09-2012 06:09 AM
question on making BLAST db rdu Bioinformatics 4 01-12-2011 11:45 PM
question on making libraries from pcr products seqgirl123 Sample Prep / Library Generation 3 02-07-2010 07:55 AM

Reply
 
Thread Tools
Old 12-12-2013, 10:19 PM   #1
lucer105
Member
 
Location: Columbia, Missouri

Join Date: Nov 2013
Posts: 12
Question Question for making count table for edgeR

I am trying to make a count table for following edgeR..... but after running the edgeR, my marker gene didn't change as it suppose to change. (It suppose to change, confirmed by qPCR, and cuffdiff, and load .bam file in IGV, and it is huge change). I checked the count table content I made and the reads in the table is weird. I got 0, 0 reads for control and 1, 7 reads for treatment.... Anyway, I think I made the table wrong, would anyone check what I did? I check it over and over again and keep changing settings but still didn't hit the core of the problem.

a=ctrl, b=ctrl, a&b are biological replicates; c=treatment, d=treatment, c&d are biological replicates



library(GenomicFeatures)
library(GenomicRanges)
library(Rsamtools)
library(edgeR)
library(org.Mm.eg.db)

txdb=makeTranscriptDbFromUCSC(genome="mm9",tablename="refGene")
ex_by_gene=exonsBy(txdb,"gene")

####read .bam file in object
readsa=readBamGappedAlignments("/home/vagrant/genedata/a_accepted_hits.bam")
readsb=readBamGappedAlignments("/home/vagrant/genedata/b_accepted_hits.bam")
readsc=readBamGappedAlignments("/home/vagrant/genedata/c_accepted_hits.bam")
readsd=readBamGappedAlignments("/home/vagrant/genedata/d_accepted_hits.bam")
#####count reads
countsa = countOverlaps(ex_by_gene,readsa)
countsb = countOverlaps(ex_by_gene,readsb)
countsc = countOverlaps(ex_by_gene,readsc)
countsd = countOverlaps(ex_by_gene,readsd)

#### making count table
countTableabcd = data.frame(conditiona=countsa,conditionb=countsb,conditionc=countsc,conditiond=countsd,stringsAsFactors=FALSE)
rownames(countTableabcd)=names(ex_by_gene)
write.table(countTableabcd,file="countTableabcd.txt",sep="\t")
#### eliminate empty rows
x <- rowSums(countTableabcd==0)!=ncol(countTableabcd)
newCountTableabcd <- countTableabcd[x,]
write.table(newCountTableabcd,file="newcountTableabcd.txt",sep="\t")
lucer105 is offline   Reply With Quote
Old 12-13-2013, 12:37 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Just use featureCounts or htseq-counts to count the reads. I've never been a fan of the R methods, particularly since they usually involve reading the whole file into memory.
dpryan is offline   Reply With Quote
Old 12-13-2013, 02:34 AM   #3
bruce01
Senior Member
 
Location: .

Join Date: Mar 2011
Posts: 157
Default

You can use featureCounts in R, if the OP wants to do it all in that wonderful enviRonment=)
bruce01 is offline   Reply With Quote
Old 12-13-2013, 07:01 AM   #4
lucer105
Member
 
Location: Columbia, Missouri

Join Date: Nov 2013
Posts: 12
Default

Sounds like it is easy to fix.

But I just wondering, why did my method failed? At least for counting that single gene....
lucer105 is offline   Reply With Quote
Reply

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 08:42 PM.


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