RNA analysis of Zebra fish Gonads

This document describes contains some basic information from the RNA-seq data from Zebra fish. All commands needed to run the analysis that have been done. Based on these results I can perform downstream analysis of the results.

Data

12 individuals from 3 different treatments (eg 4 replicates/treatment) were sequenced using Proton technology. The output per library is shown in the table below. The column ‘Number of reads’ is the actual macine output, whereas ‘Useful reads’ are the number of reads that are mapping to areas of the genome annotated as genes and hence part of the gene expression data.

Lib Name	Treatment	Replicate	Number of reads	Useful reads
up-87-1	Control	1	19741198	8153032
up-87-2	Control	2	22141691	9051586
up-87-3	Control	3	19379328	8847111
up-87-4	Control	4	20509164	5244292
up-87-5	EE2-3	1	21248836	6671255
up-87-6	EE2-3	2	19692093	9482823
up-87-7	EE2-3	3	19737758	7617444
up-87-8	EE2-3	4	14064762	6310603
up-87-9	EE2-10	1	18336786	10896593
up-87-10	EE2-10	2	14963248	8864047
up-87-11	EE2-10	3	13291577	7756863
up-87-12	EE2-10	4	21512379	7836155

All this data were mapped to the Danio rerio v. 9 and the ensembl annotation were used for extracting gene expression data.

Quality control and mapping

The average quality or reads looks okay, but the mean level is slightly lower than commonly seen for Illumina reads. I think this expected with this technology as this is seen in all projects that I have been involved in that use proton data. I hence see no need to filter reads prior to mapping reads to the references sequence.

Control 1 Control 2 Control 3 Control 4 EE2-3 1 EE2-3 2 EE2-3 3 EE2-3 4 EE2-10 1 EE2-10 2 EE2-10 3 EE2-10 4

The program Star with the following settings.

genome=/home/thomask/Zebra092/nobackup/Danio_rerio/Ensembl/Zv9/Sequence/WholeGenomeFasta/genome.fa
gtf=/home/thomask/Zebra092/nobackup/Danio_rerio/Ensembl/Zv9/Annotation/Genes/genes.gtf
dir=/home/thomask/Zebra1481/private/Reference
Mapping=/home/thomask/Zebra1481/private/Mapping

cd /proj/b2014092/private/up87/

for file in *.fastq
    do
    star --readFilesIn ${file} --genomeDir ${dir} --outFileNamePrefix ~/Zebra1481/private/Mapping/${file} --outSAMtype BAM SortedByCoordinate --runThreadN 16 
    done

Prior to this the reference genome were indexed using with the following code

genome=/home/thomask/Zebra092/nobackup/Danio_rerio/Ensembl/Zv9/Sequence/WholeGenomeFasta/genome.fa
gtf=/home/thomask/Zebra092/nobackup/Danio_rerio/Ensembl/Zv9/Annotation/Genes/genes.gtf
dir=/home/thomask/Zebra1481/private/Reference

cd ${dir}
star --runThreadN 16 --runMode genomeGenerate --genomeDir ${dir} --genomeFastaFiles ${genome} --sjdbGTFfile ${gtf} --sjdbOverhang 100

The mapped reads were converted to read counts per gene with the script htseq-count as follows

gtf=/home/thomask/Zebra092/nobackup/Danio_rerio/Ensembl/Zv9/Annotation/Genes/genes.gtf

cd /home/thomask/Zebra1481/private/Mapping
for file in *.bam
    do
    htseq-count -f bam ${file} ${gtf} > Counts_${file}
    done

From the count data additional quality control steps were performed. First of all correlations among samples and replicates were performed. The average correlation among samples were 0.91 on the unfiltered data and 0.97 when filtering the genes annotated as non-coding. The shift in average correlation is largely due to the control replicate 4 as seen in the plot ‘Correlations among samples, all genes’. This would suggest that filtering the non-coding genes will lead to a data set with less noise.

The same plot, but with non-coding genes clearly shows that sample control 4 is showing a better correlation to the other samples once the non-coding genes. It also highlights that the control 1 does over the coding genes show low correlations to all but one of the other samples.

## Loading required package: limma

## Warning: package 'biomaRt' was built under R version 3.2.2

## Creating a generic function for 'nchar' from package 'base' in package 'S4Vectors'

As genes that have 0 or very few genes mapping to them is not meaningful to analyse for differential gene expression, one often filter genes that are lowly expressed. Below I filter all genes that show less than 1 read per million in more than 4 samples.

This drops 12086 genes from the 33737 that are found full data set. For the data set with only coding genes the corresponding figures are 6352 and 26812. In both data sets we hence end up with similar number of genes after filtering the lowly expressed once, suggesting that very few non-coding genes have high enough expression levels to be retained in the analysis.

All_DGE.keep14 <- rowSums(cpm(All_DGE)>1) >= 4 # create filtered that needs at least 1 cpm in 4 libraries (number of replicates in each treatment) for the genes to be retained in analysis
All_DGE.nc.rm.keep14 <- rowSums(cpm(All_DGE.nc.rm)>1) >= 4 # create filtered that needs at least 1 cpm in 4 libraries (number of replicates in each treatment) for the genes to be retained in analysis
table(All_DGE.keep14) # how many genes are retained after filtering

## All_DGE.keep14
## FALSE  TRUE 
## 12086 21651

table(All_DGE.nc.rm.keep14) # how many genes are retained after filtering

## All_DGE.nc.rm.keep14
## FALSE  TRUE 
##  6352 20460

All_DGE14 <- All_DGE[All_DGE.keep14,] # do the actual filtering
All_DGE.nc.rm14 <- All_DGE.nc.rm[All_DGE.nc.rm.keep14,] # do the actual filtering
levelplot(cor(All_DGE14$counts))

levelplot(cor(All_DGE.nc.rm14$counts))

All_DGE14 <- calcNormFactors(All_DGE14)
All_DGE14 <- estimateGLMCommonDisp(All_DGE14, designModel, verbose=TRUE)

## Disp = 0.17522 , BCV = 0.4186

All_DGE14 <- estimateGLMTrendedDisp(All_DGE14, designModel)
All_DGE14 <- estimateGLMTagwiseDisp(All_DGE14, designModel)

All_DGE.nc.rm14 <- calcNormFactors(All_DGE.nc.rm14)
All_DGE.nc.rm14 <- estimateGLMCommonDisp(All_DGE.nc.rm14, designModel, verbose=TRUE)

## Disp = 0.1739 , BCV = 0.417

All_DGE.nc.rm14 <- estimateGLMTrendedDisp(All_DGE.nc.rm14, designModel)
All_DGE.nc.rm14 <- estimateGLMTagwiseDisp(All_DGE.nc.rm14, designModel)


plotBCV(All_DGE14)

plotBCV(All_DGE.nc.rm14)

All_DGE14_2 <- estimateGLMCommonDisp(All_DGE14, designModel2, verbose=TRUE)

## Disp = 0.09667 , BCV = 0.3109

All_DGE14_2 <- estimateGLMTrendedDisp(All_DGE14_2, designModel2)
All_DGE14_2 <- estimateGLMTagwiseDisp(All_DGE14_2, designModel2)

All_DGE14.fit_2 <- glmFit(All_DGE14_2, designModel2)
All_DGE14.fitEE2_3_2 <- glmLRT(All_DGE14.fit_2, coef = "treatmentEE2_3")
All_DGE14.fitEE2_10_2 <- glmLRT(All_DGE14.fit_2, coef = "treatmentEE2_10")
Table.All_DGE14EE2_3_2 <- topTags(All_DGE14.fitEE2_3_2, n=nrow(All_DGE14))$table
Table.All_DGE14EE2_10_2 <- topTags(All_DGE14.fitEE2_10_2, n=nrow(All_DGE14))$table
table(resultAll_DGE14EE2_3_2 <- decideTestsDGE(All_DGE14.fitEE2_3_2))

## 
##    -1     0     1 
##   104 21233   314

topTags(All_DGE14.fitEE2_3_2)

## Coefficient:  treatmentEE2_3 
##                        logFC    logCPM       LR       PValue          FDR
## ENSDARG00000081077 -3.812903 3.3600861 55.08928 1.151772e-13 2.493701e-09
## ENSDARG00000083266 -2.626477 2.6333895 33.67777 6.503960e-09 5.821616e-05
## ENSDARG00000010425  1.970942 6.1158920 33.10263 8.742055e-09 5.821616e-05
## ENSDARG00000017355  1.786352 6.3312637 31.83307 1.680077e-08 5.821616e-05
## ENSDARG00000094343 -3.724469 1.6574796 31.67200 1.825350e-08 5.821616e-05
## ENSDARG00000003020 -4.259110 0.5645728 31.62265 1.872335e-08 5.821616e-05
## ENSDARG00000086157 -2.027017 3.4880740 31.61245 1.882191e-08 5.821616e-05
## ENSDARG00000026864  1.979879 4.9310027 30.87380 2.753641e-08 7.267176e-05
## ENSDARG00000089970 -2.264131 4.5987832 30.69408 3.020857e-08 7.267176e-05
## ENSDARG00000095009 -2.972865 3.4748892 30.30152 3.698352e-08 7.543707e-05

table(resultAll_DGE14EE2_10_2 <- decideTestsDGE(All_DGE14.fitEE2_10_2))

## 
##    -1     0     1 
##   191 21087   373

topTags(All_DGE14.fitEE2_10_2)

## Coefficient:  treatmentEE2_10 
##                        logFC    logCPM       LR       PValue          FDR
## ENSDARG00000063276 -9.152203 0.3178874 49.75589 1.741148e-12 2.103961e-08
## ENSDARG00000084628 -2.139205 6.4617906 49.54016 1.943523e-12 2.103961e-08
## ENSDARG00000058370  2.073519 5.7853730 40.68043 1.792737e-10 1.293818e-06
## ENSDARG00000087368 -6.274065 3.8797329 40.03263 2.497557e-10 1.351865e-06
## ENSDARG00000080372 -2.634253 3.9352533 38.14269 6.575566e-10 2.847352e-06
## ENSDARG00000003020 -4.553787 0.5645728 36.93602 1.220700e-09 4.404898e-06
## ENSDARG00000080499 -2.196263 3.5362893 36.35339 1.645931e-09 5.090864e-06
## ENSDARG00000086774 -2.904867 4.8129059 35.82169 2.162269e-09 5.851911e-06
## ENSDARG00000089970 -2.316311 4.5987832 35.58127 2.446290e-09 5.884957e-06
## ENSDARG00000085829 -2.566893 4.0341587 35.12465 3.092619e-09 6.695830e-06

All_DGE14.fit <- glmFit(All_DGE14, designModel)
All_DGE14.fitEE2_3 <- glmLRT(All_DGE14.fit, coef = "treatmentEE2_3")
All_DGE14.fitEE2_10 <- glmLRT(All_DGE14.fit, coef = "treatmentEE2_10")
Table.All_DGE14EE2_3 <- topTags(All_DGE14.fitEE2_3, n=nrow(All_DGE14))$table
Table.All_DGE14EE2_10 <- topTags(All_DGE14.fitEE2_10, n=nrow(All_DGE14))$table
table(resultAll_DGE14EE2_3 <- decideTestsDGE(All_DGE14.fitEE2_3))

## 
##    -1     0     1 
##     7 21640     4

topTags(All_DGE14.fitEE2_3)

## Coefficient:  treatmentEE2_3 
##                        logFC     logCPM       LR       PValue          FDR
## ENSDARG00000081077 -3.288345  3.3600861 40.20076 2.291587e-10 4.961514e-06
## ENSDARG00000083266 -2.401551  2.6333895 29.36339 6.000053e-08 6.495357e-04
## ENSDARG00000095618  1.976902  2.7687710 26.05853 3.312222e-07 2.390431e-03
## ENSDARG00000075193  1.871694  3.4215179 24.15760 8.876533e-07 4.804645e-03
## ENSDARG00000091399 -5.037778 -0.1574564 22.59628 1.998715e-06 7.339875e-03
## ENSDARG00000092458 -3.941008  0.8175363 22.56261 2.034052e-06 7.339875e-03
## ENSDARG00000085168 -6.851176  6.4469853 21.64316 3.283782e-06 1.015674e-02
## ENSDARG00000086157 -1.566132  3.4880740 19.48327 1.014845e-05 2.441617e-02
## ENSDARG00000003687  2.558923  1.3867636 19.48309 1.014944e-05 2.441617e-02
## ENSDARG00000083434 -2.311802  2.1263820 18.15711 2.034084e-05 4.403994e-02

table(resultAll_DGE14EE2_10 <- decideTestsDGE(All_DGE14.fitEE2_10))

## 
##    -1     0 
##     1 21650

topTags(All_DGE14.fitEE2_10)

## Coefficient:  treatmentEE2_10 
##                         logFC     logCPM       LR       PValue         FDR
## ENSDARG00000085168  -8.088604  6.4469853 27.03828 1.994657e-07 0.004318632
## ENSDARG00000083042   3.449144  1.2876101 19.85616 8.349346e-06 0.063932314
## ENSDARG00000085829  -2.225485  4.0341587 19.40161 1.059175e-05 0.063932314
## ENSDARG00000063276  -7.164567  0.3178874 18.75954 1.482796e-05 0.063932314
## ENSDARG00000081218 -12.773096 13.1612192 18.31732 1.869993e-05 0.063932314
## ENSDARG00000080404  -2.209441  2.4327117 18.28056 1.906429e-05 0.063932314
## ENSDARG00000082850  -1.720103  4.5786217 18.02625 2.178797e-05 0.063932314
## ENSDARG00000080372  -1.957950  3.9352533 17.87235 2.362286e-05 0.063932314
## ENSDARG00000070000  -4.058072  4.4499515 17.31693 3.163566e-05 0.076104843
## ENSDARG00000080499  -1.750753  3.5362893 16.52946 4.789999e-05 0.103708261

All_DGE.nc.rm14.fit <- glmFit(All_DGE.nc.rm14, designModel)
All_DGE.nc.rm14.fitEE2_3 <- glmLRT(All_DGE.nc.rm14.fit, coef = "treatmentEE2_3")
All_DGE.nc.rm14.fitEE2_10 <- glmLRT(All_DGE.nc.rm14.fit, coef = "treatmentEE2_10")
Table.All_DGE.nc.rm14EE2_3 <- topTags(All_DGE.nc.rm14.fitEE2_3, n=nrow(All_DGE.nc.rm14))$table
Table.All_DGE.nc.rm14EE2_10 <- topTags(All_DGE.nc.rm14.fitEE2_10, n=nrow(All_DGE.nc.rm14))$table
table(resultAll_DGE.nc.rm14EE2_3 <- decideTestsDGE(All_DGE.nc.rm14.fitEE2_3))

## 
##    -1     0     1 
##     2 20457     1

topTags(All_DGE.nc.rm14.fitEE2_3)

## Coefficient:  treatmentEE2_3 
##                        logFC     logCPM       LR       PValue         FDR
## ENSDARG00000095618  1.974458  2.7682458 26.12925 3.193102e-07 0.006533087
## ENSDARG00000091399 -5.040822 -0.1575644 22.67616 1.917318e-06 0.015061744
## ENSDARG00000092458 -3.942801  0.8172499 22.40458 2.208467e-06 0.015061744
## ENSDARG00000003687  2.556425  1.3864546 19.34142 1.093091e-05 0.055911591
## ENSDARG00000041339  5.427560  1.7224635 17.71788 2.562076e-05 0.104840145
## ENSDARG00000062788  4.223931  3.9358022 16.98773 3.762222e-05 0.112064325
## ENSDARG00000070242  4.178573  3.3506308 16.95181 3.834068e-05 0.112064325
## ENSDARG00000058556  3.748632  3.4474707 16.12635 5.925373e-05 0.138706434
## ENSDARG00000070000 -3.848570  4.4497926 15.85823 6.826910e-05 0.138706434
## ENSDARG00000026925  4.178863  4.7479317 15.84352 6.880175e-05 0.138706434

table(resultAll_DGE.nc.rm14EE2_10 <- decideTestsDGE(All_DGE.nc.rm14.fitEE2_10))

## 
##     0 
## 20460

topTags(All_DGE.nc.rm14.fitEE2_10)

## Coefficient:  treatmentEE2_10 
##                        logFC    logCPM       LR       PValue       FDR
## ENSDARG00000063276 -7.169694 0.3182241 18.81848 1.437671e-05 0.2941475
## ENSDARG00000070000 -4.059054 4.4497926 17.34293 3.120577e-05 0.3192351
## ENSDARG00000034058 -3.790933 3.2150098 15.54090 8.073931e-05 0.3426614
## ENSDARG00000045835  6.744297 4.6691785 15.26711 9.332759e-05 0.3426614
## ENSDARG00000087368 -5.346894 3.8798049 14.86499 1.154867e-04 0.3426614
## ENSDARG00000003687  2.221157 1.3864546 14.74638 1.229835e-04 0.3426614
## ENSDARG00000092022  2.195444 4.2116451 14.60121 1.328295e-04 0.3426614
## ENSDARG00000037861  1.464917 3.6483854 14.30059 1.558161e-04 0.3426614
## ENSDARG00000033683 -3.700305 4.7099720 14.26634 1.586772e-04 0.3426614
## ENSDARG00000086739 -2.217343 0.9191311 14.16475 1.674787e-04 0.3426614

Below is the same analysis but omitting the sample Control 1 that seem to show a very different pattern.

All_DGEfiltered <- All_DGE[,c(2:12)] # remove the samples
designfiltered <- design[c(2:12),]
treatmentfiltered <- factor(designfiltered$Treatment)
designModelfiltered <- model.matrix(~treatmentfiltered)

All_DGE.nc.rmfiltered <- All_DGEfiltered
All_DGE.nc.rmfiltered$counts <- All_DGE.nc.rmfiltered$counts[!(row.names(All_DGE.nc.rmfiltered$counts)%in%listncRNA[,1]),]

All_DGE.keep13filtered <- rowSums(cpm(All_DGEfiltered)>1) >= 3 # create filtered that needs at least 1 cpm in 4 libraries (number of replicates in each treatment) for the genes to be retained in analysis
All_DGE.nc.rm.keep13filtered <- rowSums(cpm(All_DGE.nc.rmfiltered)>1) >= 3 # create filtered that needs at least 1 cpm in 4 libraries (number of replicates in each treatment) for the genes to be retained in analysis
table(All_DGE.keep13filtered) # how many genes are retained after filtering

## All_DGE.keep13filtered
## FALSE  TRUE 
## 11674 22063

table(All_DGE.nc.rm.keep13filtered) # how many genes are retained after filtering

## All_DGE.nc.rm.keep13filtered
## FALSE  TRUE 
##  6005 20807

All_DGE13filtered <- All_DGEfiltered[All_DGE.keep13filtered,] # do the actual filtering
All_DGE.nc.rm13filtered <- All_DGE.nc.rmfiltered[All_DGE.nc.rm.keep13filtered,] # do the actual filtering

plotMDS(All_DGE13filtered, dim.plot=c(1,2),labels=labels[c(2:12)], xlim=c(-4,6))

plotMDS(All_DGE.nc.rm13filtered, dim.plot=c(1,2),xlim=c(-4,6))

levelplot(cor(All_DGE13filtered$counts))

levelplot(cor(All_DGE.nc.rm13filtered$counts))

All_DGE13filtered <- calcNormFactors(All_DGE13filtered)
All_DGE13filtered <- estimateGLMCommonDisp(All_DGE13filtered, designModelfiltered, verbose=TRUE)

## Disp = 0.12291 , BCV = 0.3506

All_DGE13filtered <- estimateGLMTrendedDisp(All_DGE13filtered, designModelfiltered)
All_DGE13filtered <- estimateGLMTagwiseDisp(All_DGE13filtered, designModelfiltered)

All_DGE.nc.rm13filtered <- calcNormFactors(All_DGE.nc.rm13filtered)
All_DGE.nc.rm13filtered <- estimateGLMCommonDisp(All_DGE.nc.rm13filtered, designModelfiltered, verbose=TRUE)

## Disp = 0.12392 , BCV = 0.352

All_DGE.nc.rm13filtered <- estimateGLMTrendedDisp(All_DGE.nc.rm13filtered, designModelfiltered)
All_DGE.nc.rm13filtered <- estimateGLMTagwiseDisp(All_DGE.nc.rm13filtered, designModelfiltered)

All_DGE13filtered.fit <- glmFit(All_DGE13filtered, designModelfiltered)
All_DGE13filtered.fitEE2_3 <- glmLRT(All_DGE13filtered.fit, coef = "treatmentfilteredEE2_3")
All_DGE13filtered.fitEE2_10 <- glmLRT(All_DGE13filtered.fit, coef = "treatmentfilteredEE2_10")
Table.All_DGE13filteredEE2_3 <- topTags(All_DGE13filtered.fitEE2_3, n=nrow(All_DGE13filtered))$table
Table.All_DGE13filteredEE2_10 <- topTags(All_DGE13filtered.fitEE2_10, n=nrow(All_DGE13filtered))
table(resultAll_DGE13filteredEE2_3 <- decideTestsDGE(All_DGE13filtered.fitEE2_3))

## 
##    -1     0     1 
##     6 22048     9

topTags(All_DGE13filtered.fitEE2_3, n = 16)

## Coefficient:  treatmentfilteredEE2_3 
##                         logFC     logCPM       LR       PValue
## ENSDARG00000081077  -3.259998  3.2903327 36.34908 1.649571e-09
## ENSDARG00000085168  -7.244447  6.6105761 29.55016 5.448806e-08
## ENSDARG00000083266  -2.326151  2.5375884 25.65639 4.079421e-07
## ENSDARG00000092458  -4.182849  0.8985077 25.25233 5.029865e-07
## ENSDARG00000041339   7.412487  1.8527827 22.58908 2.006216e-06
## ENSDARG00000081218 -11.715476 13.3276119 22.53318 2.065456e-06
## ENSDARG00000097996   6.792520  0.6288484 22.43087 2.178447e-06
## ENSDARG00000091399  -5.263067 -0.1244900 22.31618 2.312499e-06
## ENSDARG00000095618   1.929696  2.8820710 21.78112 3.055921e-06
## ENSDARG00000073820   5.316315  1.8743178 21.64067 3.288048e-06
## ENSDARG00000062788   5.233288  4.0704738 19.89996 8.160215e-06
## ENSDARG00000097157   8.596447  1.3035282 19.52685 9.919581e-06
## ENSDARG00000069503   4.461870  2.3473625 19.20727 1.172661e-05
## ENSDARG00000075193   1.754298  3.5530263 19.20687 1.172908e-05
## ENSDARG00000088811   3.062316  4.0254160 18.55743 1.648615e-05
## ENSDARG00000006427   5.600417  2.5856179 16.68330 4.416814e-05
##                             FDR
## ENSDARG00000081077 3.639448e-05
## ENSDARG00000085168 6.010851e-04
## ENSDARG00000083266 2.774348e-03
## ENSDARG00000092458 2.774348e-03
## ENSDARG00000041339 6.377582e-03
## ENSDARG00000081218 6.377582e-03
## ENSDARG00000097996 6.377582e-03
## ENSDARG00000091399 6.377582e-03
## ENSDARG00000095618 7.254421e-03
## ENSDARG00000073820 7.254421e-03
## ENSDARG00000062788 1.636717e-02
## ENSDARG00000097157 1.823798e-02
## ENSDARG00000069503 1.848419e-02
## ENSDARG00000075193 1.848419e-02
## ENSDARG00000088811 2.424893e-02
## ENSDARG00000006427 6.090511e-02

table(resultAll_DGE13filteredEE2_10 <- decideTestsDGE(All_DGE13filtered.fitEE2_10))

## 
##    -1     0     1 
##     2 22060     1

topTags(All_DGE13filtered.fitEE2_10)

## Coefficient:  treatmentfilteredEE2_10 
##                         logFC     logCPM       LR       PValue
## ENSDARG00000085168  -8.488074  6.6105761 36.56610 1.475757e-09
## ENSDARG00000081218 -13.173566 13.3276119 25.95907 3.487329e-07
## ENSDARG00000097996   6.941902  0.6288484 23.64886 1.156143e-06
## ENSDARG00000083042   3.796417  1.4065348 18.76849 1.475851e-05
## ENSDARG00000080404  -2.263103  2.3943284 17.26530 3.250696e-05
## ENSDARG00000082850  -1.652810  4.4890515 16.51254 4.832936e-05
## ENSDARG00000082628   5.649366  6.1223891 16.05600 6.149642e-05
## ENSDARG00000063276  -6.961693  0.1743958 15.95148 6.498708e-05
## ENSDARG00000083431  -2.770505  7.9692558 15.46942 8.385103e-05
## ENSDARG00000080372  -1.876902  3.8522501 15.24261 9.454575e-05
##                             FDR
## ENSDARG00000085168 3.255963e-05
## ENSDARG00000081218 3.847047e-03
## ENSDARG00000097996 8.502662e-03
## ENSDARG00000083042 8.140427e-02
## ENSDARG00000080404 1.434402e-01
## ENSDARG00000082850 1.777151e-01
## ENSDARG00000082628 1.792262e-01
## ENSDARG00000063276 1.792262e-01
## ENSDARG00000083431 2.035962e-01
## ENSDARG00000080372 2.035962e-01

All_DGE.nc.rm13filtered.fit <- glmFit(All_DGE.nc.rm13filtered, designModelfiltered)
All_DGE.nc.rm13filtered.fitEE2_3 <- glmLRT(All_DGE.nc.rm13filtered.fit, coef = "treatmentfilteredEE2_3")
All_DGE.nc.rm13filtered.fitEE2_10 <- glmLRT(All_DGE.nc.rm13filtered.fit, coef = "treatmentfilteredEE2_10")
Table.All_DGE.nc.rm13filteredEE2_3 <- topTags(All_DGE.nc.rm13filtered.fitEE2_3, n=nrow(All_DGE.nc.rm13filtered))$table
Table.All_DGE.nc.rm13filteredEE2_10 <- topTags(All_DGE.nc.rm13filtered.fitEE2_10, n=nrow(All_DGE.nc.rm13filtered))$table
table(resultAll_DGE.nc.rm13filteredEE2_3 <- decideTestsDGE(All_DGE.nc.rm13filtered.fitEE2_3))

## 
##    -1     0     1 
##     2 20799     6

topTags(All_DGE.nc.rm13filtered.fitEE2_3)

## Coefficient:  treatmentfilteredEE2_3 
##                        logFC     logCPM       LR       PValue         FDR
## ENSDARG00000092458 -4.192857  0.8989275 25.40946 4.636390e-07 0.009646936
## ENSDARG00000091399 -5.274667 -0.1238424 22.76679 1.828990e-06 0.014329419
## ENSDARG00000041339  7.393485  1.8446449 22.53262 2.066048e-06 0.014329419
## ENSDARG00000095618  1.921188  2.8800594 21.67142 3.235768e-06 0.014595547
## ENSDARG00000073820  5.304972  1.8712922 21.51683 3.507365e-06 0.014595547
## ENSDARG00000062788  5.207482  4.0533591 19.88044 8.243934e-06 0.028588588
## ENSDARG00000069503  4.444508  2.3419183 19.09502 1.243691e-05 0.036967841
## ENSDARG00000088811  3.055496  4.0237194 18.55239 1.652974e-05 0.042991793
## ENSDARG00000006427  5.574731  2.5729826 16.62898 4.545132e-05 0.105078391
## ENSDARG00000009544  5.339240  2.8029574 16.33786 5.299453e-05 0.106496225

    table(resultAll_DGE.nc.rm13filteredEE2_10 <- decideTestsDGE(All_DGE.nc.rm13filtered.fitEE2_10))

## 
##     0 
## 20807

topTags(All_DGE.nc.rm13filtered.fitEE2_10)

## Coefficient:  treatmentfilteredEE2_10 
##                        logFC     logCPM       LR       PValue       FDR
## ENSDARG00000063276 -6.972995  0.1750804 16.03918 6.204529e-05 0.9990506
## ENSDARG00000016713 -3.435634  0.2128911 15.05291 1.045389e-04 0.9990506
## ENSDARG00000088811  2.610579  4.0237194 14.10608 1.727840e-04 0.9990506
## ENSDARG00000076761 -1.188794  4.9600651 13.74334 2.095628e-04 0.9990506
## ENSDARG00000010662  1.143213  5.2883990 13.25931 2.712295e-04 0.9990506
## ENSDARG00000052564  3.237355 -0.1653769 12.46131 4.154685e-04 0.9990506
## ENSDARG00000089687  1.343595  5.7110479 12.01383 5.280724e-04 0.9990506
## ENSDARG00000003687  2.196638  1.5092463 12.00695 5.300256e-04 0.9990506
## ENSDARG00000089762 -1.387607  2.2977225 11.96382 5.424362e-04 0.9990506
## ENSDARG00000037861  1.411788  3.7600213 11.50895 6.926184e-04 0.9990506