In order to investigate the effect of short read mapper used for analysis I started by mapping the short reads with the program subreadalign that can be accessed and started from within R. The mapped reads were then converted to gene counts using functionality from the subread package. These two steps is all that is needed for creating a count table like the one that the sequence center created, using the mapping program tophat and the software htseq-count to convert mapped reads to counts.

Mapping the short reads and creating gene count data

library(Rsubread) # For short read alignment
library(limma) # For detection of differential gene expression
library(edgeR) # For detection of differential gene expression

setwd("/proj/b2015161/private/BILS/reads")
options(digits=4)
targets <- readTargets('../references/targetfile.txt', sep = ' ')
treat <- factor(targets$Treatment)
design <- model.matrix(~0 + treat)

buildindex(basename = "Human", reference = '/pica/data/uppnex/reference/biodata/genomes/Hsapiens/GRCh37/seq/GRCh37.fa') # create index for mapping

align(index = '../alignments/Human', readfile1 = targets$Read1, readfile2 = targets$Read2, input_format = 'gzFASTQ', output_format = 'BAM', output_file = targets$BAM, tieBreakHamming=TRUE, unique=TRUE, indels=5, nthreads = 16)
 # map reads

fc <- featureCounts(files=targets$BAM, annot.ext = "/proj/b2015161/private/BILS/references/Homo_sapiens.GRCh37.75.gtf",isGTFAnnotationFile = TRUE, isPairedEnd = TRUE, strandSpecific = 2, requireBothEndsMapped = FALSE, GTF.featureType="exon", GTF.attrType='gene_id') # count reads

save(list=ls(), file = "Analysis.rdata") # Save the objects outside R

The created objects are below subjected to quality control and filtered for lowly expressed genes.

## Loading required package: limma

## Detection of differential gene expression Perform the detection of differential gene expression for the different contrasts using Limma. The summary of significant genes is seen after commands starting with table, for example table(IEIPW.fit.results) shows 1852 and 1585 down and up-regulated genes for that contrast. The toptable and topTag commands show the 10 genes sorted by adjusted p-values.

design11new <- model.matrix(~0 + treat)
y <- dge.list.filt
y11 <- y[,-11]
y11 <- calcNormFactors(y11)
plotMDS(y11, labels=treat, top=100, gene.selection="common", main = 'MDS plot omitting sample 11', xlim = c(-3, 5), ylim = c(-2,2))

y11.voom <- voom(y11, design = design11new)
fit2 <- lmFit(y11.voom,design11new)

IEIW <- makeContrasts(treatInputEmpty - treatInputWig1, levels = design11new)
IEIW.fit <- contrasts.fit(fit2, IEIW)
IEIW.fit  <- eBayes(IEIW.fit )
IEIW.fit.results <- decideTests(IEIW.fit) 

IEIPE <- makeContrasts(treatInputEmpty - treatIPEmpty, levels = design11new)
IEIPE.fit <- contrasts.fit(fit2, IEIPE)
IEIPE.fit  <- eBayes(IEIPE.fit )
IEIPE.fit.results <- decideTests(IEIPE.fit) 

IEIPW <- makeContrasts(treatInputEmpty - treatIPwig1, levels = design11new)
IEIPW.fit <- contrasts.fit(fit2, IEIPW)
IEIPW.fit  <- eBayes(IEIPW.fit )
IEIPW.fit.results <- decideTests(IEIPW.fit) 

IWIPE <- makeContrasts(treatInputWig1 - treatIPEmpty, levels = design11new)
IWIPE.fit <- contrasts.fit(fit2, IWIPE)
IWIPE.fit  <- eBayes(IWIPE.fit )
IWIPE.fit.results <- decideTests(IWIPE.fit) 

IWIPW <- makeContrasts(treatInputWig1 - treatIPwig1, levels = design11new)
IWIPW.fit <- contrasts.fit(fit2, IWIPW)
IWIPW.fit  <- eBayes(IWIPW.fit )
IWIPW.fit.results <- decideTests(IWIPW.fit) 

IPEIPW <- makeContrasts(treatIPEmpty - treatIPwig1, levels = design11new)
IPEIPW.fit <- contrasts.fit(fit2, IPEIPW)
IPEIPW.fit  <- eBayes(IPEIPW.fit )
IPEIPW.fit.results <- decideTests(IPEIPW.fit) 




table(IEIW.fit.results)
## IEIW.fit.results
##    -1     0 
##     1 13787
topTable(IEIW.fit)
##                          GeneID Length      logFC   AveExpr          t
## ENSG00000172667 ENSG00000172667   9646 -7.7735324 10.188358 -14.027247
## ENSG00000168209 ENSG00000168209   2058  1.0724830  4.229288   5.384567
## ENSG00000139318 ENSG00000139318   4632  0.8260184  4.471247   4.680587
## ENSG00000117152 ENSG00000117152   4836  0.9252529  6.709790   3.649108
## ENSG00000118473 ENSG00000118473  10493  0.6191134  5.237341   3.353563
## ENSG00000171855 ENSG00000171855    840 -0.8232917  6.815111  -3.269578
## ENSG00000186635 ENSG00000186635   9047  0.5104452  5.656320   2.899807
## ENSG00000125740 ENSG00000125740   5553  1.3077836  2.619127   4.086173
## ENSG00000170345 ENSG00000170345   3238  1.3657393  2.592377   4.314694
## ENSG00000123610 ENSG00000123610   1577  0.6835853  6.023115   2.786416
##                      P.Value    adj.P.Val         B
## ENSG00000172667 1.752064e-08 0.0002415746 -3.068738
## ENSG00000168209 2.042565e-04 0.9999913156 -3.313805
## ENSG00000139318 6.294147e-04 0.9999913156 -3.412008
## ENSG00000117152 3.683093e-03 0.9999913156 -3.610157
## ENSG00000118473 6.236511e-03 0.9999913156 -3.716651
## ENSG00000171855 7.251317e-03 0.9999913156 -3.739523
## ENSG00000186635 1.413016e-02 0.9999913156 -3.874800
## ENSG00000125740 1.715540e-03 0.9999913156 -3.878023
## ENSG00000170345 1.160542e-03 0.9999913156 -3.902401
## ENSG00000123610 1.734398e-02 0.9999913156 -3.921700
write.fit(IEIW.fit , results = IEIW.fit.results, file="InputEmpty_InputWig.txt", digits=30, dec=",", adjust = "BH")
table(IEIPE.fit.results)
## IEIPE.fit.results
##    -1     0     1 
##   823 12468   497
topTable(IEIPE.fit)
##                          GeneID Length     logFC  AveExpr          t
## ENSG00000125534 ENSG00000125534   1094  3.115539 6.732394  15.154109
## ENSG00000169372 ENSG00000169372   4276 -3.463410 4.426350 -14.144895
## ENSG00000162627 ENSG00000162627   2088 -2.412483 6.165305 -12.192369
## ENSG00000198839 ENSG00000198839   4014 -2.205874 5.242400 -10.879479
## ENSG00000099284 ENSG00000099284   1980 -2.385552 3.112131 -10.092204
## ENSG00000142173 ENSG00000142173   5386  1.885736 9.524135   9.824737
## ENSG00000116209 ENSG00000116209   4156 -1.945561 8.175756  -9.191037
## ENSG00000108423 ENSG00000108423   3320 -1.966690 3.221831  -9.189460
## ENSG00000151552 ENSG00000151552   3059 -2.088349 5.442071  -9.013825
## ENSG00000087263 ENSG00000087263   5806 -1.778703 7.184488  -8.923925
##                      P.Value    adj.P.Val         B
## ENSG00000125534 7.640511e-09 0.0001053474 10.212078
## ENSG00000169372 1.602205e-08 0.0001104560  9.516481
## ENSG00000162627 7.764500e-08 0.0003568564  8.342711
## ENSG00000198839 2.553170e-07 0.0008800777  7.302206
## ENSG00000099284 5.532328e-07 0.0015255948  6.430517
## ENSG00000142173 7.276264e-07 0.0016720856  6.297705
## ENSG00000116209 1.427586e-06 0.0024646866  5.710689
## ENSG00000108423 1.430047e-06 0.0024646866  5.642257
## ENSG00000151552 1.734998e-06 0.0024835411  5.557206
## ENSG00000087263 1.917618e-06 0.0024835411  5.455910
write.fit(IEIPE.fit , results = IEIPE.fit.results, file="InputEmpty_IPEmpty.txt", digits=30, dec=",", adjust = "BH")
table(IEIPW.fit.results)
## IEIPW.fit.results
##    -1     0     1 
##  1852 10351  1585
topTable(IEIPW.fit)
##                          GeneID Length     logFC   AveExpr         t
## ENSG00000125534 ENSG00000125534   1094  3.841126  6.732394  20.39872
## ENSG00000162627 ENSG00000162627   2088 -2.391360  6.165305 -13.56478
## ENSG00000172667 ENSG00000172667   9646 -7.212755 10.188358 -13.34761
## ENSG00000142173 ENSG00000142173   5386  2.096139  9.524135  12.16914
## ENSG00000169372 ENSG00000169372   4276 -2.640122  4.426350 -11.84749
## ENSG00000072195 ENSG00000072195  17910  2.163511  3.677818  11.73031
## ENSG00000198839 ENSG00000198839   4014 -2.010428  5.242400 -11.05004
## ENSG00000248527 ENSG00000248527    681  4.935322  6.130576  11.05001
## ENSG00000105329 ENSG00000105329   3238  1.824995  6.432629  10.69488
## ENSG00000164136 ENSG00000164136   7860 -1.673335  5.589365 -10.07802
##                      P.Value    adj.P.Val         B
## ENSG00000125534 2.999063e-10 4.135108e-06 12.745946
## ENSG00000162627 2.506612e-08 1.368197e-04  9.445309
## ENSG00000172667 2.976930e-08 1.368197e-04  9.079679
## ENSG00000142173 7.922005e-08 2.677375e-04  8.339786
## ENSG00000169372 1.049784e-07 2.677375e-04  8.099470
## ENSG00000072195 1.165089e-07 2.677375e-04  7.873515
## ENSG00000198839 2.172747e-07 3.744842e-04  7.525376
## ENSG00000248527 2.172812e-07 3.744842e-04  7.435023
## ENSG00000105329 3.047652e-07 4.669003e-04  7.223189
## ENSG00000164136 5.612458e-07 7.595691e-04  6.657347
write.fit(IEIPW.fit , results = IEIPW.fit.results, file="InputEmpty_IPWig1.txt", digits=30, dec=",", adjust = "BH")
table(IWIPE.fit.results)
## IWIPE.fit.results
##    -1     0     1 
##  1013 12129   646
topTable(IWIPE.fit)
##                          GeneID Length     logFC  AveExpr          t
## ENSG00000169372 ENSG00000169372   4276 -3.596559 4.426350 -14.553525
## ENSG00000125534 ENSG00000125534   1094  2.923762 6.732394  14.263820
## ENSG00000162627 ENSG00000162627   2088 -2.395453 6.165305 -12.115923
## ENSG00000198839 ENSG00000198839   4014 -2.237491 5.242400 -11.038274
## ENSG00000142173 ENSG00000142173   5386  2.114231 9.524135  10.953908
## ENSG00000087263 ENSG00000087263   5806 -2.047754 7.184488 -10.289432
## ENSG00000108423 ENSG00000108423   3320 -2.111680 3.221831  -9.728966
## ENSG00000099284 ENSG00000099284   1980 -2.266479 3.112131  -9.733375
## ENSG00000074266 ENSG00000074266   7016 -1.929657 4.874800  -9.570642
## ENSG00000116209 ENSG00000116209   4156 -1.997627 8.175756  -9.441505
##                      P.Value    adj.P.Val        B
## ENSG00000169372 1.180514e-08 0.0001009790 9.746919
## ENSG00000125534 1.464738e-08 0.0001009790 9.724361
## ENSG00000162627 8.295993e-08 0.0003812839 8.297046
## ENSG00000198839 2.196912e-07 0.0006560249 7.443219
## ENSG00000142173 2.378971e-07 0.0006560249 7.279853
## ENSG00000087263 4.537281e-07 0.0010426671 6.783922
## ENSG00000108423 8.038164e-07 0.0013853776 6.143616
## ENSG00000099284 8.001257e-07 0.0013853776 6.139814
## ENSG00000074266 9.493140e-07 0.0014218866 6.113579
## ENSG00000116209 1.089064e-06 0.0014218866 5.963294
write.fit(IWIPE.fit , results = IWIPE.fit.results, file="InputWig1_IPEmpty.txt", digits=30, dec=",", adjust = "BH")
table(IWIPW.fit.results)
## IWIPW.fit.results
##    -1     0     1 
##  1993 10104  1691
topTable(IWIPW.fit)
##                          GeneID Length     logFC  AveExpr          t
## ENSG00000125534 ENSG00000125534   1094  3.649349 6.732394  19.449400
## ENSG00000162627 ENSG00000162627   2088 -2.374329 6.165305 -13.481669
## ENSG00000142173 ENSG00000142173   5386  2.324634 9.524135  13.402609
## ENSG00000169372 ENSG00000169372   4276 -2.773271 4.426350 -12.307134
## ENSG00000198839 ENSG00000198839   4014 -2.042045 5.242400 -11.227427
## ENSG00000109610 ENSG00000109610   2128  2.533445 3.914346  11.075646
## ENSG00000105329 ENSG00000105329   3238  1.800323 6.432629  10.553167
## ENSG00000248527 ENSG00000248527    681  4.699973 6.130576  10.559554
## ENSG00000172340 ENSG00000172340   3090 -1.904865 6.628010 -10.114544
## ENSG00000123595 ENSG00000123595   2039 -1.841310 6.724185  -9.853249
##                      P.Value    adj.P.Val         B
## ENSG00000125534 5.060412e-10 6.977296e-06 12.439462
## ENSG00000162627 2.676342e-08 1.309588e-04  9.405989
## ENSG00000142173 2.849408e-08 1.309588e-04  9.209129
## ENSG00000169372 7.034587e-08 2.424822e-04  8.444814
## ENSG00000198839 1.841175e-07 4.874356e-04  7.684993
## ENSG00000109610 2.121130e-07 4.874356e-04  7.377635
## ENSG00000105329 3.497403e-07 6.027775e-04  7.102541
## ENSG00000248527 3.475658e-07 6.027775e-04  7.031283
## ENSG00000172340 5.408598e-07 6.854789e-04  6.696125
## ENSG00000123595 7.064782e-07 6.854789e-04  6.444504
write.fit(IWIPW.fit , results = IWIPW.fit.results, file="InputWig1_IPWig.txt", digits=30, dec=",", adjust = "BH")
table(IPEIPW.fit.results)
## IPEIPW.fit.results
##     0 
## 13788
topTable(IPEIPW.fit)
##                          GeneID Length      logFC   AveExpr         t
## ENSG00000171855 ENSG00000171855    840 -1.4762232  6.815111 -5.165231
## ENSG00000172667 ENSG00000172667   9646 -4.3706310 10.188358 -7.222277
## ENSG00000163170 ENSG00000163170   3156 -1.7248105  5.903324 -4.736378
## ENSG00000097033 ENSG00000097033   6827 -0.8750023  8.038958 -4.869088
## ENSG00000023697 ENSG00000023697   2848 -1.2461308  3.939642 -5.806371
## ENSG00000171792 ENSG00000171792   2580 -0.8853525  5.156047 -4.464726
## ENSG00000144566 ENSG00000144566   3785 -0.9273984  6.113092 -4.309108
## ENSG00000215414 ENSG00000215414    741 -0.8798991  6.273938 -4.241769
## ENSG00000147533 ENSG00000147533   2474 -0.7985502  6.423738 -4.027670
## ENSG00000107537 ENSG00000107537   3381 -0.7603143  5.327365 -4.006207
##                      P.Value adj.P.Val         B
## ENSG00000171855 2.878769e-04 0.9389942 -3.800804
## ENSG00000172667 1.494781e-05 0.2061004 -3.840138
## ENSG00000163170 5.742562e-04 0.9389942 -3.840391
## ENSG00000097033 4.624960e-04 0.9389942 -3.863082
## ENSG00000023697 1.076174e-04 0.4946097 -3.865630
## ENSG00000171792 9.011010e-04 0.9389942 -3.886878
## ENSG00000144566 1.171590e-03 0.9389942 -3.894474
## ENSG00000215414 1.313784e-03 0.9389942 -3.902591
## ENSG00000147533 1.897987e-03 0.9389942 -3.941311
## ENSG00000107537 1.969879e-03 0.9389942 -3.944081
write.fit(IPEIPW.fit , results = IPEIPW.fit.results, file="IPEmpty_IPWig.txt", digits=30, dec=",", adjust = "BH")

Do the same analysis, but instead using the package edgeR. The two packages differs mainly in the way they estimate the expected variance in expression. The package most suitable for analysis depends on the project and data, but a general pattern seem to be that edgeR tend to have higher proportions of false positives, hence genes that according to the analysis show significant difference among samples, but actually are not different. The Limma pipeline is on the other hand more conservative and might miss genes that actually are differentially expressed among treatments, but tend to show good balance between false positives and false negatives.

isexpr <- rowSums(cpm(fc$counts) > 2) >= 3 # create filter to remove all lowly expressed genes
dtemp <- fc$counts[isexpr,] # Do the actual filtering
d <- DGEList(counts = dtemp[,-11], group=treat)
d <- calcNormFactors(d)
d <- estimateCommonDisp(d, verbose=T)
## Disp = 0.08319 , BCV = 0.2884
d <- estimateTagwiseDisp(d)
de.tgw.InpEmpInpW <- exactTest(d, pair = c('InputEmpty', 'InputWig1'))
de.tgw.InpEmpIPE <- exactTest(d, pair = c('InputEmpty', 'IPEmpty'))
de.tgw.InpEmpIPW <- exactTest(d, pair = c('InputEmpty', 'IPwig1'))
de.tgw.InpWIPE <- exactTest(d, pair = c('InputWig1', 'IPEmpty'))
de.tgw.InpWIPW <- exactTest(d, pair = c('InputWig1', 'IPwig1'))
de.tgw.IPEIPW <- exactTest(d, pair = c('IPEmpty', 'IPwig1'))

table(dt.InpEmpInpW <- decideTestsDGE(de.tgw.InpEmpInpW, p.value=0.05))
## 
##     0     1 
## 13787     1
topTags(de.tgw.InpEmpInpW)
## Comparison of groups:  InputWig1-InputEmpty 
##                      logFC    logCPM       PValue          FDR
## ENSG00000172667  7.7800887 12.306791 3.333628e-38 4.596407e-34
## ENSG00000168209 -1.0488545  4.328088 1.211521e-04 7.191771e-01
## ENSG00000170345 -1.3804419  2.789281 1.564789e-04 7.191771e-01
## ENSG00000125740 -1.2883977  2.853504 1.080100e-03 1.000000e+00
## ENSG00000139318 -0.8273821  4.537784 1.289155e-03 1.000000e+00
## ENSG00000182393  1.5275807  0.956089 1.610157e-03 1.000000e+00
## ENSG00000117152 -0.8985653  6.857214 2.312903e-03 1.000000e+00
## ENSG00000166530 -2.7247388  2.085801 3.686065e-03 1.000000e+00
## ENSG00000157551 -1.3168753  1.708522 3.851858e-03 1.000000e+00
## ENSG00000171855  0.8391425  7.092331 4.695681e-03 1.000000e+00
table(decideTestsDGE(de.tgw.InpEmpIPE, p.value=0.05))
## 
##    -1     0     1 
##   709 12175   904
topTags(de.tgw.InpEmpIPE)
## Comparison of groups:  IPEmpty-InputEmpty 
##                     logFC   logCPM       PValue          FDR
## ENSG00000169372  3.462902 5.216281 7.023879e-29 9.684525e-25
## ENSG00000125534 -3.117814 7.556400 1.291638e-20 8.904552e-17
## ENSG00000162627  2.418004 6.643013 2.330737e-17 1.071207e-13
## ENSG00000230673  2.940244 7.054806 2.521535e-15 7.985705e-12
## ENSG00000198839  2.208159 5.626204 2.895889e-15 7.985705e-12
## ENSG00000238072  2.829830 6.031524 3.079824e-14 7.077436e-11
## ENSG00000099284  2.373717 3.487516 5.039549e-13 9.926472e-10
## ENSG00000183291  2.168949 8.983088 8.537350e-13 1.165530e-09
## ENSG00000223485  2.435668 4.425625 8.585209e-13 1.165530e-09
## ENSG00000248527 -3.836341 7.401819 8.973650e-13 1.165530e-09
table(dt.InpEmpIPE <- decideTestsDGE(de.tgw.InpEmpIPE, p.value=0.05))
## 
##    -1     0     1 
##   709 12175   904
topTags(de.tgw.InpEmpIPE)
## Comparison of groups:  IPEmpty-InputEmpty 
##                     logFC   logCPM       PValue          FDR
## ENSG00000169372  3.462902 5.216281 7.023879e-29 9.684525e-25
## ENSG00000125534 -3.117814 7.556400 1.291638e-20 8.904552e-17
## ENSG00000162627  2.418004 6.643013 2.330737e-17 1.071207e-13
## ENSG00000230673  2.940244 7.054806 2.521535e-15 7.985705e-12
## ENSG00000198839  2.208159 5.626204 2.895889e-15 7.985705e-12
## ENSG00000238072  2.829830 6.031524 3.079824e-14 7.077436e-11
## ENSG00000099284  2.373717 3.487516 5.039549e-13 9.926472e-10
## ENSG00000183291  2.168949 8.983088 8.537350e-13 1.165530e-09
## ENSG00000223485  2.435668 4.425625 8.585209e-13 1.165530e-09
## ENSG00000248527 -3.836341 7.401819 8.973650e-13 1.165530e-09
table(dt.InpEmpIPW <- decideTestsDGE(de.tgw.InpEmpIPW, p.value=0.05))
## 
##    -1     0     1 
##  1805 10231  1752
topTags(de.tgw.InpEmpIPW)
## Comparison of groups:  IPwig1-InputEmpty 
##                     logFC    logCPM       PValue          FDR
## ENSG00000125534 -3.840993  7.556400 2.711620e-39 3.738782e-35
## ENSG00000172667  7.316872 12.306791 2.612298e-35 1.800918e-31
## ENSG00000248527 -4.977544  7.401819 1.483309e-25 6.817289e-22
## ENSG00000169372  2.643952  5.216281 4.604673e-20 1.587231e-16
## ENSG00000162627  2.395053  6.643013 4.656188e-19 1.194034e-15
## ENSG00000198899 -4.243556 10.240415 5.195970e-19 1.194034e-15
## ENSG00000185885 -3.129355  9.345009 3.485423e-16 6.865288e-13
## ENSG00000142173 -2.092929  9.879650 9.598872e-16 1.654366e-12
## ENSG00000161956  2.820210  5.631380 1.676988e-15 2.569145e-12
## ENSG00000109610 -2.303777  4.298497 2.584879e-15 3.564032e-12
table(dt.InpWIPE <- decideTestsDGE(de.tgw.InpWIPE, p.value=0.05))
## 
##    -1     0     1 
##   887 11808  1093
topTags(de.tgw.InpWIPE)
## Comparison of groups:  IPEmpty-InputWig1 
##                     logFC   logCPM       PValue          FDR
## ENSG00000169372  3.608220 5.216281 7.188037e-31 9.910865e-27
## ENSG00000125534 -2.941265 7.556400 7.832703e-19 5.399866e-15
## ENSG00000162627  2.402222 6.643013 3.663452e-17 1.683723e-13
## ENSG00000230673  2.999630 7.054806 7.836602e-16 2.701277e-12
## ENSG00000198839  2.242722 5.626204 1.165705e-15 3.214547e-12
## ENSG00000223485  2.705084 4.425625 3.243190e-15 7.452852e-12
## ENSG00000182004  3.031265 6.312544 4.652464e-15 9.164025e-12
## ENSG00000197632  3.099726 8.911263 5.642327e-15 9.724550e-12
## ENSG00000183291  2.362537 8.983088 8.488754e-15 1.300477e-11
## ENSG00000087263  2.078172 7.446512 3.130317e-13 4.316081e-10
table(dt.InpWIPW <- decideTestsDGE(de.tgw.InpWIPW, p.value=0.05))
## 
##   -1    0    1 
## 1919 9968 1901
Allsign.InpWIPW <- topTags(de.tgw.InpWIPW, n = 1919 + 1901)
Allsign.InpWIPW.up <- Allsign.InpWIPW[Allsign.InpWIPW$table$logFC>0,]
length(Allsign.InpWIPW.up$table$logFC)
## [1] 1901
table(dt.IPEIPW <- decideTestsDGE(de.tgw.IPEIPW, p.value=0.05))
## 
##     0     1 
## 13784     4
topTags(de.tgw.IPEIPW)
## Comparison of groups:  IPwig1-IPEmpty 
##                     logFC    logCPM       PValue          FDR
## ENSG00000229991  3.263247  3.739362 1.130816e-11 1.559169e-07
## ENSG00000172667  3.972725 12.306791 4.581713e-11 3.158633e-07
## ENSG00000163170  1.840479  6.198130 2.495040e-06 1.146721e-02
## ENSG00000171855  1.529618  7.092331 1.112188e-05 3.833713e-02
## ENSG00000128965 -1.886567  2.874621 3.270113e-05 9.017664e-02
## ENSG00000213630  1.741377  3.017537 5.586058e-05 1.213462e-01
## ENSG00000023697  1.253029  4.038106 6.160600e-05 1.213462e-01
## ENSG00000166924 -2.062982  1.575933 2.031494e-04 3.501279e-01
## ENSG00000125740 -1.515291  2.853504 4.530841e-04 6.941249e-01
## ENSG00000121211  1.451095  3.128092 6.116670e-04 6.942759e-01

Below is the same analysis as described above, but for the count data matrix that was delivered from the sequencing facility eg. obained with the softwares Tophat, Cufflinks and HTseq-count.

d.th <- read.table('count_table.txt', header = T)
isexpr <- rowSums(cpm(d.th) > 2) >= 3 # create filter to remove all lowly expressed genes
d.th <- d.th[isexpr,] # Do the actual filtering
d.th <- DGEList(counts = d.th[,-11], group=treat)
d.th <- calcNormFactors(d.th)
y.th.dge <- d.th

y.th.dge.voom <- voom(y.th.dge, design = design11new)
fit.th <- lmFit(y.th.dge.voom, design11new)

IEIW.th <- makeContrasts(treatInputEmpty - treatInputWig1, levels = design11new)
IEIW.fit.th <- contrasts.fit(fit.th, IEIW)
IEIW.fit.th  <- eBayes(IEIW.fit.th )
IEIW.fit.th.results <- decideTests(IEIW.fit.th) 

IEIPE.th <- makeContrasts(treatInputEmpty - treatIPEmpty, levels = design11new)
IEIPE.fit.th <- contrasts.fit(fit.th, IEIPE)
IEIPE.fit.th  <- eBayes(IEIPE.fit.th )
IEIPE.fit.th.results <- decideTests(IEIPE.fit.th) 

IEIPW.th <- makeContrasts(treatInputEmpty - treatIPwig1, levels = design11new)
IEIPW.fit.th <- contrasts.fit(fit.th, IEIPW)
IEIPW.fit.th  <- eBayes(IEIPW.fit.th )
IEIPW.fit.th.results <- decideTests(IEIPW.fit.th) 

IWIPE.th <- makeContrasts(treatInputWig1 - treatIPEmpty, levels = design11new)
IWIPE.fit.th <- contrasts.fit(fit.th, IWIPE.th)
IWIPE.fit.th  <- eBayes(IWIPE.fit.th)
IWIPE.fit.th.results <- decideTests(IWIPE.fit.th) 

IWIPW.th <- makeContrasts(treatInputWig1 - treatIPwig1, levels = design11new)
IWIPW.fit.th <- contrasts.fit(fit.th, IWIPW)
IWIPW.fit.th  <- eBayes(IWIPW.fit.th )
IWIPW.fit.th.results <- decideTests(IWIPW.fit.th)

IPEIPW.th <- makeContrasts(treatIPEmpty - treatIPwig1, levels = design11new)
IPEIPW.fit.th <- contrasts.fit(fit.th, IPEIPW)
IPEIPW.fit.th  <- eBayes(IPEIPW.fit.th )
IPEIPW.fit.th.results <- decideTests(IPEIPW.fit.th)


table(IEIW.fit.th.results)
## IEIW.fit.th.results
##    -1     0 
##     1 12919
topTable(IEIW.fit.th, adjust="BH")
##                      logFC  AveExpr          t      P.Value    adj.P.Val
## ENSG00000172667 -7.6014767 9.958393 -13.505572 2.000376e-08 0.0002584486
## ENSG00000168209  1.1156790 4.234876   5.548780 1.452698e-04 0.9384431274
## ENSG00000139318  0.8379184 4.503332   4.658588 6.120829e-04 0.9999977411
## ENSG00000117152  0.9236464 6.751177   3.752327 2.946046e-03 0.9999977411
## ENSG00000171855 -0.8272344 6.881223  -3.294251 6.727912e-03 0.9999977411
## ENSG00000118473  0.6080325 5.281167   3.256130 7.210924e-03 0.9999977411
## ENSG00000125740  1.2873578 2.638734   4.317988 1.092043e-03 0.9999977411
## ENSG00000123610  0.7031495 6.063551   2.829159 1.571824e-02 0.9999977411
## ENSG00000175197  0.6220436 6.637415   2.743375 1.838180e-02 0.9999977411
## ENSG00000170345  1.3617356 2.605426   4.048522 1.744305e-03 0.9999977411
##                         B
## ENSG00000172667 -3.033400
## ENSG00000168209 -3.271916
## ENSG00000139318 -3.400441
## ENSG00000117152 -3.558099
## ENSG00000171855 -3.715948
## ENSG00000118473 -3.737970
## ENSG00000125740 -3.826043
## ENSG00000123610 -3.894895
## ENSG00000175197 -3.927305
## ENSG00000170345 -3.928455
write.fit(IEIW.fit.th , results = IEIW.fit.th.results, file="InputEmpty_InputWig1_th.txt", digits=30, dec=",", adjust = "BH")
table(IEIPE.fit.th.results)
## IEIPE.fit.th.results
##    -1     0     1 
##   834 11697   389
topTable(IEIPE.fit.th)
##                     logFC  AveExpr          t      P.Value    adj.P.Val
## ENSG00000125534  3.106428 6.656176  14.916763 6.729792e-09 5.887829e-05
## ENSG00000169372 -3.501662 4.536208 -14.511096 9.114286e-09 5.887829e-05
## ENSG00000162627 -2.300675 6.603339 -11.329797 1.330955e-07 5.731982e-04
## ENSG00000198839 -2.264621 5.319481 -10.958428 1.896898e-07 6.126980e-04
## ENSG00000099284 -2.447185 3.197093 -10.426900 3.206728e-07 8.286186e-04
## ENSG00000230673 -3.456785 3.180767 -10.100349 4.476606e-07 9.639625e-04
## ENSG00000151846 -2.126098 5.855516  -9.817332 6.020533e-07 1.079192e-03
## ENSG00000142173  1.898211 9.550370   9.719257 6.682302e-07 1.079192e-03
## ENSG00000074266 -1.807898 4.921036  -9.456484 8.873559e-07 1.273849e-03
## ENSG00000128394 -1.667316 5.664676  -9.230723 1.137888e-06 1.282773e-03
##                         B
## ENSG00000125534 10.422245
## ENSG00000169372 10.082455
## ENSG00000162627  7.917949
## ENSG00000198839  7.602792
## ENSG00000099284  6.953157
## ENSG00000230673  6.562441
## ENSG00000151846  6.553330
## ENSG00000142173  6.397421
## ENSG00000074266  6.186792
## ENSG00000128394  5.961508
write.fit(IEIPE.fit.th , results = IEIPE.fit.th.results, file="InputEmpty_IPEmpty_th.txt", digits=30, dec=",", adjust = "BH")
table(IEIPW.fit.th.results)
## IEIPW.fit.th.results
##   -1    0    1 
## 1792 9678 1450
topTable(IEIPW.fit.th)
##                     logFC  AveExpr         t      P.Value    adj.P.Val
## ENSG00000125534  3.886828 6.656176  20.27056 2.207242e-10 2.851757e-06
## ENSG00000172667 -6.913305 9.958393 -12.66804 4.010585e-08 1.748346e-04
## ENSG00000162627 -2.217805 6.603339 -12.30519 5.489929e-08 1.748346e-04
## ENSG00000142173  2.124637 9.550370  12.11520 6.491854e-08 1.748346e-04
## ENSG00000169372 -2.648386 4.536208 -12.06871 6.766046e-08 1.748346e-04
## ENSG00000198839 -2.024408 5.319481 -10.91828 1.972166e-07 3.346588e-04
## ENSG00000072195  2.108774 3.724948  11.07922 1.688557e-07 3.346588e-04
## ENSG00000105329  1.894794 6.405344  10.86740 2.072191e-07 3.346588e-04
## ENSG00000164136 -1.693613 5.602915 -10.24027 3.876205e-07 5.447184e-04
## ENSG00000109610  2.284028 3.890680  10.15841 4.216086e-07 5.447184e-04
##                         B
## ENSG00000125534 13.023608
## ENSG00000172667  8.874659
## ENSG00000162627  8.776575
## ENSG00000142173  8.518442
## ENSG00000169372  8.500336
## ENSG00000198839  7.614157
## ENSG00000072195  7.577859
## ENSG00000105329  7.575131
## ENSG00000164136  6.997421
## ENSG00000109610  6.756192
write.fit(IEIPW.fit.th , results = IEIPW.fit.th.results, file="InputEmpty_IPWig1_th.txt", digits=30, dec=",", adjust = "BH")

table(IWIPE.fit.th.results)
## IWIPE.fit.th.results
##    -1     0     1 
##  1063 11317   540
topTable(IWIPE.fit.th)
##                     logFC  AveExpr          t      P.Value    adj.P.Val
## ENSG00000169372 -3.573327 4.536208 -14.759388 7.563303e-09 9.289783e-05
## ENSG00000125534  2.889451 6.656176  13.919885 1.438047e-08 9.289783e-05
## ENSG00000198839 -2.265788 5.319481 -10.978884 1.859749e-07 5.735727e-04
## ENSG00000162627 -2.211952 6.603339 -10.905585 1.996620e-07 5.735727e-04
## ENSG00000142173  2.121289 9.550370  10.797017 2.219709e-07 5.735727e-04
## ENSG00000151846 -2.211035 5.855516 -10.205313 4.017588e-07 7.512424e-04
## ENSG00000074266 -1.956266 4.921036 -10.192643 4.070199e-07 7.512424e-04
## ENSG00000099284 -2.281672 3.197093  -9.922324 5.389520e-07 8.452425e-04
## ENSG00000266501 -2.306136 2.853423  -9.838382 5.887912e-07 8.452425e-04
## ENSG00000119820 -1.905582 6.539471  -9.600261 7.592278e-07 8.667468e-04
##                         B
## ENSG00000169372 10.189374
## ENSG00000125534  9.791031
## ENSG00000198839  7.609621
## ENSG00000162627  7.544196
## ENSG00000142173  7.357339
## ENSG00000151846  6.917458
## ENSG00000074266  6.894948
## ENSG00000099284  6.512405
## ENSG00000266501  6.347446
## ENSG00000119820  6.329909
write.fit(IWIPE.fit.th , results = IWIPE.fit.th.results, file="InputWig1_IPempty_th.txt", digits=30, dec=",", adjust = "BH")

table(IWIPW.fit.th.results)
## IWIPW.fit.th.results
##   -1    0    1 
## 1920 9465 1535
topTable(IWIPW.fit.th)
##                     logFC  AveExpr          t      P.Value    adj.P.Val
## ENSG00000125534  3.669851 6.656176  19.212310 4.031278e-10 5.208411e-06
## ENSG00000142173  2.347714 9.550370  13.288937 2.385763e-08 1.541203e-04
## ENSG00000169372 -2.720051 4.536208 -12.346050 5.297107e-08 2.281287e-04
## ENSG00000162627 -2.129082 6.603339 -11.830430 8.382235e-08 2.707462e-04
## ENSG00000109610  2.566498 3.890680  11.474043 1.162900e-07 3.004934e-04
## ENSG00000198839 -2.025575 5.319481 -10.942909 1.925615e-07 4.146491e-04
## ENSG00000105329  1.876883 6.405344  10.765464 2.289483e-07 4.225732e-04
## ENSG00000074266 -1.742531 4.921036 -10.062507 4.655700e-07 5.549904e-04
## ENSG00000119782 -2.304478 4.334525  -9.927661 5.359409e-07 5.549904e-04
## ENSG00000114942 -1.930030 5.758868  -9.859987 5.755058e-07 5.549904e-04
##                         B
## ENSG00000125534 12.723842
## ENSG00000142173  9.405346
## ENSG00000169372  8.747402
## ENSG00000162627  8.427206
## ENSG00000109610  7.914555
## ENSG00000198839  7.658131
## ENSG00000105329  7.502074
## ENSG00000074266  6.827224
## ENSG00000119782  6.669908
## ENSG00000114942  6.637663
write.fit(IWIPW.fit.th , results = IWIPW.fit.th.results, file="InputWig1_IPWig1_th.txt", digits=30, dec=",", adjust = "BH")

table(IPEIPW.fit.th.results)
## IPEIPW.fit.th.results
##     0 
## 12920
topTable(IPEIPW.fit.th)
##                      logFC  AveExpr         t      P.Value adj.P.Val
## ENSG00000171855 -1.4420283 6.881223 -5.045110 3.233552e-04 0.9113719
## ENSG00000163170 -1.6508790 6.560287 -4.855031 4.414146e-04 0.9113719
## ENSG00000172667 -4.2345144 9.958393 -6.945716 1.898241e-05 0.2452528
## ENSG00000097033 -0.8602341 8.087934 -4.729462 5.436608e-04 0.9113719
## ENSG00000023697 -1.2143693 3.993409 -5.501536 1.563543e-04 0.9113719
## ENSG00000171792 -0.8704547 5.161537 -4.146943 1.468605e-03 0.9113719
## ENSG00000089063 -0.7226579 6.828379 -4.021183 1.830048e-03 0.9113719
## ENSG00000215414 -0.8972707 5.853095 -3.956997 2.048945e-03 0.9113719
## ENSG00000147364 -1.0802656 5.199331 -4.058327 1.714579e-03 0.9113719
## ENSG00000169756 -0.9921721 4.788562 -4.573527 7.062909e-04 0.9113719
##                         B
## ENSG00000171855 -3.944951
## ENSG00000163170 -3.954380
## ENSG00000172667 -3.968888
## ENSG00000097033 -4.002187
## ENSG00000023697 -4.009657
## ENSG00000171792 -4.043342
## ENSG00000089063 -4.046614
## ENSG00000215414 -4.054754
## ENSG00000147364 -4.055667
## ENSG00000169756 -4.066746
write.fit(IPEIPW.fit.th , results = IPEIPW.fit.th.results, file="IPEmpty_IPWig1_th.txt", digits=30, dec=",", adjust = "BH")


# edgeR             
d.th <- estimateCommonDisp(d.th, verbose=T)
## Disp = 0.07707 , BCV = 0.2776
d.th <- estimateTagwiseDisp(d.th)

de.th.tgw.InpEmpInpW <- exactTest(d.th, pair = c('InputEmpty', 'InputWig1'))
de.th.tgw.InpEmpIPE <- exactTest(d.th, pair = c('InputEmpty', 'IPEmpty'))
de.th.tgw.InpEmpIPW <- exactTest(d.th, pair = c('InputEmpty', 'IPwig1'))
de.th.tgw.InpWIPE <- exactTest(d.th, pair = c('InputWig1', 'IPEmpty'))
de.th.tgw.InpWIPW <- exactTest(d.th, pair = c('InputWig1', 'IPwig1'))
de.th.tgw.IPEIPW <- exactTest(d.th, pair = c('IPEmpty', 'IPwig1'))

table(dt.th.InpEmpInpW <- decideTestsDGE(de.th.tgw.InpEmpInpW, p.value=0.05))
## 
##     0     1 
## 12919     1
topTags(de.th.tgw.InpEmpInpW)
## Comparison of groups:  InputWig1-InputEmpty 
##                      logFC     logCPM       PValue          FDR
## ENSG00000172667  7.6114640 12.0268643 1.710634e-35 2.210139e-31
## ENSG00000168209 -1.0934748  4.3296236 4.410470e-05 2.849164e-01
## ENSG00000170345 -1.3748118  2.7957053 1.860548e-04 8.012762e-01
## ENSG00000125740 -1.2488131  2.8650467 6.965599e-04 1.000000e+00
## ENSG00000139318 -0.8399987  4.5680984 8.670483e-04 1.000000e+00
## ENSG00000232561  2.6726201  0.6026312 9.976134e-04 1.000000e+00
## ENSG00000117152 -0.8995946  6.8916780 1.820369e-03 1.000000e+00
## ENSG00000182393  1.5073256  1.0415123 2.222620e-03 1.000000e+00
## ENSG00000203836  3.5199434  3.9060960 3.126553e-03 1.000000e+00
## ENSG00000171855  0.8408431  7.1562351 4.369952e-03 1.000000e+00
table(dt.th.tgw.InpEmpIPE <- decideTestsDGE(de.th.tgw.InpEmpIPE, p.value=0.05))
## 
##    -1     0     1 
##   537 11416   967
topTags(de.th.tgw.InpEmpIPE)
## Comparison of groups:  IPEmpty-InputEmpty 
##                     logFC   logCPM       PValue          FDR
## ENSG00000169372  3.497497 5.318430 1.299192e-30 1.678556e-26
## ENSG00000125534 -3.107520 7.480775 5.164929e-21 3.336544e-17
## ENSG00000230673  3.402008 3.932917 2.512902e-19 1.082223e-15
## ENSG00000198839  2.268150 5.705942 1.787563e-16 5.773829e-13
## ENSG00000162627  2.305746 7.011228 3.384680e-16 8.746013e-13
## ENSG00000238072  2.933869 6.214244 1.987418e-15 4.279574e-12
## ENSG00000099284  2.433501 3.571900 1.051960e-14 1.941617e-11
## ENSG00000151846  2.112800 6.183278 5.076542e-14 8.198616e-11
## ENSG00000223485  2.520562 4.474278 7.479566e-14 1.073733e-10
## ENSG00000197632  2.732914 8.847440 8.599917e-14 1.111109e-10
table(dt.th.InpEmpIPW <- decideTestsDGE(de.th.tgw.InpEmpIPW, p.value=0.05))
## 
##   -1    0    1 
## 1600 9561 1759
topTags(de.th.tgw.InpEmpIPW)
## Comparison of groups:  IPwig1-InputEmpty 
##                     logFC    logCPM       PValue          FDR
## ENSG00000125534 -3.886421  7.480775 4.843075e-41 6.257253e-37
## ENSG00000172667  7.068171 12.026864 2.919337e-32 1.885892e-28
## ENSG00000248527 -4.958177  7.427534 8.924122e-22 3.843322e-18
## ENSG00000169372  2.649756  5.318430 5.824553e-21 1.881331e-17
## ENSG00000198899 -4.302465 10.172689 2.326540e-19 6.011780e-16
## ENSG00000185885 -3.394676  8.408570 5.206057e-19 1.121038e-15
## ENSG00000230673  2.996976  3.932917 4.531128e-17 8.363167e-14
## ENSG00000162627  2.211002  7.011228 6.999511e-17 1.130421e-13
## ENSG00000142173 -2.119629  9.906423 4.126640e-16 5.924021e-13
## ENSG00000109610 -2.288673  4.281133 6.800895e-16 8.786757e-13
table(dt.th.InpWIPE <- decideTestsDGE(de.th.tgw.InpWIPE, p.value=0.05))
## 
##    -1     0     1 
##   679 11090  1151
topTags(de.th.tgw.InpWIPE)
## Comparison of groups:  IPEmpty-InputWig1 
##                     logFC   logCPM       PValue          FDR
## ENSG00000169372  3.586875 5.318430 6.429400e-32 8.306785e-28
## ENSG00000197632  3.328251 8.847440 4.341833e-19 2.762529e-15
## ENSG00000125534 -2.905013 7.480775 6.414542e-19 2.762529e-15
## ENSG00000223485  2.824222 4.474278 1.163877e-16 3.759323e-13
## ENSG00000198839  2.273232 5.705942 1.589763e-16 4.107948e-13
## ENSG00000230673  3.026928 3.932917 5.101383e-16 1.098498e-12
## ENSG00000183291  2.420663 9.030230 1.780250e-15 3.285834e-12
## ENSG00000162627  2.223271 7.011228 3.281501e-15 5.299624e-12
## ENSG00000151846  2.208773 6.183278 3.943816e-15 5.661566e-12
## ENSG00000179766  3.901487 4.379430 9.410110e-15 1.215786e-11
table(dt.th.InpWIPW <- decideTestsDGE(de.th.tgw.InpWIPW, p.value=0.05))
## 
##   -1    0    1 
## 1726 9316 1878
topTags(de.th.tgw.InpWIPW)
## Comparison of groups:  IPwig1-InputWig1 
##                     logFC    logCPM       PValue          FDR
## ENSG00000125534 -3.683915  7.480775 1.105879e-37 1.428796e-33
## ENSG00000169372  2.739130  5.318430 4.068472e-22 2.628233e-18
## ENSG00000248527 -4.669473  7.427534 5.608738e-20 2.415496e-16
## ENSG00000109610 -2.576787  4.281133 1.627193e-19 5.255833e-16
## ENSG00000142173 -2.342656  9.906423 4.719358e-19 1.021874e-15
## ENSG00000224411 -3.220785  3.769173 4.745543e-19 1.021874e-15
## ENSG00000198899 -4.070600 10.172689 7.616634e-18 1.405813e-14
## ENSG00000185885 -3.190953  8.408570 2.958393e-17 4.777805e-14
## ENSG00000123975  2.603033  8.539223 3.775211e-16 4.962558e-13
## ENSG00000119782  2.329777  4.719924 3.840989e-16 4.962558e-13
table(dt.th.IPEIPW <- decideTestsDGE(de.th.tgw.IPEIPW, p.value=0.05))
## 
##    -1     0     1 
##     1 12916     3
topTags(de.th.tgw.IPEIPW)
## Comparison of groups:  IPwig1-IPEmpty 
##                     logFC    logCPM       PValue          FDR
## ENSG00000172667  3.905469 12.026864 2.287829e-10 2.955875e-06
## ENSG00000128965 -2.038372  2.901964 3.904474e-06 1.879289e-02
## ENSG00000163170  1.747967  6.806657 4.363674e-06 1.879289e-02
## ENSG00000171855  1.497440  7.156235 1.458872e-05 4.712156e-02
## ENSG00000023697  1.221894  4.089261 6.957131e-05 1.614744e-01
## ENSG00000125740 -1.606418  2.865047 7.498813e-05 1.614744e-01
## ENSG00000145736  1.352190  2.641134 2.631911e-04 4.544163e-01
## ENSG00000106538 -1.772663  2.636618 2.813723e-04 4.544163e-01
## ENSG00000165233 -1.453091  3.701268 4.916764e-04 6.755388e-01
## ENSG00000121211  1.551315  3.206208 5.228629e-04 6.755388e-01

Comparing the results between the different methods shows that the analysis is largely congruent between data sets and methods, but there are certainly some differences. As an example below is seen by taking the 1000 smallest p-values in the contrast Input Empty vs. IP Empty mapped with subreadalign and analysed with Limma and EdgeR.

table(row.names(toptable(IEIPE.fit, n = 1000))%in% row.names(topTags(de.tgw.InpEmpIPE, n = 1000))) # Compares the significant genes in the limma analysis to the results with edgeR.
## 
## FALSE  TRUE 
##   221   779

The correponding values for the tophat ht-seq data are:

table(row.names(toptable(IEIPE.fit.th, n = 1000))%in% row.names(topTags(de.th.tgw.InpEmpIPE, n = 1000)))
## 
## FALSE  TRUE 
##   184   816

And comparing the effect of the mapper also shows similar patterns with around 80% of common genes in the 1000 with lowest p-values in both the Limma and EdgeR analysis.

table(row.names(toptable(IEIPE.fit.th, n = 1000))%in% row.names(toptable(IEIPE.fit, n = 1000)))
## 
## FALSE  TRUE 
##   153   847
table(row.names(topTags(de.tgw.InpEmpIPE, n = 1000))%in% row.names(topTags(de.th.tgw.InpEmpIPE, n = 1000)))
## 
## FALSE  TRUE 
##   193   807

Additional filtering and creating new list

From the created lists of significant genes the up-regulated genes in the IP treatment is extracted and the common once that is found in both the empty and the wig contrasts are subtracted to create a list of 1087 genes. This gene list is hence genes that are significantly up-regulated between InputWIG1 - IPWIG1 and not significantly upregulated in InputEmpty - IPEmpty. Since the list of 1087 is fairly long I have also done the same subtractions but setting the FDR rate to 0.01 instead (gives 628 genes). All steps has been done using the subread alignment data and edgeR differential expression analysis.

# p-val = 0.05 that gives 709 + 904 sign. in the first contrast and 1919 + 1901 for the second
Allsign.InpEmpIPE <- topTags(de.tgw.InpEmpIPE, n =709 + 904)
Allsign.InpEmpIPE.Up <- Allsign.InpEmpIPE[Allsign.InpEmpIPE$table$logFC>0,]
Allsign.InpWIPW <- topTags(de.tgw.InpWIPW, n = 1919 + 1901)
Allsign.InpWIPW.up <- Allsign.InpWIPW[Allsign.InpWIPW$table$logFC>0,]
Unique.up.IPW <- Allsign.InpWIPW.up[!row.names(Allsign.InpWIPW.up) %in% row.names(Allsign.InpEmpIPE.Up),]
write.table(Unique.up.IPW, file = 'Unique.up.IPW_pval_0.05.txt', sep = '\t')

# p-val = 0.01
table(dt.InpEmpIPE <- decideTestsDGE(de.tgw.InpEmpIPE, p.value=0.01))
## 
##    -1     0     1 
##   277 12999   512
table(dt.InpEmpIPW <- decideTestsDGE(de.tgw.InpEmpIPW, p.value=0.01))
## 
##    -1     0     1 
##  1020 11705  1063
Allsign.InpEmpIPE <- topTags(de.tgw.InpEmpIPE, n = 277 + 512)
Allsign.InpEmpIPE.Up <- Allsign.InpEmpIPE[Allsign.InpEmpIPE$table$logFC>0,]
Allsign.InpWIPW <- topTags(de.tgw.InpWIPW, n = 1020 + 1063)
Allsign.InpWIPW.up <- Allsign.InpWIPW[Allsign.InpWIPW$table$logFC>0,]
Unique.up.IPW <- Allsign.InpWIPW.up[!row.names(Allsign.InpWIPW.up) %in% row.names(Allsign.InpEmpIPE.Up),]

write.table(Unique.up.IPW, file = 'Unique.up.IPW_pval_0.01.txt')