1 Major analysis

1.1 Loading packages

library(limma)
library(edgeR)
library(magrittr)
library(janitor)
library(dplyr)
library(tibble)
library(rtracklayer)
library(venn)
library(UpSetR)

1.2 Change the working dir

Change the working directory to where the project is located.

setwd("~/Documents/Kelly_2020/Research/Six_M_Work_2020/Jo_RNAseq/Kelly_JoRNAseq_DEG_Parthway_R_codes")

1.3 Data loading

1.3.1 Input the reference

load("count_jo_Bos_taurus.ARS-UCD1.2.99.rda")
load("annotation_jo_Bos_taurus.ARS-UCD1.2.99.rda")

Hereford_ARSgtf_EST <- ebi_anno_gtf_gene_as_df

1.3.2 Input counts

Load the counts for genes from the featureCounts to identify the DEGs.

genes <- countsensembl$counts 

1.4 Make master DGEList

x <- list()

# Set the count
x$counts <- genes%>%apply(.,MARGIN = 2,FUN = as.numeric)%>%
  set_rownames(rownames(genes))

colnames(x$counts) <-  gsub(".R1.all.ARS.sorted.bam", "", colnames(x$counts ))

x$anno <- subset(Hereford_ARSgtf_EST, gene_id %in% rownames(x$counts))%>%subset(feature %in% "gene")

x$genes <- data.frame(gene_id = x$anno$gene_id,
                      chr = x$anno$chr,
                      gene_name = x$anno$gene_name,
                      gene_biotype = x$anno$gene_biotype)
#set expression level cutoff, can change from 0.5 to 1. 
sel <- rowSums(cpm(x$counts) > 1) >= 3
x$counts <- x$counts[sel,]
x$genes <- x$genes[sel,]
x$anno <- x$anno[sel,]
sampleinfo<- read.csv("sampleinfo.csv")
sampleinfo
##    rna_id Subject Phenotype  Tissue_type
## 1   546FB     546    horned    forebrain
## 2   546FS     546    horned frontal_skin
## 3   546HB     546    horned      hornbud
## 4   546MB     546    horned     midbrain
## 5   618FB     618    horned    forebrain
## 6   618FS     618    horned frontal_skin
## 7   618HB     618    horned      hornbud
## 8   618MB     618    horned     midbrain
## 9   667FB     667    polled    forebrain
## 10  667FS     667    polled frontal_skin
## 11  667HB     667    polled      hornbud
## 12  667MB     667    polled     midbrain
## 13  668FB     668    horned    forebrain
## 14  668FS     668    horned frontal_skin
## 15  668HB     668    horned      hornbud
## 16  668MB     668    horned     midbrain
## 17  698FB     698    polled    forebrain
## 18  698FS     698    polled frontal_skin
## 19  698HB     698    polled      hornbud
## 20  698MB     698    polled     midbrain
## 21  709FB     709    polled    forebrain
## 22  709FS     709    polled frontal_skin
## 23  709HB     709    polled      hornbud
## 24  709MB     709    polled     midbrain
## 25  736FB     736    horned    forebrain
## 26  736FS     736    horned frontal_skin
## 27  736HB     736    horned      hornbud
## 28  736MB     736    horned     midbrain
x$counts <- x$counts[,match(sampleinfo[,1],colnames(x$counts))]
colnames(x$counts)
##  [1] "546FB" "546FS" "546HB" "546MB" "618FB" "618FS" "618HB" "618MB" "667FB"
## [10] "667FS" "667HB" "667MB" "668FB" "668FS" "668HB" "668MB" "698FB" "698FS"
## [19] "698HB" "698MB" "709FB" "709FS" "709HB" "709MB" "736FB" "736FS" "736HB"
## [28] "736MB"
Subject <- as.factor(sampleinfo[,2])

phenotype <- as.factor(sampleinfo[,3])
tissues <- as.factor(sampleinfo[,4])
Treat <- factor(paste(phenotype,tissues,sep="."))

x <- new("DGEList", x)
dim(x)
## [1] 20297    28

1.5 Multi-level analysis

1.5.0.1 Normalization

x <- calcNormFactors(x, method="TMM")

1.5.0.2 Analysis

#set design matrix
design <- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)

xtreat <- estimateDisp(x,design)

#plot biological coefficient variation
plotBCV(xtreat)

sqrt(xtreat$common.dispersion)
## [1] 0.3368676
vtreat <- voomWithQualityWeights(xtreat,design=design,normalize.method = "none",
                            plot=T,col=as.numeric(Treat))

cols <- rep("red",28)

cols[phenotype=="horned"] <- "blue"

#mds plot with batch included
plotMDS(vtreat,label=tissues,col=cols,dim.plot=c(1,2),main="horned and polled RNA-seq MDSplot")

cols <- rep("red",28)

cols[tissues=="forebrain"] <- "blue"
cols[tissues=="frontal_skin"] <- "black"
cols[tissues=="hornbud"] <- "red"
cols[tissues=="midbrain"] <- "green"

#mds plot with batch included
plotMDS(vtreat,label=phenotype,col=cols,dim.plot=c(1,2),main="horned and polled RNA-seq MDSplot(colored by tissues)")

cols <- rep("red",28)

cols[Subject=="546"] <- "blue"
cols[Subject=="618"] <- "black"
cols[Subject=="667"] <- "red"
cols[Subject=="668"] <- "green"
cols[Subject=="698"] <- "yellow"
cols[Subject=="709"] <- "orange"
cols[Subject=="736"] <- "grey"

#mds plot with batch included
plotMDS(vtreat,label=tissues,col=cols,dim.plot=c(1,2),main="horned and polled RNA-seq MDSplot(colored by tissues)")

#### Comparisions & results

# estimate the correlation between measurements made on the same subject
corfit <- duplicateCorrelation(vtreat,design,block=sampleinfo$Subject)
corfit$consensus
## [1] 0.254882
vfittreat <- lmFit(vtreat,design,block=sampleinfo$Subject,correlation=corfit$consensus)

cm <- makeContrasts(
polledvshornedForforebrain = polled.forebrain-horned.forebrain, 
polledvshornedForfrontal_skin = polled.frontal_skin-horned.frontal_skin, 
polledvshornedForhornbud = polled.hornbud-horned.hornbud, 
polledvshornedFormidbrain = polled.midbrain-horned.midbrain, 
forebrainvsfrontal_skinForpolled = polled.forebrain-polled.frontal_skin, 
forebrainvsfrontal_skinForhorned = horned.forebrain-horned.frontal_skin,
forebrainvshornbudForpolled = polled.forebrain-polled.hornbud, 
forebrainvshornbudForhorned = horned.forebrain-horned.hornbud,
midbrainvsforebrainForpolled = polled.midbrain-polled.forebrain, 
midbrainvsforebrainForhorned = horned.midbrain-horned.forebrain,
hornbudvsfrontal_skinForpolled = polled.hornbud-polled.frontal_skin, 
hornbudvsfrontal_skinForhorned = horned.hornbud-horned.frontal_skin,
midbrainvsfrontal_skinForpolled = polled.midbrain-polled.frontal_skin, 
midbrainvsfrontal_skinForhorned = horned.midbrain-horned.frontal_skin,
midbrainvshornbudForpolled = polled.midbrain-polled.hornbud, 
midbrainvshornbudForhorned = horned.midbrain-horned.hornbud,levels=design)

vfittreat2 <- contrasts.fit(vfittreat, cm)
vfittreat2 <- eBayes(vfittreat2)

# outcome of each hypothesis test
resulttreat <- decideTests(vfittreat2,p.value=0.05, lfc = 1)

# Venn diagram showing numbers of genes significant in each comparison 
#venn(as.data.frame(resulttreat)%>%abs(), zcolor = "style")

summary(decideTests(vfittreat2,p.value=0.05,lfc = 1)) 
##        polledvshornedForforebrain polledvshornedForfrontal_skin
## Down                         1480                          1122
## NotSig                      15703                         13285
## Up                           3114                          5890
##        polledvshornedForhornbud polledvshornedFormidbrain
## Down                       1214                      1451
## NotSig                    12711                     16897
## Up                         6372                      1949
##        forebrainvsfrontal_skinForpolled forebrainvsfrontal_skinForhorned
## Down                                161                              294
## NotSig                            20014                            17607
## Up                                  122                             2396
##        forebrainvshornbudForpolled forebrainvshornbudForhorned
## Down                           114                         740
## NotSig                       20181                       15069
## Up                               2                        4488
##        midbrainvsforebrainForpolled midbrainvsforebrainForhorned
## Down                              3                           36
## NotSig                        20291                        20200
## Up                                3                           61
##        hornbudvsfrontal_skinForpolled hornbudvsfrontal_skinForhorned
## Down                                0                             43
## NotSig                          20295                          20200
## Up                                  2                             54
##        midbrainvsfrontal_skinForpolled midbrainvsfrontal_skinForhorned
## Down                               174                             302
## NotSig                           19961                           16135
## Up                                 162                            3860
##        midbrainvshornbudForpolled midbrainvshornbudForhorned
## Down                          176                        772
## NotSig                      20017                      14002
## Up                            104                       5523
1.5.0.2.1 Polled vs horned
summary(decideTests(vfittreat2,p.value=0.05,lfc = 1))[,1:4]
##        polledvshornedForforebrain polledvshornedForfrontal_skin
## Down                         1480                          1122
## NotSig                      15703                         13285
## Up                           3114                          5890
##        polledvshornedForhornbud polledvshornedFormidbrain
## Down                       1214                      1451
## NotSig                    12711                     16897
## Up                         6372                      1949
venn(as.data.frame(resulttreat)[,1:4]%>%abs(), zcolor = "style")

# transfer the data format
Upset_input <- c(
  polledvshornedForforebrain = 4594,
  polledvshornedForfrontal_skin = 7012,
  polledvshornedForhornbud = 7586,
  polledvshornedFormidbrain = 3400,
  "polledvshornedForforebrain&polledvshornedForfrontal_skin" = 3232,
  "polledvshornedForforebrain&polledvshornedForhornbud" = 3377,
  "polledvshornedForforebrain&polledvshornedFormidbrain" = 2173,
  "polledvshornedForfrontal_skin&polledvshornedForhornbud" = 6015,
  "polledvshornedForfrontal_skin&polledvshornedFormidbrain" = 2268,
  "polledvshornedForhornbud&polledvshornedFormidbrain" = 2315,
  "polledvshornedForhornbud&polledvshornedForforebrain&polledvshornedForfrontal_skin" = 2980,
  "polledvshornedForfrontal_skin&polledvshornedFormidbrain&polledvshornedForforebrain" = 1491,
  "polledvshornedForhornbud&polledvshornedFormidbrain&polledvshornedForfrontal_skin" = 2016,
  "polledvshornedForhornbud&polledvshornedFormidbrain&polledvshornedForforebrain" = 1558,
  "polledvshornedForhornbud&polledvshornedFormidbrain&polledvshornedForfrontal_skin&polledvshornedForforebrain" = 1353)


# plot
upset(fromExpression(Upset_input), 
      nintersects = NA, 
      nsets = 4, 
      mb.ratio = c(0.6, 0.4),
      number.angles = 0, 
      text.scale = 1.1, 
      point.size = 2.8, 
      line.size = 1, 
      main.bar.color = "black")

###### venn for forebrainvsfrontal_skin

venn(as.data.frame(resulttreat)[,c(5,6)]%>%abs(), zcolor = "style")

1.5.0.2.1.1 venn for forebrainvshornbud
venn(as.data.frame(resulttreat)[,c(7,8)]%>%abs(), zcolor = "style")

###### venn for midbrainvsforebrain

venn(as.data.frame(resulttreat)[,c(9,10)]%>%abs(), zcolor = "style")

###### venn for hornbudvsfrontal_skin

venn(as.data.frame(resulttreat)[,c(11,12)]%>%abs(), zcolor = "style")

###### venn for midbrainvsfrontal_skin

venn(as.data.frame(resulttreat)[,c(13,14)]%>%abs(), zcolor = "style")

###### venn for midbrainvshornbud

venn(as.data.frame(resulttreat)[,c(15,16)]%>%abs(), zcolor = "style")

1.5.0.2.2 Between tissues
1.5.0.2.2.1 polled
summary(decideTests(vfittreat2,p.value=0.05,lfc = 1))[,c(5,7,9,11,13,15)]
##        forebrainvsfrontal_skinForpolled forebrainvshornbudForpolled
## Down                                161                         114
## NotSig                            20014                       20181
## Up                                  122                           2
##        midbrainvsforebrainForpolled hornbudvsfrontal_skinForpolled
## Down                              3                              0
## NotSig                        20291                          20295
## Up                                3                              2
##        midbrainvsfrontal_skinForpolled midbrainvshornbudForpolled
## Down                               174                        176
## NotSig                           19961                      20017
## Up                                 162                        104
venn(as.data.frame(resulttreat)[,c(5,7,9,11,13,15)]%>%abs(), zcolor = "style")

# transfer the data format
Upset_input <- c(
  forebrainvsfrontal_skinForpolled = 283,
  forebrainvshornbudForpolled = 116,
  midbrainvsforebrainForpolled = 6,
  hornbudvsfrontal_skinForpolled = 2,
  midbrainvsfrontal_skinForpolled = 336, 
  midbrainvshornbudForpolled = 280,
  "forebrainvsfrontal_skinForpolled&forebrainvshornbudForpolled" = 68,
  "forebrainvsfrontal_skinForpolled&midbrainvsforebrainForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&hornbudvsfrontal_skinForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled" = 106,
  "forebrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 58,
  "forebrainvshornbudForpolled&midbrainvsforebrainForpolled" = 3,
  "forebrainvshornbudForpolled&hornbudvsfrontal_skinForpolled" = 0,
  "forebrainvshornbudForpolled&midbrainvsfrontal_skinForpolled" = 41,
  "forebrainvshornbudForpolled&midbrainvshornbudForpolled" = 56,
  "midbrainvsforebrainForpolled&hornbudvsfrontal_skinForpolled" = 0,
  "midbrainvsforebrainForpolled&midbrainvsfrontal_skinForpolled" = 6,
  "midbrainvsforebrainForpolled&midbrainvshornbudForpolled" = 3,
  "hornbudvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled" = 1,
  "hornbudvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 0,
  "midbrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 120,
  "forebrainvsfrontal_skinForpolled&forebrainvshornbudForpolled&midbrainvsforebrainForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&forebrainvshornbudForpolled&hornbudvsfrontal_skinForpolled" = 0,
  "forebrainvsfrontal_skinForpolled&forebrainvshornbudForpolled&midbrainvsfrontal_skinForpolled" = 34,
  "forebrainvsfrontal_skinForpolled&forebrainvshornbudForpolled&midbrainvshornbudForpolled" = 39,
  "forebrainvsfrontal_skinForpolled&midbrainvsforebrainForpolled&hornbudvsfrontal_skinForpolled" = 0,
  "forebrainvsfrontal_skinForpolled&midbrainvsforebrainForpolled&midbrainvsfrontal_skinForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&midbrainvsforebrainForpolled&midbrainvshornbudForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&hornbudvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled" = 1,
  "forebrainvsfrontal_skinForpolled&hornbudvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 0,
  "forebrainvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 50,
  "forebrainvshornbudForpolled&midbrainvsforebrainForpolled&hornbudvsfrontal_skinForpolled" = 0,
  "forebrainvshornbudForpolled&midbrainvsforebrainForpolled&midbrainvsfrontal_skinForpolled" = 3,
  "forebrainvshornbudForpolled&midbrainvsforebrainForpolled&midbrainvshornbudForpolled" = 2,
  "forebrainvshornbudForpolled&hornbudvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled" = 0,
  "forebrainvshornbudForpolled&hornbudvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 0,
  "forebrainvshornbudForpolled&midbrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" =36,
  "midbrainvsforebrainForpolled&hornbudvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled" = 0,
  "midbrainvsforebrainForpolled&hornbudvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 0,
  "midbrainvsforebrainForpolled&midbrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 3,
  "hornbudvsfrontal_skinForpolled&midbrainvsfrontal_skinForpolled&midbrainvshornbudForpolled" = 0)

# plot
upset(fromExpression(Upset_input), 
      nintersects = NA, 
      nsets = 6, 
      mb.ratio = c(0.6, 0.4),
      number.angles = 0, 
      text.scale = 1.1, 
      point.size = 2.8, 
      line.size = 1, 
      main.bar.color = "red")

1.5.0.2.2.2 horned
summary(decideTests(vfittreat2,p.value=0.05,lfc = 1))[,c(6,8,10,12,14,16)]
##        forebrainvsfrontal_skinForhorned forebrainvshornbudForhorned
## Down                                294                         740
## NotSig                            17607                       15069
## Up                                 2396                        4488
##        midbrainvsforebrainForhorned hornbudvsfrontal_skinForhorned
## Down                             36                             43
## NotSig                        20200                          20200
## Up                               61                             54
##        midbrainvsfrontal_skinForhorned midbrainvshornbudForhorned
## Down                               302                        772
## NotSig                           16135                      14002
## Up                                3860                       5523
venn(as.data.frame(resulttreat)[,c(6,8,10,12,14,16)]%>%abs(), zcolor = "style")

# transfer the data format
Upset_input <- c(
  forebrainvsfrontal_skinForhorned = 2690,
  forebrainvshornbudForhorned = 5228,
  midbrainvsforebrainForhorned = 97,
  hornbudvsfrontal_skinForhorned = 97,
  midbrainvsfrontal_skinForhorned = 4162, 
  midbrainvshornbudForhorned = 6295,
  "forebrainvsfrontal_skinForhorned&forebrainvshornbudForhorned" = 2292,
  "forebrainvsfrontal_skinForhorned&midbrainvsforebrainForhorned" = 18,
  "forebrainvsfrontal_skinForhorned&hornbudvsfrontal_skinForhorned" = 18,
  "forebrainvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned" = 2381,
  "forebrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 2359,
  "forebrainvshornbudForhorned&midbrainvsforebrainForhorned" = 33,
  "forebrainvshornbudForhorned&hornbudvsfrontal_skinForhorned" = 60,
  "forebrainvshornbudForhorned&midbrainvsfrontal_skinForhorned" = 3100,
  "forebrainvshornbudForhorned&midbrainvshornbudForhorned" = 4737,
  "midbrainvsforebrainForhorned&hornbudvsfrontal_skinForhorned" = 1,
  "midbrainvsforebrainForhorned&midbrainvsfrontal_skinForhorned" = 64,
  "midbrainvsforebrainForhorned&midbrainvshornbudForhorned" = 66,
  "hornbudvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned" = 20,
  "hornbudvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 58,
  "midbrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 3800,
  
  "forebrainvsfrontal_skinForhorned&forebrainvshornbudForhorned&midbrainvsforebrainForhorned" = 13,
  "forebrainvsfrontal_skinForhorned&forebrainvshornbudForhorned&hornbudvsfrontal_skinForhorned" = 4,
  "forebrainvsfrontal_skinForhorned&forebrainvshornbudForhorned&midbrainvsfrontal_skinForhorned" = 2075,
  "forebrainvsfrontal_skinForhorned&forebrainvshornbudForhorned&midbrainvshornbudForhorned" = 2187,
  "forebrainvsfrontal_skinForhorned&midbrainvsforebrainForhorned&hornbudvsfrontal_skinForhorned" = 0,
  "forebrainvsfrontal_skinForhorned&midbrainvsforebrainForhorned&midbrainvsfrontal_skinForhorned" = 6,
  "forebrainvsfrontal_skinForhorned&midbrainvsforebrainForhorned&midbrainvshornbudForhorned" = 9,
  "forebrainvsfrontal_skinForhorned&hornbudvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned" = 13,
  "forebrainvsfrontal_skinForhorned&hornbudvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 4,
  "forebrainvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 2221,
  "forebrainvshornbudForhorned&midbrainvsforebrainForhorned&hornbudvsfrontal_skinForhorned" = 1,
  "forebrainvshornbudForhorned&midbrainvsforebrainForhorned&midbrainvsfrontal_skinForhorned" = 16,
  "forebrainvshornbudForhorned&midbrainvsforebrainForhorned&midbrainvshornbudForhorned" = 14,
  "forebrainvshornbudForhorned&hornbudvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned" = 6,
  "forebrainvshornbudForhorned&hornbudvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 49,
  "forebrainvshornbudForhorned&midbrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" =3068,
  "midbrainvsforebrainForhorned&hornbudvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned" = 0,
  "midbrainvsforebrainForhorned&hornbudvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 0,
  "midbrainvsforebrainForhorned&midbrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 58,
  "hornbudvsfrontal_skinForhorned&midbrainvsfrontal_skinForhorned&midbrainvshornbudForhorned" = 2)

# plot
upset(fromExpression(Upset_input), 
      nintersects = NA, 
      nsets = 6, 
      mb.ratio = c(0.6, 0.4),
      number.angles = 0, 
      text.scale = 1.1, 
      point.size = 2.8, 
      line.size = 1, 
      main.bar.color = "blue")

1.5.0.4 Save results in CSV files

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForforebrain",p.value=0.05,lfc = 1)%>%remove_rownames(),file="polledvshornedForforebrain-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForfrontal_skin",p.value=0.05,lfc = 1)%>%remove_rownames(),file="polledvshornedForfrontal_skin-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForhornbud",p.value=0.05,lfc = 1)%>%remove_rownames(),file="polledvshornedForhornbud-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedFormidbrain",p.value=0.05,lfc = 1)%>%remove_rownames(),file="polledvshornedFormidbrain-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvsfrontal_skinForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="forebrainvsfrontal_skinForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvsfrontal_skinForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="forebrainvsfrontal_skinForhorned-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvshornbudForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="forebrainvshornbudForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvshornbudForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="forebrainvshornbudForhorned-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsforebrainForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvsforebrainForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsforebrainForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvsforebrainForhorned-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="hornbudvsfrontal_skinForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="hornbudvsfrontal_skinForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="hornbudvsfrontal_skinForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="hornbudvsfrontal_skinForhorned-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsfrontal_skinForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvsfrontal_skinForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsfrontal_skinForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvsfrontal_skinForhorned-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvshornbudForpolled",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvshornbudForpolled-ARS_jo-DEG_p0.05_lfc1.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvshornbudForhorned",p.value=0.05,lfc = 1)%>%remove_rownames(),file="midbrainvshornbudForhorned-ARS_jo-DEG_p0.05_lfc1.csv")
write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForforebrain")%>%remove_rownames(),file="polledvshornedForforebrain-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForfrontal_skin")%>%remove_rownames(),file="polledvshornedForfrontal_skin-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedForhornbud")%>%remove_rownames(),file="polledvshornedForhornbud-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="polledvshornedFormidbrain")%>%remove_rownames(),file="polledvshornedFormidbrain-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvsfrontal_skinForpolled")%>%remove_rownames(),file="forebrainvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvsfrontal_skinForhorned")%>%remove_rownames(),file="forebrainvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvshornbudForpolled")%>%remove_rownames(),file="forebrainvshornbudForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="forebrainvshornbudForhorned")%>%remove_rownames(),file="forebrainvshornbudForhorned-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsforebrainForpolled")%>%remove_rownames(),file="midbrainvsforebrainForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsforebrainForhorned")%>%remove_rownames(),file="midbrainvsforebrainForhorned-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="hornbudvsfrontal_skinForpolled")%>%remove_rownames(),file="hornbudvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="hornbudvsfrontal_skinForhorned")%>%remove_rownames(),file="hornbudvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsfrontal_skinForpolled")%>%remove_rownames(),file="midbrainvsfrontal_skinForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvsfrontal_skinForhorned")%>%remove_rownames(),file="midbrainvsfrontal_skinForhorned-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvshornbudForpolled")%>%remove_rownames(),file="midbrainvshornbudForpolled-ARS_jo-DEG_all.csv")

write.csv(topTable(vfittreat2,n=Inf,coef="midbrainvshornbudForhorned")%>%remove_rownames(),file="midbrainvshornbudForhorned-ARS_jo-DEG_all.csv")

2 Appendix

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## 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  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] UpSetR_1.4.0         venn_1.9             rtracklayer_1.48.0  
##  [4] GenomicRanges_1.40.0 GenomeInfoDb_1.24.2  IRanges_2.22.2      
##  [7] S4Vectors_0.26.1     BiocGenerics_0.34.0  tibble_3.0.3        
## [10] dplyr_1.0.2          janitor_2.0.1        magrittr_1.5        
## [13] edgeR_3.30.3         limma_3.44.3        
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5                  locfit_1.5-9.4             
##  [3] lubridate_1.7.9             lattice_0.20-41            
##  [5] Rsamtools_2.4.0             Biostrings_2.56.0          
##  [7] digest_0.6.25               R6_2.4.1                   
##  [9] plyr_1.8.6                  evaluate_0.14              
## [11] ggplot2_3.3.2               pillar_1.4.6               
## [13] zlibbioc_1.34.0             rlang_0.4.7                
## [15] Matrix_1.2-18               rmarkdown_2.3              
## [17] labeling_0.3                splines_4.0.2              
## [19] statmod_1.4.34              BiocParallel_1.22.0        
## [21] stringr_1.4.0               RCurl_1.98-1.2             
## [23] munsell_0.5.0               DelayedArray_0.14.1        
## [25] compiler_4.0.2              xfun_0.16                  
## [27] pkgconfig_2.0.3             htmltools_0.5.0            
## [29] tidyselect_1.1.0            SummarizedExperiment_1.18.2
## [31] gridExtra_2.3               GenomeInfoDbData_1.2.3     
## [33] bookdown_0.20               matrixStats_0.56.0         
## [35] XML_3.99-0.5                crayon_1.3.4               
## [37] GenomicAlignments_1.24.0    bitops_1.0-6               
## [39] grid_4.0.2                  gtable_0.3.0               
## [41] lifecycle_0.2.0             scales_1.1.1               
## [43] stringi_1.4.6               farver_2.0.3               
## [45] XVector_0.28.0              snakecase_0.11.0           
## [47] ellipsis_0.3.1              admisc_0.8                 
## [49] generics_0.0.2              vctrs_0.3.2                
## [51] tools_4.0.2                 Biobase_2.48.0             
## [53] glue_1.4.1                  purrr_0.3.4                
## [55] yaml_2.2.1                  colorspace_1.4-1           
## [57] knitr_1.29