Hi! 

I'm doing a differential gene expression analysis on two time points for a human iPS cell line, using 100bp PE RNA-seq data mapped with tophat and aggregated to generate a count table using htseq-count in intersection-strict mode. While  the differentiation process changes the phenotype of the cell line significantly, I'm still not sure I would expect >15% of genes in the genome to change significantly. 

I'm especially worried about the fact that after I filter for lowly expressed genes, ~50% of those for which testing is carried out are called significant. Naturally, I can filter based on cpm/logFC/detection by both tools before coming to a "significant" DG table, but I'm wondering if there is something intrinsically wrong with my data/samples, and if so - what? Or am I just making a very silly mistake in my code somewhere?

Thanks in advance for your help,
Darya Vanichkina



edgeR BCV plot (Disp = 0.0189 , BCV = 0.1375 )

http://i.imgur.com/Pf32qeBh.png

edgeR smearplot

http://i.imgur.com/bs6JyoEh.png

DESeq MAplot
http://i.imgur.com/Ds9nfwth.png

## Sample code used
# edgeR -----
row.names(counttable.sam) <- counttable.sam$Gene
counttable.sam$Gene <- NULL

y <- DGEList(counts=counttable.sam[1:6], group=rep(1:2,each=3), genes=row.names(counttable.sam))

keep <- rowSums(cpm(y)>1) >= 1
y <- y[keep,]
dim(y)
# 16002    6
# Kept, out of 62890 genes in reference (gencode 17). Reads were counted to gene features using htseq-count in intersection-strict mode

y$samples$lib.size <- colSums(y$counts)
y <- calcNormFactors(y)
plotMDS(y)
y <- estimateCommonDisp(y, verbose=TRUE)
# Disp = 0.0189 , BCV = 0.1375 
y <- estimateTagwiseDisp(y)
plotBCV(y)
et <- exactTest(y)
summary(de <- decideTestsDGE(et))
# Out of 16002 tags for which the counts are not too low, 6263 are upregulated and 6235 are down-regulated
detags <- rownames(y)[as.logical(de)]
plotSmear(et, de.tags=detags)
abline(h=c(-1, 1), col="blue")


# DESeq -----
# Filter the lowly expressed tags
rs <- rowSums (counttable.sam)
use <- (rs > 10)
table(use)
# use
# FALSE  TRUE 
# 41704 21186 
counttable.sam <- counttable.sam[ use, ]
conds <- factor(c("day0", "day0", "day0", "day34", "day34", "day34"))

d <- newCountDataSet(counttable.sam, conds)
d <- estimateSizeFactors(d)
sizeFactors(d)
#  d0_1   d0_2   d0_3  d34_1  d34_2  d34_3 
# 2.0770099 1.4882383 0.7514035 0.4643261 1.3976540 0.6846997
d <- estimateDispersions(d, sharingMode="maximum")
plotDispEsts(d) 
DESeqRes <- nbinomTest(d, "day0", "day34")
DESeqRes$GeneShort <- substr(DESeqRes$id, 1, 15)
DESeqResSig <- subset(DESeqRes, padj < 0.01)
# 10163/21186 genes were called significant by DESeq ...
plotMA(DESeqRes)



> sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq_1.12.0         lattice_0.20-15      locfit_1.5-9.1       Biobase_2.20.1      
 [5] BiocGenerics_0.6.0   ggplot2_0.9.3.1      BiocInstaller_1.10.2 edgeR_3.2.4         
 [9] limma_3.16.6         biomaRt_2.16.0      

loaded via a namespace (and not attached):
 [1] annotate_1.38.0      AnnotationDbi_1.22.6 colorspace_1.2-2     DBI_0.2-7           
 [5] dichromat_2.0-0      digest_0.6.3         genefilter_1.42.0    geneplotter_1.38.0  
 [9] grid_3.0.1           gtable_0.1.2         IRanges_1.18.2       labeling_0.2        
[13] MASS_7.3-27          munsell_0.4.2        plyr_1.8             proto_0.3-10        
[17] RColorBrewer_1.0-5   RCurl_1.95-4.1       reshape2_1.2.2       RSQLite_0.11.4      
[21] scales_0.2.3         splines_3.0.1        stats4_3.0.1         stringr_0.6.2       
[25] survival_2.37-4      tools_3.0.1          XML_3.95-0.2         xtable_1.7-1   



____________________
Darya Vanichkina
PhD Student in Bioinformatics
IMB Science Ambassador
Institute for Molecular Biosciences (IMB)
University of Queensland Building 80
St Lucia Queensland, 4072 Australia
d.vanichkina@gmail.com
d.vanichkina@imb.uq.edu.au




	[[alternative HTML version deleted]]

