Histone Acetyl RNA-seq DGE

Source: https://github.com/markziemann/histone_acetyl_rnaseq

Introduction

In this study we are looking at the effect of knock-down of two genes which are required for the supply of acetyl units. It is hypothesised that knock-down of these genes could lead to problems with the response to exercise. We have both ctrl and kd under either sedentary or exercise conditions.

The focus of this report is to look at the overall structure of the dataset, understand whether the knock-down was successful and identify any outliers.

Mark’s notes: the ctrl vs kd comparisons to be conducted as paired analysis. Need to include new contrasts for Sed vs Ex (unpaired). Using the unpaired analysis will reduce statistical power by a bit, but we’ll see.

suppressPackageStartupMessages({
    library("zoo")
    library("tidyverse")
    library("reshape2")
    library("DESeq2")
    library("gplots")
    library("fgsea")
    library("MASS")
    library("mitch")
    library("eulerr")
    library("limma")
    library("topconfects")
    library("kableExtra")
    library("vioplot")
    library("beeswarm")
})

Import read counts

Importing RNA-seq data

tmp <- read.table("3col.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
xx$geneid = NULL
xx <- round(xx)
xx <- xx[,which(colnames(xx)!="test")]
xx[1:6,1:6]

##                             5176140_S26_L001 5176140_S26_L002 5176141_S27_L001
## ENSMUSG00000000001.5 Gnai3                23               17               29
## ENSMUSG00000000003.16 Pbsn                 0                0                0
## ENSMUSG00000000028.16 Cdc45                1                1                2
## ENSMUSG00000000031.17 H19                441              448              281
## ENSMUSG00000000037.18 Scml2                0                0                0
## ENSMUSG00000000049.12 Apoh                 0                1                3
##                             5176141_S27_L002 5176142_S28_L001 5176142_S28_L002
## ENSMUSG00000000001.5 Gnai3                11               17               13
## ENSMUSG00000000003.16 Pbsn                 0                0                0
## ENSMUSG00000000028.16 Cdc45                1                1                1
## ENSMUSG00000000031.17 H19                289              355              384
## ENSMUSG00000000037.18 Scml2                0                0                0
## ENSMUSG00000000049.12 Apoh                 0                1                0

dim(xx)

## [1] 55357   112

Fix the sample names.

They are duplicated for lane 1 and 2, which I will aggregate.

labels <- unique(sapply(strsplit(colnames(xx),"_"),"[[",1))
l <- lapply(labels,function(x) { rowSums(xx[,grep(x,colnames(xx))]) } )
ll <- do.call(cbind,l)
colnames(ll) <- labels
ll <- as.data.frame(ll[,order(colnames(ll))])
write.table(ll,file="counts.tsv",sep="\t",quote=FALSE)
rpm <- apply(ll,2, function(x) { x / sum(x) * 1000000 } )
write.table(rpm,file="rpm.tsv",sep="\t",quote=FALSE)

head(ll)

##                             5176140 5176141 5176142 5176143 5176144 5176145
## ENSMUSG00000000001.5 Gnai3       40      40      30      43      37      26
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       2       3       2       0       1       1
## ENSMUSG00000000031.17 H19       889     570     739     686     536     639
## ENSMUSG00000000037.18 Scml2       0       0       0       0       0       0
## ENSMUSG00000000049.12 Apoh        1       3       1       0       0       0
##                             5176146 5176147 5176148 5176149 5176150 5176151
## ENSMUSG00000000001.5 Gnai3       28      40      45      40      37      60
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       3       1       2       5       0       6
## ENSMUSG00000000031.17 H19       554     901     441     614     945     410
## ENSMUSG00000000037.18 Scml2       0       0       1       1       0       1
## ENSMUSG00000000049.12 Apoh        0       0       0       1       0       4
##                             5176152 5176153 5176154 5176155 5176156 5176157
## ENSMUSG00000000001.5 Gnai3       56      28      19      47      29      27
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       6       7       3       7       4       7
## ENSMUSG00000000031.17 H19       235     479     736     284     377     580
## ENSMUSG00000000037.18 Scml2       1       1       0       0       0       0
## ENSMUSG00000000049.12 Apoh        0       0       0       0       0       0
##                             5176158 5176159 5176160 5176161 5176162 5176163
## ENSMUSG00000000001.5 Gnai3       39      30      34      45      37      24
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       6       3       2       2       3       3
## ENSMUSG00000000031.17 H19       628     754     327     612     656     883
## ENSMUSG00000000037.18 Scml2       1       1       2       0       0       0
## ENSMUSG00000000049.12 Apoh        0       0       0       0       0       0
##                             5176164 5176165 5176166 5176167 5176168 5176169
## ENSMUSG00000000001.5 Gnai3       29      31      58      39      55      63
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       1       7       2       3      12       5
## ENSMUSG00000000031.17 H19       674     839    1451    1102     889     707
## ENSMUSG00000000037.18 Scml2       0       1       0       0       0       0
## ENSMUSG00000000049.12 Apoh        2       1       1       0       1       0
##                             5176170 5176171 5176172 5176173 5176174 5176175
## ENSMUSG00000000001.5 Gnai3       37      38      32      39      30      36
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       3       3       4       2       1       4
## ENSMUSG00000000031.17 H19       899     937     968     512     850     426
## ENSMUSG00000000037.18 Scml2       0       0       0       0       1       0
## ENSMUSG00000000049.12 Apoh        0       0       0       1       2       0
##                             5176176 5176177 5176178 5176179 5176180 5176181
## ENSMUSG00000000001.5 Gnai3       42      39      48      40      67      22
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       5       4       1       3      13       2
## ENSMUSG00000000031.17 H19       998    1025     879     644     273     544
## ENSMUSG00000000037.18 Scml2       2       2       0       0       0       0
## ENSMUSG00000000049.12 Apoh        0       0       0       0       0       0
##                             5176182 5176183 5176184 5176185 5176186 5176187
## ENSMUSG00000000001.5 Gnai3       37      31      30      45      38      49
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       3       1       3       4       7       7
## ENSMUSG00000000031.17 H19       617     642     977     696     881    1169
## ENSMUSG00000000037.18 Scml2       0       1       1       0       0       0
## ENSMUSG00000000049.12 Apoh        0       2       2       0       1       0
##                             5176188 5176189 5176190 5176191 5176192 5176193
## ENSMUSG00000000001.5 Gnai3       54      47      41      30      34      56
## ENSMUSG00000000003.16 Pbsn        0       0       0       0       0       0
## ENSMUSG00000000028.16 Cdc45       8       5       5       4       2       5
## ENSMUSG00000000031.17 H19       642     581     646     793     697     828
## ENSMUSG00000000037.18 Scml2       0       0       0       0       0       0
## ENSMUSG00000000049.12 Apoh        0       0       1       0       1       2
##                             5176194 5176195
## ENSMUSG00000000001.5 Gnai3       33      42
## ENSMUSG00000000003.16 Pbsn        0       0
## ENSMUSG00000000028.16 Cdc45       5       5
## ENSMUSG00000000031.17 H19       988     892
## ENSMUSG00000000037.18 Scml2       0       0
## ENSMUSG00000000049.12 Apoh        0       0

dim(xx)

## [1] 55357   112

dim(ll)

## [1] 55357    56

xx <- ll

Make a sample sheet.

ss <- read.table("samplesheet.tsv")

# fix the capitalisation
ss$ex_sed <- gsub("SED","Sed",ss$ex_sed)

colnames(xx) == ss$Sample_ID

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

rownames(ss) <- ss$Sample_ID

ss$animal_id <- gsub("Right","",gsub("Left","",ss$Orginal_ID))

ss %>% kbl(caption = "Sample sheet") %>%
  kable_paper("hover", full_width = F)

Sample sheet

	Sample_ID	Orginal_ID	percent	num_reads	Gbases	target_gene	ex_sed	construct	animal_id
5176140	5176140	1Left	0.3412	5860429	1.465107	ATP-CL	Sed	kd	1
5176141	5176141	2Left	0.2909	4996480	1.249120	ATP-CL	Sed	kd	2
5176142	5176142	3Left	0.3103	5329693	1.332423	ATP-CL	Sed	kd	3
5176143	5176143	5Left	0.3037	5216332	1.304083	ACSS2	Sed	kd	5
5176144	5176144	6Left	0.2452	4211539	1.052885	ACSS2	Sed	kd	6
5176145	5176145	8Left	0.3060	5255836	1.313959	ACSS2	Ex	kd	8
5176146	5176146	11Left	0.2811	4828155	1.207039	ATP-CL	Ex	kd	11
5176147	5176147	13Left	0.3101	5326258	1.331564	ACSS2	Ex	kd	13
5176148	5176148	17Left	0.3048	5235225	1.308806	ATP-CL	Sed	kd	17
5176149	5176149	18Left	0.2875	4938081	1.234520	ATP-CL	Ex	kd	18
5176150	5176150	23Left	0.3212	5516910	1.379228	ACSS2	Sed	kd	23
5176151	5176151	21Left	0.2896	4974151	1.243538	ACSS2	Ex	kd	21
5176152	5176152	31Left	0.2603	4470896	1.117724	ACSS2	Ex	kd	31
5176153	5176153	36Left	0.3191	5480841	1.370210	ATP-CL	Ex	kd	36
5176154	5176154	38Left	0.2763	4745711	1.186428	ACSS2	Ex	kd	38
5176155	5176155	42Left	0.2947	5061748	1.265437	ATP-CL	Sed	kd	42
5176156	5176156	43Left	0.3414	5863864	1.465966	ATP-CL	Ex	kd	43
5176157	5176157	44Left	0.3028	5200873	1.300218	ATP-CL	Ex	kd	44
5176158	5176158	45Left	0.2799	4807544	1.201886	ACSS2	Ex	kd	45
5176159	5176159	46Left	0.3113	5346869	1.336717	ACSS2	Sed	kd	46
5176160	5176160	51Left	0.2578	4427956	1.106989	ATP-CL	Sed	kd	51
5176161	5176161	55Left	0.3013	5175109	1.293777	ACSS2	Sed	kd	55
5176162	5176162	57Left	0.3502	6015013	1.503753	ATP-CL	Sed	kd	57
5176163	5176163	59Left	0.3226	5540957	1.385239	ATP-CL	Ex	kd	59
5176164	5176164	60Left	0.3099	5322822	1.330706	ATP-CL	Ex	kd	60
5176165	5176165	61Left	0.3006	5163086	1.290772	ACSS2	Sed	kd	61
5176166	5176166	63Left	0.3385	5814054	1.453514	ACSS2	Ex	kd	63
5176167	5176167	64Left	0.3145	5401832	1.350458	ACSS2	Sed	kd	64
5176168	5176168	1Right	0.3798	6523420	1.630855	ATP-CL	Sed	ctrl	1
5176169	5176169	2Right	0.3058	5252401	1.313100	ATP-CL	Sed	ctrl	2
5176170	5176170	3Right	0.3551	6099175	1.524794	ATP-CL	Sed	ctrl	3
5176171	5176171	5Right	0.3466	5953179	1.488295	ACSS2	Sed	ctrl	5
5176172	5176172	6Right	0.3413	5862147	1.465537	ACSS2	Sed	ctrl	6
5176173	5176173	8Right	0.3493	5999554	1.499889	ACSS2	Ex	ctrl	8
5176174	5176174	11Right	0.2878	4943234	1.235809	ATP-CL	Ex	ctrl	11
5176175	5176175	13Right	0.2869	4927776	1.231944	ACSS2	Ex	ctrl	13
5176176	5176176	17Right	0.4182	7182976	1.795744	ATP-CL	Sed	ctrl	17
5176177	5176177	18Right	0.3428	5887911	1.471978	ATP-CL	Ex	ctrl	18
5176178	5176178	23Right	0.3309	5683517	1.420879	ACSS2	Sed	ctrl	23
5176179	5176179	21Right	0.3609	6198795	1.549699	ACSS2	Ex	ctrl	21
5176180	5176180	31Right	0.3253	5587332	1.396833	ACSS2	Ex	ctrl	31
5176181	5176181	36Right	0.2955	5075489	1.268872	ATP-CL	Ex	ctrl	36
5176182	5176182	38Right	0.3686	6331050	1.582762	ACSS2	Ex	ctrl	38
5176183	5176183	42Right	0.2869	4927776	1.231944	ATP-CL	Sed	ctrl	42
5176184	5176184	43Right	0.3843	6600712	1.650178	ATP-CL	Ex	ctrl	43
5176185	5176185	44Right	0.3362	5774550	1.443637	ATP-CL	Ex	ctrl	44
5176186	5176186	45Right	0.3223	5535804	1.383951	ACSS2	Ex	ctrl	45
5176187	5176187	46Right	0.3214	5520346	1.380086	ACSS2	Sed	ctrl	46
5176188	5176188	51Right	0.3238	5561568	1.390392	ATP-CL	Sed	ctrl	51
5176189	5176189	55Right	0.3069	5271295	1.317824	ACSS2	Sed	ctrl	55
5176190	5176190	57Right	0.3376	5798596	1.449649	ATP-CL	Sed	ctrl	57
5176191	5176191	59Right	0.3315	5693823	1.423456	ATP-CL	Ex	ctrl	59
5176192	5176192	60Right	0.2904	4987892	1.246973	ATP-CL	Ex	ctrl	60
5176193	5176193	61Right	0.3326	5712716	1.428179	ACSS2	Sed	ctrl	61
5176194	5176194	63Right	0.3239	5563286	1.390821	ACSS2	Ex	ctrl	63
5176195	5176195	64Right	0.2913	5003350	1.250838	ACSS2	Sed	ctrl	64

QC analysis

Here I’ll look at a few different quality control measures.

par(mar=c(5,8,3,1))
barplot(colSums(ll),horiz=TRUE,las=1,xlab="num reads",col=ss$cols)

sums <- colSums(ll)
sums <- sums[order(sums)]
barplot(sums,horiz=TRUE,las=1,xlab="num reads")
abline(v=15000000,col="red")

MDS plot for all samples

Multidimensional scaling plot to show the variation between all samples, very similar to PCA.

Firstly with the data before aggregating technical replicates.

mds <- cmdscale(dist(t(xx)))

plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=19, cex=4 ,col="gray")
text(mds, labels=rownames(mds) )

That doesn’t make much sense.

Now I will colour the samples by group.

pch <- ss$ex_sed
pch <- gsub("Sed","15",pch)
pch <- gsub("Ex","17",pch)
pch <- as.numeric(pch)

col <- paste(ss$target_gene,ss$construct)
col <- gsub("ATP-CL ctrl","pink",col)
col <- gsub("ATP-CL kd","red",col)
col <- gsub("ACSS2 ctrl","lightblue",col)
col <- gsub("ACSS2 kd","blue",col)

plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=pch, cex=4 ,col=col)

legend("bottomleft", legend=c("ATP-CL ctrl", "ATP-CL kd", "ACSS2 ctrl", "ACSS2 kd"),
       pch=19,col=c("pink", "red", "lightblue", "blue"), cex=1)

legend("bottomright", legend=c("Sed", "Ex"),
       pch=c(15,17), cex=1)

Now I’ll separate ACSS2 and ATP-CL experiments into different charts

lightblue = ctrl

pink = kd

circle = sedentary

square = exercise

I will also make a barplot and boxplot of the ATP-CL gene expression, encoded by Acly. It will show whether there has been a measurable decrease in Acly expression between ctrl and kd groups.

ss_atp <- ss[which(ss$target_gene=="ATP-CL"),]
xx_atp <- xx[,which(colnames(xx) %in% rownames(ss_atp))]

rpm_atp <- apply(xx_atp,2, function(x) { x / sum(x) * 1000000 } )

acly <- unlist(rpm_atp[grep("Acly",rownames(rpm_atp)),,drop=TRUE])
names(acly) <- paste(names(acly), ss_atp$construct)

barplot(acly,horiz=TRUE,las=1,main="Acly expression",xlab="reads per million")

ctrl <- acly[grep("ctrl",names(acly))]
kd <- acly[grep("kd",names(acly))]
mylist <- list("ctrl"=ctrl,"kd"=kd)
boxplot(mylist,col="white",ylab="reads per million",main="Acly expression")
beeswarm(mylist,pch=19,add=TRUE)
myp <- signif(t.test(kd,ctrl)$p.value,3)
mtext(paste("p=",myp))

ratio <- kd/ctrl
ratio <- ratio[order(ratio)]
barplot(ratio,horiz=TRUE,las=1,xlab="relative quantification (kd/ctrl)",main="ATP-CL")

# shapes for exercise or sedentary
shapes <- as.numeric(factor(ss_atp$ex_sed))+14
# colours for knock out or control
colours <- as.numeric(factor(ss_atp$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_atp)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=shapes, cex=4 ,col=colours)
mtext("ATP-CL study")
text(mds, labels=ss_atp$animal_id , cex=1)

legend("topright",
  legend = c("ctrl sed", "ctrl ex", "kd sed" , "kd ex"),
  col = c("lightblue","lightblue","pink","pink"),
  pch = c(16,15,16,15),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Look at the effect of the ATP-CL KD in the Sed animals.

ss_atp_sed <- ss_atp[which(ss_atp$ex_sed == "Sed"),]

xx_atp_sed <- xx_atp[,which(colnames(xx_atp) %in% rownames(ss_atp_sed))]

colours <- as.numeric(factor(ss_atp_sed$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_atp_sed)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=16, cex=4 ,col=colours)
mtext("Effect of ATP-CL KD in Sed animals")
text(mds, labels=ss_atp_sed$animal_id , cex=1)

legend("topright",
  legend = c("ctrl Sed", "kd Sed"),
  col = c("lightblue","pink"),
  pch = c(16,16),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Look at the effect of the ATP-CL KD in the Ex animals.

ss_atp_ex <- ss_atp[which(ss_atp$ex_sed == "Ex"),]

xx_atp_ex <- xx_atp[,which(colnames(xx_atp) %in% rownames(ss_atp_ex))]

colours <- as.numeric(factor(ss_atp_ex$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_atp_sed)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=15, cex=4 ,col=colours)
mtext("Effect of ATP-CL KD in Ex animals")
text(mds, labels=ss_atp_ex$animal_id , cex=1)

legend("topright",
  legend = c("ctrl Ex", "kd Ex"),
  col = c("lightblue","pink"),
  pch = c(15,15),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Now look at the ACSS2 study.

I will also make a barplot of the ACSS2 gene expression, encoded by Acss2. It will show whether there has been a measurable decrease in Acss2 expression between ctrl and kd groups.

ss_acs <- ss[which(ss$target_gene=="ACSS2"),]
xx_acs <- xx[,which(colnames(xx) %in% rownames(ss_acs))]

rpm_acs <- apply(xx_acs,2, function(x) { x / sum(x) * 1000000 } )

acss <- unlist(rpm_acs[grep("Acss2$",rownames(rpm_acs)),,drop=TRUE])
names(acss) <- paste(names(acss) , ss_acs$construct)

barplot(acss,horiz=TRUE,las=1,main="Acss2 expression",xlab="reads per million")

ctrl <- acss[grep("ctrl",names(acss))]
kd <- acss[grep("kd",names(acss))]
mylist <- list("ctrl"=ctrl,"kd"=kd)
boxplot(mylist,col="white",ylab="reads per million",main="Acss2 expression")
beeswarm(mylist,pch=19,add=TRUE)
myp <- signif(t.test(kd,ctrl)$p.value,3)
mtext(paste("p=",myp))

ratio <- kd/ctrl
ratio <- ratio[order(ratio)]
barplot(ratio,horiz=TRUE,las=1,xlab="relative quantification (kd/ctrl)",main="Acss2")

# shapes for exercise or sedentary
shapes <- as.numeric(factor(ss_acs$ex_sed))+14
# colours for knock out or control
colours <- as.numeric(factor(ss_acs$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_acs)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=shapes, cex=4 ,col=colours)
mtext("ACSS2 study")
text(mds, labels=ss_acs$animal_id , cex=1)

legend("topright",
  legend = c("ctrl sed", "ctrl ex", "kd sed" , "kd ex"),
  col = c("lightblue","lightblue","pink","pink"),
  pch = c(16,15,16,15),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Look at the effect of the ACSS2 KD in the Sed animals.

ss_acs_sed <- ss_acs[which(ss_acs$ex_sed == "Sed"),]

xx_acs_sed <- xx_acs[,which(colnames(xx_acs) %in% rownames(ss_acs_sed))]

colours <- as.numeric(factor(ss_acs_sed$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_acs_sed)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=16, cex=4 ,col=colours)
mtext("Effect of ACSS2 KD in Sed animals")
text(mds, labels=ss_acs_sed$animal_id , cex=1)

legend("topright",
  legend = c("ctrl Sed", "kd Sed"),
  col = c("lightblue","pink"),
  pch = c(16,16),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Look at the effect of the ACSS2 KD in the Ex animals.

ss_acs_ex <- ss_acs[which(ss_acs$ex_sed == "Ex"),]

xx_acs_ex <- xx_acs[,which(colnames(xx_acs) %in% rownames(ss_acs_ex))]

colours <- as.numeric(factor(ss_acs_ex$construct))
colours <- gsub("2","pink",gsub("1","lightblue",as.character(colours)))

mds <- cmdscale(dist(t(xx_acs_sed)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=15, cex=4 ,col=colours)
mtext("Effect of ACSS2 KD in Ex animals")
text(mds, labels=ss_acs_ex$animal_id , cex=1)

legend("topright",
  legend = c("ctrl Ex", "kd Ex"),
  col = c("lightblue","pink"),
  pch = c(15,15),
  pt.cex = 2,
  cex = 1.2,
  text.col = "black",
  horiz = F ,
  inset = c(0.1, 0.1))

Correlation heatmap

Search for outliers with a correlation heatmap.

5176172 is potentially an outlier.

heatmap.2(cor(xx),trace="n",main="Pearson correlation heatmap",margin=c(8,8))

heatmap.2(cor(xx_atp),trace="n",main="Cor - ATP-CL",margin=c(8,8))

heatmap.2(cor(xx_acs),trace="n",main="Cor - ACSS2",margin=c(8,8))

Set up the different datasets for differential expression analysis

Don’t forget to remove poorly detected genes from the matrix with a threshold of 10 reads per sample on average.

There are 4 contrasts to set up.

1. Effect of ATP-CL KD in Sed animals. (ss_atp_sed)

2. Effect of ATP-CL KD in Ex animals. (ss_atp_ex)

3. Effect of ACSS2 KD in Sed animals. (ss_acs_sed)

4. Effect of ACSS2 KD in Ex animals. (ss_acs_ex)

There may be other contrasts we will add in future.

DGE 1 Effect of ATP-CL KD in Sed animals

dim(xx_atp_sed)

## [1] 55357    14

xx_atp_sed <- xx_atp_sed[which(rowMeans(xx_atp_sed)>=10),]
dim(xx_atp_sed)

## [1] 9731   14

ss_atp_sed$construct <- factor(ss_atp_sed$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_atp_sed , colData = ss_atp_sed,
  design = ~ construct )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000052305.7 Hbb-bs	1991.61520	2.7445520	0.6585414	4.167623	0.0000308	0.1314082
ENSMUSG00000073940.4 Hbb-bt	298.94874	2.7360817	0.6619947	4.133087	0.0000358	0.1314082
ENSMUSG00000069917.8 Hba-a2	937.22103	2.8861833	0.7031687	4.104539	0.0000405	0.1314082
ENSMUSG00000069919.8 Hba-a1	895.64500	2.8292112	0.7071852	4.000665	0.0000632	0.1536638
ENSMUSG00000025270.14 Alas2	32.23192	2.8753901	0.7345283	3.914608	0.0000906	0.1762311
ENSMUSG00000032754.15 Slc8b1	15.52746	-0.8772522	0.2553896	-3.434956	0.0005926	0.9611794
ENSMUSG00000035443.11 Thyn1	30.79490	0.5978610	0.1879002	3.181801	0.0014636	0.9998239
ENSMUSG00000031448.13 Adprhl1	80.54420	0.6183200	0.2156025	2.867870	0.0041325	0.9998239
ENSMUSG00000030878.12 Cdr2	10.26062	1.0041048	0.3504475	2.865208	0.0041674	0.9998239
ENSMUSG00000097392.5 Thoc2l	51.38669	0.4773075	0.1690408	2.823623	0.0047484	0.9998239
ENSMUSG00000083012.10 Fam220a	32.46395	0.8736736	0.3120714	2.799595	0.0051167	0.9998239
ENSMUSG00000027102.5 Hoxd8	120.35253	-0.4902413	0.1751403	-2.799135	0.0051240	0.9998239
ENSMUSG00000028580.16 Pum1	71.18988	-0.3325650	0.1202664	-2.765236	0.0056882	0.9998239
ENSMUSG00000062683.12 Atp5g2	28.39540	-0.8140740	0.2984961	-2.727252	0.0063864	0.9998239
ENSMUSG00000034800.16 Zfp661	13.08517	-0.7300091	0.2745900	-2.658542	0.0078480	0.9998239
ENSMUSG00000055491.15 Pprc1	33.78304	0.4933753	0.1872039	2.635496	0.0084014	0.9998239
ENSMUSG00000022837.15 Iqcb1	11.44752	0.8078495	0.3092115	2.612611	0.0089853	0.9998239
ENSMUSG00000035372.3 1810055G02Rik	12.01980	0.7490477	0.2921792	2.563659	0.0103575	0.9998239
ENSMUSG00000100826.7 Snhg14	30.51278	0.5609355	0.2191422	2.559687	0.0104766	0.9998239
ENSMUSG00000044876.16 Zfp444	18.12292	-0.6051971	0.2408291	-2.512974	0.0119718	0.9998239

dge1 <- dge
d1up <- rownames(subset(dge1,padj <= 0.05 & log2FoldChange > 0))
d1dn <- rownames(subset(dge1,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1,file="dge1.tsv",quote=FALSE,sep="\t")

Now paired.

dim(xx_atp_sed)

## [1] 9731   14

xx_atp_sed <- xx_atp_sed[which(rowMeans(xx_atp_sed)>=10),]
dim(xx_atp_sed)

## [1] 9731   14

ss_atp_sed$construct <- factor(ss_atp_sed$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_atp_sed , colData = ss_atp_sed,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000068614.8 Actc1	240.02972	-0.8076848	0.1585609	-5.093846	0.0000004	0.0034143
ENSMUSG00000025270.14 Alas2	32.23192	2.0342894	0.5660663	3.593730	0.0003260	0.9999466
ENSMUSG00000076613.5 Ighg2b	39.96923	-10.7411244	3.0090239	-3.569637	0.0003575	0.9999466
ENSMUSG00000020061.19 Mybpc1	1362.91124	0.6417473	0.1848772	3.471208	0.0005181	0.9999466
ENSMUSG00000073940.4 Hbb-bt	298.94874	1.8167623	0.5433036	3.343917	0.0008260	0.9999466
ENSMUSG00000113337.2 Gm19220	1733.35603	0.2757742	0.0839421	3.285292	0.0010188	0.9999466
ENSMUSG00000052305.7 Hbb-bs	1991.61520	1.8129798	0.5530722	3.278017	0.0010454	0.9999466
ENSMUSG00000069917.8 Hba-a2	937.22103	1.8750593	0.5833897	3.214077	0.0013086	0.9999466
ENSMUSG00000027102.5 Hoxd8	120.35253	-0.5118406	0.1593011	-3.213039	0.0013134	0.9999466
ENSMUSG00000096921.2 Gm26822	1354.43786	0.3876625	0.1251925	3.096531	0.0019580	0.9999466
ENSMUSG00000069919.8 Hba-a1	895.64500	1.7488079	0.5761464	3.035353	0.0024025	0.9999466
ENSMUSG00000029167.14 Ppargc1a	300.21147	0.4816588	0.1618032	2.976820	0.0029126	0.9999466
ENSMUSG00000035578.17 Iqcg	15.39691	1.0852294	0.3985490	2.722951	0.0064702	0.9999466
ENSMUSG00000031448.13 Adprhl1	80.54420	0.5352995	0.2090013	2.561226	0.0104303	0.9999466
ENSMUSG00000062683.12 Atp5g2	28.39540	-0.8246798	0.3241333	-2.544262	0.0109509	0.9999466
ENSMUSG00000000628.11 Hk2	692.60973	0.3237669	0.1281849	2.525781	0.0115442	0.9999466
ENSMUSG00000019997.12 Ccn2	46.67322	0.5839434	0.2331967	2.504081	0.0122770	0.9999466
ENSMUSG00002076161.1 7SK	11.52781	1.3024816	0.5216579	2.496812	0.0125315	0.9999466
ENSMUSG00000025754.12 Agbl1	17.65449	0.8661382	0.3479771	2.489066	0.0128079	0.9999466
ENSMUSG00000032754.15 Slc8b1	15.52746	-0.8842827	0.3555770	-2.486895	0.0128863	0.9999466

dge1p <- dge
d1pup <- rownames(subset(dge1p,padj <= 0.05 & log2FoldChange > 0))
d1pdn <- rownames(subset(dge1p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1p,file="dge1p.tsv",quote=FALSE,sep="\t")

Paired and excluding poor ATP-CL knockdown samples.

ss_atp_sed <- ss_atp[which(ss_atp$ex_sed == "Sed"),]
xx_atp_sed <- xx_atp[,which(colnames(xx_atp) %in% rownames(ss_atp_sed))]
ss_atp_sed$construct <- factor(ss_atp_sed$construct,levels=c("ctrl","kd"))
rpm_atp_sed <- apply(xx_atp_sed,2,function(x) { x / sum(x) }) * 1000000
ss_atp_sed$acly <- rpm_atp_sed[grep("Acly",rownames(rpm_atp_sed) ),]
animals <- unique(ss_atp_sed$animal_id)
animals_kd <- t( sapply(animals,function(a) {
  ss_atp_sed[ss_atp_sed$animal_id==a,"acly"]
} ) )
animals_kd2 <- animals_kd[,1] / animals_kd[,2]
# threshold
animals_kd2 <- animals_kd2[order(animals_kd2)]
barplot(animals_kd2, ylab="foldchange (kd/ctrl)",xlab="animal ID")

animals_kd2 <- animals_kd2[animals_kd2<0.7]
ss_atp_sed <- ss_atp_sed[ss_atp_sed$animal_id %in% names(animals_kd2),]
ss_atp_sed %>%
  kbl(caption = "Contrast 1: samples to use after filtering") %>%
  kable_paper("hover", full_width = F)

Contrast 1: samples to use after filtering

	Sample_ID	Orginal_ID	percent	num_reads	Gbases	target_gene	ex_sed	construct	animal_id	acly
5176140	5176140	1Left	0.3412	5860429	1.465107	ATP-CL	Sed	kd	1	6.577757
5176141	5176141	2Left	0.2909	4996480	1.249120	ATP-CL	Sed	kd	2	12.702420
5176142	5176142	3Left	0.3103	5329693	1.332423	ATP-CL	Sed	kd	3	10.516205
5176160	5176160	51Left	0.2578	4427956	1.106989	ATP-CL	Sed	kd	51	7.898470
5176168	5176168	1Right	0.3798	6523420	1.630855	ATP-CL	Sed	ctrl	1	11.831473
5176169	5176169	2Right	0.3058	5252401	1.313100	ATP-CL	Sed	ctrl	2	24.752642
5176170	5176170	3Right	0.3551	6099175	1.524794	ATP-CL	Sed	ctrl	3	16.030406
5176188	5176188	51Right	0.3238	5561568	1.390392	ATP-CL	Sed	ctrl	51	24.058521

xx_atp_sed <- xx_atp_sed[,colnames(xx_atp_sed) %in% rownames(ss_atp_sed)]
dim(xx_atp_sed)

## [1] 55357     8

xx_atp_sed <- xx_atp_sed[which(rowMeans(xx_atp_sed)>=10),]
dim(xx_atp_sed)

## [1] 9729    8

dds <- DESeqDataSetFromMatrix(countData = xx_atp_sed , colData = ss_atp_sed,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Sed mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000065987.13 Cd209b	25.676302	-4.1179786	0.7276676	-5.659148	0.0000000	0.0001480
ENSMUSG00000029167.14 Ppargc1a	420.766985	0.8668544	0.1601355	5.413255	0.0000001	0.0003011
ENSMUSG00000061780.7 Cfd	184.582980	-1.3178833	0.2500533	-5.270410	0.0000001	0.0004414
ENSMUSG00000035697.15 Arhgap45	33.253751	-2.3199014	0.5186477	-4.472981	0.0000077	0.0187615
ENSMUSG00000005540.11 Fcer2a	9.566438	-4.3039427	1.0741385	-4.006879	0.0000615	0.1005484
ENSMUSG00000045826.10 Ptprcap	14.438982	-3.2870059	0.8207193	-4.005031	0.0000620	0.1005484
ENSMUSG00000030278.12 Cidec	26.453596	-1.6790539	0.4550862	-3.689529	0.0002247	0.3116319
ENSMUSG00000043263.14 Ifi209	13.429215	-3.1238515	0.8583067	-3.639552	0.0002731	0.3116319
ENSMUSG00000059456.14 Ptk2b	9.576665	-3.5211905	0.9752598	-3.610515	0.0003056	0.3116319
ENSMUSG00000032312.8 Csk	23.507963	-1.7928583	0.4982524	-3.598293	0.0003203	0.3116319
ENSMUSG00000025270.14 Alas2	27.454533	2.0935708	0.6089061	3.438249	0.0005855	0.4913577
ENSMUSG00000022878.6 Adipoq	35.832112	-1.4871147	0.4337014	-3.428891	0.0006061	0.4913577
ENSMUSG00000102070.2 Gm28661	27125.714689	0.4011037	0.1199262	3.344586	0.0008241	0.6167098
ENSMUSG00000037337.12 Map4k1	9.774560	-2.9659084	0.8929763	-3.321374	0.0008958	0.6224846
ENSMUSG00000038642.11 Ctss	28.597913	-1.6992384	0.5162977	-3.291199	0.0009976	0.6470516
ENSMUSG00000000628.11 Hk2	775.096111	0.5495904	0.1716692	3.201451	0.0013674	0.8314494
ENSMUSG00000054435.17 Gimap4	16.878383	-1.9900137	0.6274111	-3.171786	0.0015150	0.8670523
ENSMUSG00000025479.10 Cyp2e1	31.369413	-1.3773149	0.4387245	-3.139362	0.0016932	0.8964599
ENSMUSG00000030577.15 Cd22	18.873361	-5.0512287	1.6140410	-3.129554	0.0017507	0.8964599
ENSMUSG00000020061.19 Mybpc1	1401.692364	0.5119868	0.1651623	3.099901	0.0019359	0.9406171

dge1p <- dge
d1pup <- rownames(subset(dge1p,padj <= 0.05 & log2FoldChange > 0))
d1pdn <- rownames(subset(dge1p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1p,file="dge1pkd.tsv",quote=FALSE,sep="\t")

DGE 2 Effect of ATP-CL KD in Ex animals

dim(xx_atp_ex)

## [1] 55357    14

xx_atp_ex <- xx_atp_ex[which(rowMeans(xx_atp_ex)>=10),]
dim(xx_atp_ex)

## [1] 9350   14

ss_atp_ex$construct <- factor(ss_atp_ex$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_atp_ex , colData = ss_atp_ex,
  design = ~ construct )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000025085.17 Ablim1	80.83082	-0.5830156	0.1363528	-4.275787	0.0000190	0.1780835
ENSMUSG00000032285.16 Dnaja4	176.93117	-0.5125999	0.1501048	-3.414946	0.0006379	0.9999758
ENSMUSG00000015882.19 Lcorl	13.27752	-0.9419652	0.2802011	-3.361747	0.0007745	0.9999758
ENSMUSG00000039126.11 Prune2	19.41197	-1.0216183	0.3201392	-3.191169	0.0014170	0.9999758
ENSMUSG00000035891.17 Cerk	15.22785	-0.8450373	0.2649997	-3.188824	0.0014285	0.9999758
ENSMUSG00000071748.3 Styx-ps	20.27826	-1.0192650	0.3208714	-3.176553	0.0014904	0.9999758
ENSMUSG00000042688.17 Mapk6	85.97710	-0.4433145	0.1404235	-3.156983	0.0015941	0.9999758
ENSMUSG00000057229.12 Dmac2	69.82397	-0.4132302	0.1348957	-3.063331	0.0021889	0.9999758
ENSMUSG00000082585.4 Gm15387	15.71188	-1.2743908	0.4173266	-3.053701	0.0022604	0.9999758
ENSMUSG00000025810.10 Nrp1	32.37937	-0.6343818	0.2085990	-3.041155	0.0023567	0.9999758
ENSMUSG00000032355.17 Mlip	116.21198	0.3613513	0.1197469	3.017625	0.0025476	0.9999758
ENSMUSG00000022228.15 Zscan26	49.86120	-0.4824735	0.1613334	-2.990537	0.0027849	0.9999758
ENSMUSG00000118841.1 Rn7s2	40.22956	0.7861496	0.2633795	2.984854	0.0028371	0.9999758
ENSMUSG00000118866.1 Rn7s1	40.22956	0.7861496	0.2633795	2.984854	0.0028371	0.9999758
ENSMUSG00000001630.14 Stk38l	29.16045	-0.5777016	0.1944549	-2.970877	0.0029695	0.9999758
ENSMUSG00000044788.11 Fads6	23.09655	-0.6487560	0.2213555	-2.930833	0.0033805	0.9999758
ENSMUSG00000022748.9 Cmss1	80.79703	0.4021062	0.1377294	2.919538	0.0035055	0.9999758
ENSMUSG00000047205.13 Dusp18	50.90863	-0.7855929	0.2701439	-2.908053	0.0036369	0.9999758
ENSMUSG00000085795.9 Zfp703	48.20385	-0.5600892	0.1937655	-2.890551	0.0038457	0.9999758
ENSMUSG00000044068.8 Zrsr1	76.79049	0.4587746	0.1594409	2.877395	0.0040097	0.9999758

dge2 <- dge
d2up <- rownames(subset(dge2,padj <= 0.05 & log2FoldChange > 0))
d2dn <- rownames(subset(dge2,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2,file="dge2.tsv",quote=FALSE,sep="\t")

Now paired.

dim(xx_atp_ex)

## [1] 9350   14

xx_atp_ex <- xx_atp_ex[which(rowMeans(xx_atp_ex)>=10),]
dim(xx_atp_ex)

## [1] 9350   14

ss_atp_ex$construct <- factor(ss_atp_ex$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_atp_ex , colData = ss_atp_ex,
  design =  ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000028341.10 Nr4a3	97.32005	-0.7833945	0.1926513	-4.066386	0.0000477	0.4464422
ENSMUSG00000025085.17 Ablim1	80.83082	-0.5942138	0.1582008	-3.756072	0.0001726	0.8069095
ENSMUSG00000082585.4 Gm15387	15.71188	-1.5517443	0.4455379	-3.482856	0.0004961	0.9999912
ENSMUSG00000010064.16 Slc38a3	178.04272	-0.5192180	0.1557997	-3.332601	0.0008604	0.9999912
ENSMUSG00000071748.3 Styx-ps	20.27826	-1.0167107	0.3060107	-3.322468	0.0008923	0.9999912
ENSMUSG00000113337.2 Gm19220	1666.77920	0.4379943	0.1351313	3.241250	0.0011901	0.9999912
ENSMUSG00000093752.2 Gm20716	34.45183	0.9401689	0.3006225	3.127407	0.0017636	0.9999912
ENSMUSG00000032285.16 Dnaja4	176.93117	-0.4892699	0.1643933	-2.976215	0.0029183	0.9999912
ENSMUSG00000069917.8 Hba-a2	245.76385	0.4864669	0.1666744	2.918667	0.0035153	0.9999912
ENSMUSG00000022383.14 Ppara	26.30567	-0.9094578	0.3246230	-2.801581	0.0050853	0.9999912
ENSMUSG00000000628.11 Hk2	780.37486	-0.3639162	0.1302654	-2.793652	0.0052117	0.9999912
ENSMUSG00000047205.13 Dusp18	50.90863	-0.6901430	0.2482482	-2.780052	0.0054350	0.9999912
ENSMUSG00000118841.1 Rn7s2	40.22956	0.7232683	0.2616802	2.763940	0.0057108	0.9999912
ENSMUSG00000118866.1 Rn7s1	40.22956	0.7232683	0.2616802	2.763940	0.0057108	0.9999912
ENSMUSG00000022237.18 Ankrd33b	92.95912	-0.5184671	0.1894808	-2.736252	0.0062143	0.9999912
ENSMUSG00000025608.10 Podxl	69.86122	-0.5189870	0.1907499	-2.720771	0.0065130	0.9999912
ENSMUSG00000042569.14 Dhrs7b	50.32567	0.5263753	0.1937499	2.716777	0.0065921	0.9999912
ENSMUSG00000039126.11 Prune2	19.41197	-0.9819264	0.3705711	-2.649765	0.0080548	0.9999912
ENSMUSG00000105096.2 Gbp10	35.38626	-1.3776272	0.5209614	-2.644394	0.0081837	0.9999912
ENSMUSG00000083365.2 Gm12895	2462.34239	0.7707310	0.2921792	2.637871	0.0083428	0.9999912

dge2p <- dge
d2pup <- rownames(subset(dge2p,padj <= 0.05 & log2FoldChange > 0))
d2pdn <- rownames(subset(dge2p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2p,file="dge2p.tsv",quote=FALSE,sep="\t")

Paired and excluding poor ATP-CL knockdown samples.

ss_atp_ex <- ss_atp[which(ss_atp$ex_sed == "Ex"),]
xx_atp_ex <- xx_atp[,which(colnames(xx_atp) %in% rownames(ss_atp_ex))]
ss_atp_ex$construct <- factor(ss_atp_ex$construct,levels=c("ctrl","kd"))
rpm_atp_ex <- apply(xx_atp_ex,2,function(x) { x / sum(x) }) * 1000000
ss_atp_ex$acly <- rpm_atp_ex[grep("Acly",rownames(rpm_atp_ex) ),]
animals <- unique(ss_atp_ex$animal_id)
animals_kd <- t( sapply(animals,function(a) {
  ss_atp_ex[ss_atp_ex$animal_id==a,"acly"]
} ) )
animals_kd2 <- animals_kd[,1] / animals_kd[,2]
# threshold
animals_kd2 <- animals_kd2[order(animals_kd2)]
barplot(animals_kd2, ylab="foldchange (kd/ctrl)",xlab="animal ID")

animals_kd2 <- animals_kd2[animals_kd2<0.7]
ss_atp_ex <- ss_atp_ex[ss_atp_ex$animal_id %in% names(animals_kd2),]
ss_atp_ex %>%
  kbl(caption = "Contrast 2: samples to use after filtering") %>%
  kable_paper("hover", full_width = F)

Contrast 2: samples to use after filtering

	Sample_ID	Orginal_ID	percent	num_reads	Gbases	target_gene	ex_sed	construct	animal_id	acly
5176156	5176156	43Left	0.3414	5863864	1.465966	ATP-CL	Ex	kd	43	8.606060
5176157	5176157	44Left	0.3028	5200873	1.300218	ATP-CL	Ex	kd	44	4.231777
5176184	5176184	43Right	0.3843	6600712	1.650178	ATP-CL	Ex	ctrl	43	14.587129
5176185	5176185	44Right	0.3362	5774550	1.443637	ATP-CL	Ex	ctrl	44	7.193799

xx_atp_ex <- xx_atp_ex[,colnames(xx_atp_ex) %in% rownames(ss_atp_ex)]
dim(xx_atp_ex)

## [1] 55357     4

xx_atp_ex <- xx_atp_ex[which(rowMeans(xx_atp_ex)>=10),]
dim(xx_atp_ex)

## [1] 9475    4

dds <- DESeqDataSetFromMatrix(countData = xx_atp_ex , colData = ss_atp_ex,
  design =  ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ATP-CL KD in Ex mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000083365.2 Gm12895	3320.72911	1.8487260	0.2943122	6.281512	0.0000000	0.0000016
ENSMUSG00000083558.2 Gm12896	3320.72911	1.8487260	0.2943122	6.281512	0.0000000	0.0000016
ENSMUSG00000059741.14 Myl3	169.57857	-2.4007884	0.4256130	-5.640778	0.0000000	0.0000535
ENSMUSG00000053093.17 Myh7	81.73305	-3.6629409	0.6560636	-5.583210	0.0000000	0.0000559
ENSMUSG00002075188.1 ENSMUSG00002075188	4601.90355	1.5278415	0.2805590	5.445704	0.0000001	0.0000698
ENSMUSG00002075524.1 ENSMUSG00002075524	4601.90355	1.5278415	0.2805590	5.445704	0.0000001	0.0000698
ENSMUSG00002076173.1 ENSMUSG00002076173	4601.90355	1.5278415	0.2805590	5.445704	0.0000001	0.0000698
ENSMUSG00000098973.3 Mir6236	13930.15378	1.3857210	0.2820092	4.913743	0.0000009	0.0010583
ENSMUSG00000073867.3 AA474408	9500.70983	1.3969417	0.2993601	4.666425	0.0000031	0.0032266
ENSMUSG00000096921.2 Gm26822	1955.56286	1.5386940	0.3389079	4.540153	0.0000056	0.0052006
ENSMUSG00000035202.9 Lars2	13783.13621	1.3837694	0.3058008	4.525068	0.0000060	0.0052006
ENSMUSG00000119584.1 Rn18s-rs5	74422.23486	1.1701708	0.3013001	3.883739	0.0001029	0.0812184
ENSMUSG00000013936.13 Myl2	62.47812	-2.1061612	0.6139691	-3.430403	0.0006027	0.4144974
ENSMUSG00000027077.8 Smtnl1	98.19440	-1.6195454	0.4727164	-3.426040	0.0006124	0.4144974
ENSMUSG00000030098.16 Grip2	398.61078	1.1988244	0.3557988	3.369389	0.0007534	0.4758663
ENSMUSG00000116207.2 Nnt	33.37880	-2.6833909	0.8452096	-3.174823	0.0014993	0.8878546
ENSMUSG00000118642.2 AY036118	93.16534	1.3572950	0.4524350	2.999978	0.0027000	0.9999040
ENSMUSG00000118841.1 Rn7s2	56.09762	1.5921595	0.5348912	2.976604	0.0029146	0.9999040
ENSMUSG00000118866.1 Rn7s1	56.09762	1.5921595	0.5348912	2.976604	0.0029146	0.9999040
ENSMUSG00000028464.17 Tpm2	5478.02016	-0.8790404	0.2989910	-2.940023	0.0032819	0.9999040

dge2p <- dge
d2pup <- rownames(subset(dge2p,padj <= 0.05 & log2FoldChange > 0))
d2pdn <- rownames(subset(dge2p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2p,file="dge2pkd.tsv",quote=FALSE,sep="\t")

DGE 3 Effect of ACSS2 KD in Sed animals

dim(xx_acs_sed)

## [1] 55357    14

xx_acs_sed <- xx_acs_sed[which(rowMeans(xx_acs_sed)>=10),]
dim(xx_acs_sed)

## [1] 9476   14

ss_acs_sed$construct <- factor(ss_acs_sed$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_acs_sed , colData = ss_acs_sed,
  design = ~ construct )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 6 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ACSS2 KD in Sed mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ACSS2 KD in Sed mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000027605.19 Acss2	69.42746	-0.8668631	0.1416468	-6.119891	0.0000000	0.0000089
ENSMUSG00000025813.15 Homer2	65.86986	-0.7482197	0.1965904	-3.805983	0.0001412	0.6692053
ENSMUSG00000027452.12 Acss1	57.77032	-0.7253830	0.1974947	-3.672924	0.0002398	0.7574204
ENSMUSG00000064356.1 mt-Atp8	526.99810	-0.7406859	0.2078583	-3.563418	0.0003661	0.8671906
ENSMUSG00000041702.8 Btbd7	59.05674	-0.4495734	0.1285848	-3.496319	0.0004717	0.8940114
ENSMUSG00000003808.19 Farsa	44.43966	0.5460277	0.1622555	3.365235	0.0007648	0.9999088
ENSMUSG00000030214.8 Plbd1	45.77055	-0.5496793	0.1641406	-3.348831	0.0008115	0.9999088
ENSMUSG00000036733.17 Rbm42	57.08366	0.4523195	0.1401298	3.227861	0.0012472	0.9999088
ENSMUSG00000031358.18 Msl3	84.25154	-0.3744717	0.1183747	-3.163445	0.0015591	0.9999088
ENSMUSG00000022881.16 Rfc4	20.28920	0.7072034	0.2260305	3.128797	0.0017552	0.9999088
ENSMUSG00000040740.8 Slc25a34	30.28045	-0.9033175	0.2900841	-3.113985	0.0018458	0.9999088
ENSMUSG00000039601.17 Rcan2	110.79367	-0.4653798	0.1501885	-3.098637	0.0019441	0.9999088
ENSMUSG00000102070.2 Gm28661	27794.22808	-0.4702067	0.1525373	-3.082569	0.0020522	0.9999088
ENSMUSG00000041351.17 Rap1gap	23.68455	-0.6884797	0.2242583	-3.070030	0.0021404	0.9999088
ENSMUSG00000096606.3 Tpbgl	187.25003	0.3936662	0.1283491	3.067153	0.0021611	0.9999088
ENSMUSG00000047044.8 D030056L22Rik	10.68044	0.9605605	0.3163285	3.036592	0.0023927	0.9999088
ENSMUSG00000022995.17 Enah	71.52311	-0.5134940	0.1700666	-3.019369	0.0025330	0.9999088
ENSMUSG00000044968.17 Napepld	36.01436	-0.5155553	0.1713665	-3.008495	0.0026255	0.9999088
ENSMUSG00000046058.8 Eid2	10.92884	0.9477383	0.3159202	2.999930	0.0027004	0.9999088
ENSMUSG00000047371.8 Zfp768	70.75257	0.3599124	0.1203944	2.989444	0.0027949	0.9999088

dge3 <- dge
d3up <- rownames(subset(dge3,padj <= 0.05 & log2FoldChange > 0))
d3dn <- rownames(subset(dge3,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3,file="dge3.tsv",quote=FALSE,sep="\t")

Now paired.

dim(xx_acs_sed)

## [1] 9476   14

xx_acs_sed <- xx_acs_sed[which(rowMeans(xx_acs_sed)>=10),]
dim(xx_acs_sed)

## [1] 9476   14

ss_acs_sed$construct <- factor(ss_acs_sed$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_acs_sed , colData = ss_acs_sed,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ACSS2 KD in Sed mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ACSS2 KD in Sed mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000102070.2 Gm28661	27794.22808	-0.4675174	0.0789410	-5.922364	0.0000000	0.0000301
ENSMUSG00000069917.8 Hba-a2	279.43406	0.6008402	0.1100441	5.459994	0.0000000	0.0002256
ENSMUSG00000027605.19 Acss2	69.42746	-0.8639573	0.1785717	-4.838153	0.0000013	0.0038116
ENSMUSG00000052305.7 Hbb-bs	623.85277	0.4683109	0.0976216	4.797208	0.0000016	0.0038116
ENSMUSG00000027861.14 Casq2	48.50935	-0.9776945	0.2470343	-3.957728	0.0000757	0.1434024
ENSMUSG00000025813.15 Homer2	65.86986	-0.8013585	0.2186662	-3.664757	0.0002476	0.3910016
ENSMUSG00000116207.2 Nnt	40.82966	-0.7821841	0.2210471	-3.538541	0.0004023	0.5446607
ENSMUSG00000073940.4 Hbb-bt	103.59021	0.5525853	0.1587451	3.480961	0.0004996	0.5917983
ENSMUSG00000037940.18 Inpp4b	36.06628	-0.8991303	0.2715240	-3.311421	0.0009282	0.8729905
ENSMUSG00000078636.5 Gm7336	38.29271	0.8681289	0.2639698	3.288743	0.0010064	0.8729905
ENSMUSG00000041559.8 Fmod	47.92517	-0.6781218	0.2063178	-3.286782	0.0010134	0.8729905
ENSMUSG00000036585.17 Fgf1	66.08354	-0.5815154	0.1785826	-3.256283	0.0011288	0.8913857
ENSMUSG00000027452.12 Acss1	57.77032	-0.7124008	0.2233024	-3.190296	0.0014213	0.9996583
ENSMUSG00000040740.8 Slc25a34	30.28045	-0.9140674	0.2901894	-3.149899	0.0016333	0.9996583
ENSMUSG00000018381.16 Abi3	21.89584	2.1323222	0.6967907	3.060205	0.0022119	0.9996583
ENSMUSG00000039601.17 Rcan2	110.79367	-0.4721204	0.1568985	-3.009082	0.0026204	0.9996583
ENSMUSG00000069919.8 Hba-a1	261.25182	0.3968958	0.1319294	3.008396	0.0026263	0.9996583
ENSMUSG00000078881.11 Gm14434	11.08812	-1.3385889	0.4477214	-2.989781	0.0027918	0.9996583
ENSMUSG00000022995.17 Enah	71.52311	-0.5378175	0.1812178	-2.967796	0.0029994	0.9996583
ENSMUSG00000096606.3 Tpbgl	187.25003	0.3758852	0.1287073	2.920466	0.0034951	0.9996583

dge3p <- dge
d3pup <- rownames(subset(dge3p,padj <= 0.05 & log2FoldChange > 0))
d3pdn <- rownames(subset(dge3p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3p,file="dge3p.tsv",quote=FALSE,sep="\t")

Paired and excluding poor ACSS2 knockdown samples.

ss_acs_sed <- ss_acs[which(ss_acs$ex_sed == "Sed"),]
xx_acs_sed <- xx_acs[,which(colnames(xx_acs) %in% rownames(ss_acs_sed))]
ss_acs_sed$construct <- factor(ss_acs_sed$construct,levels=c("ctrl","kd"))
rpm_acs_sed <- apply(xx_acs_sed,2,function(x) { x / sum(x) }) * 1000000
ss_acs_sed$acly <- rpm_acs_sed[grep("Acly",rownames(rpm_acs_sed) ),]
animals <- unique(ss_acs_sed$animal_id)
animals_kd <- t( sapply(animals,function(a) {
  ss_acs_sed[ss_acs_sed$animal_id==a,"acly"]
} ) )
animals_kd2 <- animals_kd[,1] / animals_kd[,2]
# threshold
animals_kd2 <- animals_kd2[order(animals_kd2)]
barplot(animals_kd2, ylab="foldchange (kd/ctrl)",xlab="animal ID")

animals_kd2 <- animals_kd2[animals_kd2<0.7]
ss_acs_sed <- ss_acs_sed[ss_acs_sed$animal_id %in% names(animals_kd2),]
ss_acs_sed %>%
  kbl(caption = "Contrast 3: samples to use after filtering") %>%
  kable_paper("hover", full_width = F)

Contrast 3: samples to use after filtering

	Sample_ID	Orginal_ID	percent	num_reads	Gbases	target_gene	ex_sed	construct	animal_id	acly
5176159	5176159	46Left	0.3113	5346869	1.336717	ACSS2	Sed	kd	46	15.133754
5176167	5176167	64Left	0.3145	5401832	1.350458	ACSS2	Sed	kd	64	8.098522
5176187	5176187	46Right	0.3214	5520346	1.380086	ACSS2	Sed	ctrl	46	24.008242
5176195	5176195	64Right	0.2913	5003350	1.250838	ACSS2	Sed	ctrl	64	15.587650

xx_acs_sed <- xx_acs_sed[,colnames(xx_acs_sed) %in% rownames(ss_acs_sed)]
dim(xx_acs_sed)

## [1] 55357     4

xx_acs_sed <- xx_acs_sed[which(rowMeans(xx_acs_sed)>=10),]
dim(xx_acs_sed)

## [1] 9546    4

dds <- DESeqDataSetFromMatrix(countData = xx_acs_sed , colData = ss_acs_sed,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 3: Effect of ACSS2 KD in Sed mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 3: Effect of ACSS2 KD in Sed mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000053093.17 Myh7	460.75798	-2.2320897	0.3789115	-5.890794	0.0000000	0.0000367
ENSMUSG00000102070.2 Gm28661	26176.52819	-0.5823278	0.1171265	-4.971783	0.0000007	0.0031664
ENSMUSG00000026418.17 Tnni1	119.17888	-1.9191404	0.4036509	-4.754455	0.0000020	0.0063316
ENSMUSG00000091898.9 Tnnc1	151.70607	-2.4704699	0.5617535	-4.397783	0.0000109	0.0260994
ENSMUSG00000030470.16 Csrp3	111.21503	-1.5055543	0.3589080	-4.194819	0.0000273	0.0521382
ENSMUSG00000023092.17 Fhl1	1602.56690	-0.5978935	0.1536208	-3.892008	0.0000994	0.1581741
ENSMUSG00000096921.2 Gm26822	875.29562	0.5963240	0.1693524	3.521201	0.0004296	0.5858464
ENSMUSG00000113337.2 Gm19220	1356.64603	0.5374967	0.1544134	3.480893	0.0004997	0.5963203
ENSMUSG00000057003.13 Myh4	71922.56902	0.3961550	0.1163337	3.405333	0.0006608	0.6872736
ENSMUSG00000025172.4 Ankrd2	190.76675	-1.0406400	0.3077120	-3.381863	0.0007200	0.6872736
ENSMUSG00000064354.1 mt-Co2	24042.99202	0.4239946	0.1314494	3.225534	0.0012574	0.9977401
ENSMUSG00000022665.16 Ccdc80	51.57202	-1.6993278	0.5297595	-3.207734	0.0013378	0.9977401
ENSMUSG00000101111.2 Gm28437	12153.70164	-0.9314483	0.2907802	-3.203273	0.0013587	0.9977401
ENSMUSG00000101249.2 Gm29216	9916.70302	-0.3849950	0.1268330	-3.035448	0.0024018	0.9999752
ENSMUSG00000035202.9 Lars2	7180.23515	0.3673083	0.1214990	3.023140	0.0025017	0.9999752
ENSMUSG00000027861.14 Casq2	53.44469	-1.4709668	0.5051865	-2.911730	0.0035943	0.9999752
ENSMUSG00000064179.14 Tnnt1	210.69089	-1.6532412	0.5782913	-2.858838	0.0042520	0.9999752
ENSMUSG00000063954.8 H2ac19	50.50615	1.5187090	0.5323966	2.852590	0.0043365	0.9999752
ENSMUSG00000030278.12 Cidec	18.72607	-2.6910466	0.9559968	-2.814912	0.0048791	0.9999752
ENSMUSG00000079017.4 Ifi27l2a	36.79955	-1.7482650	0.6254293	-2.795304	0.0051851	0.9999752

dge3p <- dge
d3pup <- rownames(subset(dge3p,padj <= 0.05 & log2FoldChange > 0))
d3pdn <- rownames(subset(dge3p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3p,file="dge3pkd.tsv",quote=FALSE,sep="\t")

DGE 4 Effect of ACSS2 KD in Ex animals

dim(xx_acs_ex)

## [1] 55357    14

xx_acs_ex <- xx_acs_ex[which(rowMeans(xx_acs_ex)>=10),]
dim(xx_acs_ex)

## [1] 9586   14

ss_acs_ex$construct <- factor(ss_acs_ex$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_acs_ex , colData = ss_acs_ex,
  design = ~ construct )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 77 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ACSS2 KD in Ex mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ACSS2 KD in Ex mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000082368.2 Gm11225	53.815179	0.8217104	0.1901818	4.320657	0.0000156	0.0839587
ENSMUSG00000093445.8 Lrch4	25.392993	-1.0098550	0.2351565	-4.294396	0.0000175	0.0839587
ENSMUSG00000064356.1 mt-Atp8	452.241900	-0.7467557	0.1870716	-3.991818	0.0000656	0.2095137
ENSMUSG00000032280.17 Tle3	11.907164	1.4515449	0.4195112	3.460086	0.0005400	0.9998986
ENSMUSG00000036873.14 2410004B18Rik	50.647219	-0.4862178	0.1407441	-3.454623	0.0005511	0.9998986
ENSMUSG00000102275.2 Gm37144	17.966059	1.0407117	0.3068971	3.391077	0.0006962	0.9998986
ENSMUSG00000053877.13 Srcap	103.651371	0.4004359	0.1266557	3.161611	0.0015690	0.9998986
ENSMUSG00000018572.11 Phf23	28.431285	0.6132548	0.1943046	3.156152	0.0015987	0.9998986
ENSMUSG00000056569.12 Mpz	14.725827	-0.9203541	0.2944868	-3.125281	0.0017764	0.9998986
ENSMUSG00000044037.16 Als2cl	23.569078	0.8052763	0.2627649	3.064627	0.0021794	0.9998986
ENSMUSG00000053398.12 Phgdh	8.291418	2.2875045	0.7475731	3.059908	0.0022141	0.9998986
ENSMUSG00000026675.13 Hsd17b7	42.227514	-0.5582106	0.1841540	-3.031216	0.0024357	0.9998986
ENSMUSG00000043801.8 Oaz1-ps	37.221298	-0.6260365	0.2093052	-2.991023	0.0027804	0.9998986
ENSMUSG00000039903.19 Eva1c	10.727972	1.2005710	0.4051821	2.963040	0.0030462	0.9998986
ENSMUSG00000105008.2 Gm43652	28.042631	0.8329197	0.2836187	2.936758	0.0033166	0.9998986
ENSMUSG00000025153.10 Fasn	36.661576	1.1678799	0.4021245	2.904274	0.0036811	0.9998986
ENSMUSG00000035495.16 Tstd2	18.669500	-0.8030495	0.2767248	-2.901979	0.0037081	0.9998986
ENSMUSG00000113337.2 Gm19220	1689.618927	0.5417903	0.1877987	2.884953	0.0039147	0.9998986
ENSMUSG00000071847.14 Apcdd1	9.474889	1.1876643	0.4159372	2.855393	0.0042984	0.9998986
ENSMUSG00000047153.5 Khnyn	32.188335	0.5299053	0.1902388	2.785475	0.0053449	0.9998986

dge4 <- dge
d4up <- rownames(subset(dge4,padj <= 0.05 & log2FoldChange > 0))
d4dn <- rownames(subset(dge4,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4,file="dge4.tsv",quote=FALSE,sep="\t")

Now paired.

dim(xx_acs_ex)

## [1] 9586   14

xx_acs_ex <- xx_acs_ex[which(rowMeans(xx_acs_ex)>=10),]
dim(xx_acs_ex)

## [1] 9586   14

ss_acs_ex$construct <- factor(ss_acs_ex$construct,levels=c("ctrl","kd"))

dds <- DESeqDataSetFromMatrix(countData = xx_acs_ex , colData = ss_acs_ex,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 1: Effect of ACSS2 KD in Ex mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of ACSS2 KD in Ex mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000029720.10 Gm20605	30.01174	1.2205191	0.3091525	3.947952	0.0000788	0.3796991
ENSMUSG00000051985.13 Igfn1	35.55728	1.2241701	0.3180466	3.849028	0.0001186	0.3796991
ENSMUSG00000082368.2 Gm11225	53.81518	0.8378147	0.2200450	3.807470	0.0001404	0.3796991
ENSMUSG00000093445.8 Lrch4	25.39299	-1.0575579	0.2799657	-3.777455	0.0001584	0.3796991
ENSMUSG00000064356.1 mt-Atp8	452.24190	-0.7235168	0.2172954	-3.329646	0.0008696	0.9998823
ENSMUSG00000067719.6 Gm10221	13.45493	-4.1371093	1.2560411	-3.293769	0.0009885	0.9998823
ENSMUSG00000110631.2 Gm42047	11.78495	1.4252824	0.4328984	3.292418	0.0009933	0.9998823
ENSMUSG00000105008.2 Gm43652	28.04263	0.8793884	0.2696604	3.261096	0.0011098	0.9998823
ENSMUSG00000039903.19 Eva1c	10.72797	1.3296421	0.4227682	3.145085	0.0016604	0.9998823
ENSMUSG00000102275.2 Gm37144	17.96606	1.0165016	0.3255659	3.122261	0.0017947	0.9998823
ENSMUSG00000052305.7 Hbb-bs	522.27262	0.3335285	0.1079414	3.089903	0.0020022	0.9998823
ENSMUSG00000026395.17 Ptprc	38.48284	1.0680438	0.3668535	2.911364	0.0035985	0.9998823
ENSMUSG00000023064.5 Sncg	31.92579	-0.9801378	0.3384897	-2.895621	0.0037841	0.9998823
ENSMUSG00000044037.16 Als2cl	23.56908	0.8157478	0.2830158	2.882340	0.0039473	0.9998823
ENSMUSG00000053877.13 Srcap	103.65137	0.4169147	0.1518529	2.745517	0.0060416	0.9998823
ENSMUSG00000043801.8 Oaz1-ps	37.22130	-0.6378066	0.2337079	-2.729076	0.0063512	0.9998823
ENSMUSG00000035495.16 Tstd2	18.66950	-0.8501129	0.3220916	-2.639352	0.0083065	0.9998823
ENSMUSG00000056569.12 Mpz	14.72583	-0.9361332	0.3583965	-2.612005	0.0090013	0.9998823
ENSMUSG00000026675.13 Hsd17b7	42.22751	-0.5700170	0.2202220	-2.588374	0.0096430	0.9998823
ENSMUSG00000018572.11 Phf23	28.43129	0.6331481	0.2494780	2.537891	0.0111523	0.9998823

dge4p <- dge
d4pup <- rownames(subset(dge4p,padj <= 0.05 & log2FoldChange > 0))
d4pdn <- rownames(subset(dge4p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4p,file="dge4p.tsv",quote=FALSE,sep="\t")

Paired and excluding poor ACSS2 knockdown samples.

ss_acs_ex <- ss_acs[which(ss_acs$ex_sed == "Ex"),]
xx_acs_ex <- xx_acs[,which(colnames(xx_acs) %in% rownames(ss_acs_ex))]
ss_acs_ex$construct <- factor(ss_acs_ex$construct,levels=c("ctrl","kd"))
rpm_acs_ex <- apply(xx_acs_ex,2,function(x) { x / sum(x) }) * 1000000
ss_acs_ex$acss2 <- rpm_acs_ex[grep("Acss2$",rownames(rpm_acs_ex) ),]
animals <- unique(ss_acs_ex$animal_id)
animals_kd <- t( sapply(animals,function(a) {
  ss_acs_ex[ss_acs_ex$animal_id==a,"acss2"]
} ) )
animals_kd2 <- animals_kd[,1] / animals_kd[,2]
# threshold
animals_kd2 <- animals_kd2[order(animals_kd2)]
barplot(animals_kd2, ylab="foldchange (kd/ctrl)",xlab="animal ID")

animals_kd2 <- animals_kd2[animals_kd2<0.7]
ss_acs_ex <- ss_acs_ex[ss_acs_ex$animal_id %in% names(animals_kd2),]
ss_acs_ex %>%
  kbl(caption = "Contrast 4: samples to use after filtering") %>%
  kable_paper("hover", full_width = F)

Contrast 4: samples to use after filtering

	Sample_ID	Orginal_ID	percent	num_reads	Gbases	target_gene	ex_sed	construct	animal_id	acss2
5176158	5176158	45Left	0.2799	4807544	1.201886	ACSS2	Ex	kd	45	34.86514
5176166	5176166	63Left	0.3385	5814054	1.453514	ACSS2	Ex	kd	63	11.80049
5176186	5176186	45Right	0.3223	5535804	1.383951	ACSS2	Ex	ctrl	45	50.04461
5176194	5176194	63Right	0.3239	5563286	1.390821	ACSS2	Ex	ctrl	63	41.01297

xx_acs_ex <- xx_acs_ex[,colnames(xx_acs_ex) %in% rownames(ss_acs_ex)]
dim(xx_acs_ex)

## [1] 55357     4

xx_acs_ex <- xx_acs_ex[which(rowMeans(xx_acs_ex)>=10),]
dim(xx_acs_ex)

## [1] 9476    4

dds <- DESeqDataSetFromMatrix(countData = xx_acs_ex , colData = ss_acs_ex,
  design = ~ animal_id + construct )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 3: Effect of ACSS2 KD in Sed mice (paired)") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 3: Effect of ACSS2 KD in Sed mice (paired)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000113337.2 Gm19220	1800.41139	1.1283003	0.2102351	5.366850	0.0000001	0.0004354
ENSMUSG00000059741.14 Myl3	364.31913	-2.9387768	0.5852930	-5.021035	0.0000005	0.0013964
ENSMUSG00000030246.13 Ldhb	273.16726	-2.5724186	0.5278060	-4.873795	0.0000011	0.0019066
ENSMUSG00000064356.1 mt-Atp8	506.92202	-1.2961419	0.2692002	-4.814788	0.0000015	0.0019066
ENSMUSG00000102070.2 Gm28661	28770.87117	-0.9315054	0.1948820	-4.779844	0.0000018	0.0019066
ENSMUSG00000101249.2 Gm29216	10067.31050	-0.8895674	0.1994323	-4.460497	0.0000082	0.0074056
ENSMUSG00000046598.16 Bdh1	77.19495	-2.7552518	0.6354453	-4.335939	0.0000145	0.0112670
ENSMUSG00000030541.17 Idh2	295.53712	-1.5415545	0.3627629	-4.249482	0.0000214	0.0132390
ENSMUSG00000064368.1 mt-Nd6	4395.76640	-0.8758682	0.2063630	-4.244307	0.0000219	0.0132390
ENSMUSG00000064363.1 mt-Nd4	38348.06968	-0.7746223	0.1973096	-3.925923	0.0000864	0.0446997
ENSMUSG00000096921.2 Gm26822	1039.58310	0.8618494	0.2208298	3.902777	0.0000951	0.0446997
ENSMUSG00000051985.13 Igfn1	74.06022	1.8674275	0.4806591	3.885140	0.0001023	0.0446997
ENSMUSG00000064337.1 mt-Rnr1	4882.53193	-0.8424655	0.2183169	-3.858912	0.0001139	0.0446997
ENSMUSG00000064345.1 mt-Nd2	27455.05218	-0.7925189	0.2055181	-3.856200	0.0001152	0.0446997
ENSMUSG00000029467.16 Atp2a2	301.98485	-1.0611726	0.2832481	-3.746441	0.0001794	0.0584462
ENSMUSG00000064357.1 mt-Atp6	47799.61865	-0.7496823	0.2001817	-3.745008	0.0001804	0.0584462
ENSMUSG00000091898.9 Tnnc1	61.83641	-1.9867293	0.5309827	-3.741609	0.0001828	0.0584462
ENSMUSG00000022186.15 Oxct1	747.90361	-0.9082557	0.2455065	-3.699517	0.0002160	0.0652110
ENSMUSG00000064367.1 mt-Nd5	22750.95033	-0.7283573	0.2010467	-3.622826	0.0002914	0.0833411
ENSMUSG00000064370.1 mt-Cytb	43173.96045	-0.7115485	0.1998128	-3.561077	0.0003693	0.0969039

dge4p <- dge
d4pup <- rownames(subset(dge4p,padj <= 0.05 & log2FoldChange > 0))
d4pdn <- rownames(subset(dge4p,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4p,file="dge4pkd.tsv",quote=FALSE,sep="\t")

DGE 5 Effect of Sed vs Ex (ATP-CL gene ctrl)

ss5 <- subset(ss,target_gene=="ATP-CL" & construct=="ctrl")
xx5 <- xx[,colnames(xx) %in% ss5$Sample_ID ]
dim(xx5)

## [1] 55357    14

xx5 <- xx5[which(rowMeans(xx5)>=10),]
dim(xx5)

## [1] 9706   14

ss5$ex_sed <- factor(ss5$ex_sed,levels=c("Sed","Ex"))

dds <- DESeqDataSetFromMatrix(countData = xx5 , colData = ss5, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 6 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ATP-CL ctrl mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ATP-CL ctrl mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000026390.8 Marco	11.04453	-22.1606764	2.9887768	-7.414631	0.0000000	0.0000000
ENSMUSG00000032942.15 Ucp3	588.61957	0.8472575	0.1662958	5.094882	0.0000003	0.0016935
ENSMUSG00000028525.17 Pde4b	152.02135	-1.1157382	0.2470264	-4.516676	0.0000063	0.0203237
ENSMUSG00000036528.17 Ppfibp2	37.79259	1.0786502	0.2454869	4.393922	0.0000111	0.0270127
ENSMUSG00000034714.10 Ttyh2	16.18307	-1.3794754	0.3437901	-4.012552	0.0000601	0.1062209
ENSMUSG00000035967.16 Ints6l	36.41940	-0.7739467	0.1938998	-3.991477	0.0000657	0.1062209
ENSMUSG00000013846.11 St3gal1	44.35171	-0.7814295	0.1981363	-3.943899	0.0000802	0.1111579
ENSMUSG00000062115.16 Rai1	61.79445	0.6282009	0.1626194	3.863012	0.0001120	0.1254995
ENSMUSG00000022139.18 Mbnl2	543.05690	0.4768490	0.1237396	3.853648	0.0001164	0.1254995
ENSMUSG00000056938.17 Acbd4	36.19872	-0.7004054	0.1867819	-3.749857	0.0001769	0.1717335
ENSMUSG00000023805.18 Synj2	71.89707	0.7133584	0.1932938	3.690539	0.0002238	0.1974550
ENSMUSG00000065987.13 Cd209b	15.70792	-4.8458234	1.3267070	-3.652520	0.0002597	0.2100376
ENSMUSG00000000628.11 Hk2	791.74664	0.6701418	0.1857764	3.607249	0.0003095	0.2310482
ENSMUSG00000020427.12 Igfbp3	22.43974	-0.8248116	0.2337759	-3.528215	0.0004184	0.2729760
ENSMUSG00000022096.15 Hr	221.87809	0.4184070	0.1189893	3.516340	0.0004375	0.2729760
ENSMUSG00000034563.17 Ccpg1	142.90633	0.4455308	0.1269722	3.508884	0.0004500	0.2729760
ENSMUSG00000029864.8 Gstk1	62.30386	0.6672959	0.1932860	3.452375	0.0005557	0.2773670
ENSMUSG00000076613.5 Ighg2b	33.61829	-8.9378416	2.5902501	-3.450571	0.0005594	0.2773670
ENSMUSG00000040945.14 Rcc2	36.21626	-0.9109950	0.2660529	-3.424113	0.0006168	0.2773670
ENSMUSG00000022994.10 Adcy6	66.01249	0.5611707	0.1640434	3.420866	0.0006242	0.2773670

dge5 <- dge
d5up <- rownames(subset(dge5,padj <= 0.05 & log2FoldChange > 0))
d5dn <- rownames(subset(dge5,padj <= 0.05 & log2FoldChange < 0))
write.table(dge5,file="dge5.tsv",quote=FALSE,sep="\t")

DGE 6 Effect of Sed vs Ex (ATP-CL gene KO)

ss6 <- subset(ss,target_gene=="ATP-CL" & construct=="kd")
xx6 <- xx[,colnames(xx) %in% ss6$Sample_ID ]
dim(xx6)

## [1] 55357    14

xx6 <- xx6[which(rowMeans(xx6)>=10),]
dim(xx6)

## [1] 9334   14

ss6$ex_sed <- factor(ss5$ex_sed,levels=c("Sed","Ex"))

dds <- DESeqDataSetFromMatrix(countData = xx6 , colData = ss6, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 15 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ctrl mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ctrl mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000032942.15 Ucp3	530.28447	0.8559652	0.1661521	5.151697	0.0000003	0.0018680
ENSMUSG00000033502.16 Cdc14a	61.67385	0.9422541	0.1870795	5.036652	0.0000005	0.0018680
ENSMUSG00000019841.16 Rev3l	117.77683	0.7344267	0.1499399	4.898142	0.0000010	0.0025432
ENSMUSG00000035900.19 Gramd4	36.77071	-0.8892073	0.1940430	-4.582528	0.0000046	0.0090568
ENSMUSG00000046532.9 Ar	130.39356	0.6896748	0.1532965	4.498960	0.0000068	0.0107702
ENSMUSG00000021453.3 Gadd45g	97.22503	-1.2495205	0.2850109	-4.384114	0.0000116	0.0153065
ENSMUSG00000068606.7 Gm4841	81.18066	0.7790695	0.1817200	4.287198	0.0000181	0.0203843
ENSMUSG00000045092.9 S1pr1	43.72419	-0.8297693	0.2086363	-3.977109	0.0000698	0.0687641
ENSMUSG00000089703.2 Gm15833	18.23155	1.0499076	0.2699663	3.889032	0.0001006	0.0881871
ENSMUSG00000042961.14 Egflam	62.26268	0.5465912	0.1425796	3.833585	0.0001263	0.0981998
ENSMUSG00000056708.6 Ier5	30.09298	-1.1982299	0.3146992	-3.807541	0.0001404	0.0981998
ENSMUSG00000026687.15 Aldh9a1	41.82470	1.2296232	0.3242663	3.792016	0.0001494	0.0981998
ENSMUSG00000071076.9 Jund	82.89084	-0.5941419	0.1579411	-3.761794	0.0001687	0.0994709
ENSMUSG00000028788.15 Ptp4a2	1935.54115	0.2308231	0.0617770	3.736390	0.0001867	0.0994709
ENSMUSG00000073940.4 Hbb-bt	279.24607	-2.5233884	0.6759659	-3.733011	0.0001892	0.0994709
ENSMUSG00000069919.8 Hba-a1	835.55888	-2.6180922	0.7148634	-3.662367	0.0002499	0.1231673
ENSMUSG00000052305.7 Hbb-bs	1885.48347	-2.4052127	0.6680646	-3.600270	0.0003179	0.1306935
ENSMUSG00000043795.11 Prr33	322.19559	0.2985847	0.0831465	3.591066	0.0003293	0.1306935
ENSMUSG00000069917.8 Hba-a2	888.90343	-2.5021560	0.6969068	-3.590374	0.0003302	0.1306935
ENSMUSG00000020053.19 Igf1	77.83589	0.9394266	0.2617234	3.589387	0.0003315	0.1306935

dge6 <- dge
d6up <- rownames(subset(dge6,padj <= 0.05 & log2FoldChange > 0))
d6dn <- rownames(subset(dge6,padj <= 0.05 & log2FoldChange < 0))
write.table(dge6,file="dge6.tsv",quote=FALSE,sep="\t")

Now to remove samples where the KD didn’t work.

ss6 <- subset(ss,target_gene=="ATP-CL" & construct=="kd")
xx6 <- xx[,colnames(xx) %in% ss6$Sample_ID ]

rpm6 <- apply(xx6,2, function(x) { x / sum(x) } ) *1000000
acly <- rpm6[grep("Acly",rownames(rpm6)),]
barplot(acly ,horiz=TRUE, las=1)

ss6$acly <- acly

ss6 <- subset(ss6,acly<11)
xx6 <- xx[,colnames(xx) %in% ss6$Sample_ID ]

dim(xx6)

## [1] 55357     7

xx6 <- xx6[which(rowMeans(xx6)>=10),]
dim(xx6)

## [1] 9342    7

ss6$ex_sed <- factor(ss6$ex_sed,levels=c("Sed","Ex"))

dds <- DESeqDataSetFromMatrix(countData = xx6 , colData = ss6, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ctrl mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ctrl mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000021453.3 Gadd45g	94.82049	-1.5595842	0.2520694	-6.187122	0.0000000	0.0000057
ENSMUSG00000026687.15 Aldh9a1	59.84100	1.5245276	0.3202402	4.760576	0.0000019	0.0090102
ENSMUSG00000033502.16 Cdc14a	71.68470	1.1331555	0.2588444	4.377747	0.0000120	0.0373127
ENSMUSG00000019841.16 Rev3l	129.36404	0.9581417	0.2227514	4.301395	0.0000170	0.0396099
ENSMUSG00000020053.19 Igf1	100.11299	1.0381293	0.2516485	4.125314	0.0000370	0.0691217
ENSMUSG00000019947.11 Arid5b	22.13908	-1.7515093	0.4324774	-4.049944	0.0000512	0.0797050
ENSMUSG00000039910.11 Cited2	106.83630	-0.9172648	0.2311770	-3.967802	0.0000725	0.0967351
ENSMUSG00000060548.14 Tnfrsf19	12.14674	-2.1964739	0.5885622	-3.731932	0.0001900	0.2217261
ENSMUSG00000024610.16 Cd74	110.12718	1.0935391	0.2980555	3.668911	0.0002436	0.2514013
ENSMUSG00000060586.12 H2-Eb1	47.02679	1.3483527	0.3701052	3.643161	0.0002693	0.2514013
ENSMUSG00000020692.15 Nle1	32.71497	-1.2785527	0.3726159	-3.431289	0.0006007	0.4794396
ENSMUSG00000064354.1 mt-Co2	23552.70744	0.6151048	0.1798154	3.420757	0.0006245	0.4794396
ENSMUSG00000091712.3 Sec14l5	69.41885	-0.9648875	0.2854493	-3.380241	0.0007242	0.4794396
ENSMUSG00000028982.10 Slc25a33	38.07485	-1.0006130	0.2961396	-3.378855	0.0007279	0.4794396
ENSMUSG00000032942.15 Ucp3	597.43687	0.8008578	0.2381223	3.363220	0.0007704	0.4794396
ENSMUSG00000066595.10 Flvcr1	44.75355	-1.1392611	0.3416503	-3.334582	0.0008543	0.4984174
ENSMUSG00000046532.9 Ar	145.17687	0.6439110	0.1997857	3.223009	0.0012685	0.6611672
ENSMUSG00000069014.5 Gm5641	83.13644	0.8060534	0.2502048	3.221575	0.0012749	0.6611672
ENSMUSG00000038594.10 Cep85l	138.59364	0.7229856	0.2256873	3.203484	0.0013578	0.6624075
ENSMUSG00000026380.11 Tfcp2l1	54.68889	0.9532620	0.2998634	3.178987	0.0014779	0.6624075

dge6 <- dge
d6up <- rownames(subset(dge6,padj <= 0.05 & log2FoldChange > 0))
d6dn <- rownames(subset(dge6,padj <= 0.05 & log2FoldChange < 0))
write.table(dge6,file="dge6kd.tsv",quote=FALSE,sep="\t")

DGE 7 Effect of Sed vs Ex (ACSS2 gene ctrl)

ss7 <- subset(ss,target_gene=="ACSS2" & construct=="ctrl")
xx7 <- xx[,colnames(xx) %in% ss7$Sample_ID ]
dim(xx7)

## [1] 55357    14

xx7 <- xx7[which(rowMeans(xx7)>=10),]
dim(xx7)

## [1] 9611   14

ss7$ex_sed <- factor(ss7$ex_sed,levels=c("Sed","Ex"))

dds <- DESeqDataSetFromMatrix(countData = xx7 , colData = ss7, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 68 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ACSS2 ctrl mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ACSS2 ctrl mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000048490.15 Nrip1	108.40303	1.1453462	0.1682215	6.808559	0.00e+00	0.0000001
ENSMUSG00000025648.19 Pfkfb4	192.64096	0.8052716	0.1213475	6.636079	0.00e+00	0.0000002
ENSMUSG00000026313.17 Hdac4	116.81696	0.8291389	0.1527317	5.428727	1.00e-07	0.0001818
ENSMUSG00000054978.4 Kbtbd13	46.09569	-0.9427582	0.1764110	-5.344099	1.00e-07	0.0001975
ENSMUSG00000028630.10 Dyrk2	148.52172	-0.6329280	0.1189320	-5.321764	1.00e-07	0.0001975
ENSMUSG00000032060.11 Cryab	563.48777	-0.8123409	0.1576692	-5.152186	3.00e-07	0.0003711
ENSMUSG00000038594.10 Cep85l	175.73954	1.0153997	0.1974658	5.142154	3.00e-07	0.0003711
ENSMUSG00000044499.12 Hs3st5	15.18664	-1.8695151	0.3652882	-5.117918	3.00e-07	0.0003711
ENSMUSG00000031700.12 Gpt2	882.14658	0.7971687	0.1639266	4.862963	1.20e-06	0.0011808
ENSMUSG00000026675.13 Hsd17b7	37.75391	1.1809043	0.2434369	4.850967	1.20e-06	0.0011808
ENSMUSG00000032449.14 Slc25a36	124.01497	-0.5620710	0.1198576	-4.689489	2.70e-06	0.0023930
ENSMUSG00000116207.2 Nnt	38.24236	-1.1923230	0.2554549	-4.667450	3.00e-06	0.0024425
ENSMUSG00000038418.8 Egr1	10.54214	-2.2738870	0.4919891	-4.621824	3.80e-06	0.0028122
ENSMUSG00000044167.7 Foxo1	69.62727	1.0930449	0.2454335	4.453527	8.40e-06	0.0057989
ENSMUSG00000074178.4 Gm10638	14.86624	1.5717986	0.3546416	4.432076	9.30e-06	0.0059800
ENSMUSG00000041235.13 Chd7	101.68882	0.6300526	0.1426353	4.417227	1.00e-05	0.0060054
ENSMUSG00000040584.9 Abcb1a	19.91313	-1.0879983	0.2504946	-4.343400	1.40e-05	0.0079316
ENSMUSG00000022180.8 Slc7a8	32.59774	0.9410382	0.2224329	4.230661	2.33e-05	0.0124412
ENSMUSG00000048406.6 B330016D10Rik	24.51338	0.8588648	0.2051325	4.186878	2.83e-05	0.0143061
ENSMUSG00000024789.14 Jak2	239.78103	0.7129541	0.1714991	4.157188	3.22e-05	0.0154828

dge7 <- dge
d7up <- rownames(subset(dge7,padj <= 0.05 & log2FoldChange > 0))
d7dn <- rownames(subset(dge7,padj <= 0.05 & log2FoldChange < 0))
write.table(dge7,file="dge7.tsv",quote=FALSE,sep="\t")

DGE 8 Effect of Sed vs Ex (ACSS2 gene kd)

ss8 <- subset(ss,target_gene=="ACSS2" & construct=="kd")
xx8 <- xx[,colnames(xx) %in% ss8$Sample_ID ]
dim(xx8)

## [1] 55357    14

xx8 <- xx8[which(rowMeans(xx8)>=10),]
dim(xx8)

## [1] 9441   14

ss8$ex_sed <- factor(ss8$ex_sed,levels=c("Sed","Ex"))
dds <- DESeqDataSetFromMatrix(countData = xx8 , colData = ss8, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 30 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ACSS2 kd mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ACSS2 kd mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000038594.10 Cep85l	166.06460	1.0213643	0.1960430	5.209901	0.0000002	0.0017838
ENSMUSG00000028630.10 Dyrk2	137.87586	-0.5337978	0.1184105	-4.508028	0.0000065	0.0261023
ENSMUSG00000026048.17 Ercc5	37.50491	0.8656112	0.1943734	4.453342	0.0000085	0.0261023
ENSMUSG00000026313.17 Hdac4	104.67492	0.7533275	0.1713917	4.395356	0.0000111	0.0261023
ENSMUSG00000053877.13 Srcap	95.59953	0.5627091	0.1335740	4.212714	0.0000252	0.0476433
ENSMUSG00000029456.13 Acad10	11.67196	1.1475550	0.2983556	3.846266	0.0001199	0.1692938
ENSMUSG00000031482.10 Slc25a15	32.23529	-0.7281563	0.1920233	-3.792020	0.0001494	0.1692938
ENSMUSG00000027620.17 Rbm39	269.74886	0.3154454	0.0833915	3.782703	0.0001551	0.1692938
ENSMUSG00000005413.9 Hmox1	29.87911	1.5436575	0.4105281	3.760175	0.0001698	0.1692938
ENSMUSG00000007817.16 Zmiz1	92.24660	0.6469107	0.1733100	3.732680	0.0001895	0.1692938
ENSMUSG00000068742.12 Cry2	84.88694	0.6067031	0.1630927	3.719990	0.0001992	0.1692938
ENSMUSG00000044791.17 Setd2	56.40038	0.5620546	0.1526585	3.681776	0.0002316	0.1692938
ENSMUSG00000030123.16 Plxnd1	76.28548	-0.5020969	0.1364345	-3.680132	0.0002331	0.1692938
ENSMUSG00000027253.16 Lrp4	19.61840	-0.9078061	0.2532588	-3.584500	0.0003377	0.2277468
ENSMUSG00000031626.19 Sorbs2	220.27380	0.5683906	0.1635272	3.475818	0.0005093	0.2737386
ENSMUSG00000048490.15 Nrip1	99.88082	0.7183240	0.2075901	3.460300	0.0005396	0.2737386
ENSMUSG00000024011.18 Pi16	60.06747	0.7321150	0.2123999	3.446870	0.0005671	0.2737386
ENSMUSG00000074178.4 Gm10638	13.58386	1.3854317	0.4022946	3.443823	0.0005736	0.2737386
ENSMUSG00000091712.3 Sec14l5	74.85437	-0.7876107	0.2288124	-3.442168	0.0005771	0.2737386
ENSMUSG00000033088.20 Triobp	77.75599	-0.4604689	0.1338242	-3.440848	0.0005799	0.2737386

dge8 <- dge
d8up <- rownames(subset(dge8,padj <= 0.05 & log2FoldChange > 0))
d8dn <- rownames(subset(dge8,padj <= 0.05 & log2FoldChange < 0))
write.table(dge8,file="dge8.tsv",quote=FALSE,sep="\t")

Now to remove samples with RPM>30 Acss2 expression

ss8 <- subset(ss,target_gene=="ACSS2" & construct=="kd")
xx8 <- xx[,colnames(xx) %in% ss8$Sample_ID ]

rpm8 <- apply(xx8,2, function(x) { x / sum(x) } ) *1000000
acss2 <- rpm6[grep("Acss2$",rownames(rpm6)),]
barplot(acss2 ,horiz=TRUE, las=1)

ss8$acss2 <- acss2

ss8 <-subset(ss8,acss2<30)
xx8 <- xx[,colnames(xx) %in% ss8$Sample_ID ]

dim(xx8)

## [1] 55357     8

xx8 <- xx8[which(rowMeans(xx8)>=10),]
dim(xx8)

## [1] 9287    8

ss8$ex_sed <- factor(ss8$ex_sed,levels=c("Sed","Ex"))

dds <- DESeqDataSetFromMatrix(countData = xx8 , colData = ss8, design = ~ ex_sed )

## converting counts to integer mode

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])

dge[1:20,1:6] %>%
  kbl(caption = "Top gene expression differences for contrast 5: Effect of exercise in ACSS2 kd mice") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of exercise in ACSS2 kd mice

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000025006.19 Sorbs1	109.58623	1.4457691	0.1775581	8.142511	0.00e+00	0.0000000
ENSMUSG00000041075.9 Fzd7	121.49670	-0.9761301	0.1859462	-5.249531	2.00e-07	0.0004878
ENSMUSG00000010064.16 Slc38a3	144.23642	0.9823922	0.1894790	5.184702	2.00e-07	0.0004878
ENSMUSG00000046818.8 Ddit4l	408.60248	-0.8703481	0.1727187	-5.039108	5.00e-07	0.0007909
ENSMUSG00000035621.14 Midn	63.78028	1.0317688	0.2138104	4.825625	1.40e-06	0.0018880
ENSMUSG00000028584.4 Lrrc38	66.61206	-1.1070210	0.2332827	-4.745404	2.10e-06	0.0023459
ENSMUSG00000024789.14 Jak2	193.47888	0.7514293	0.1628108	4.615353	3.90e-06	0.0034316
ENSMUSG00000068606.7 Gm4841	58.98670	-1.0982733	0.2383222	-4.608354	4.10e-06	0.0034316
ENSMUSG00000000628.11 Hk2	613.59872	0.6239488	0.1427484	4.370970	1.24e-05	0.0092964
ENSMUSG00000035963.9 Odf3l2	73.98266	1.2665599	0.2932095	4.319641	1.56e-05	0.0094680
ENSMUSG00000003810.14 Mast2	226.03654	0.7699018	0.1785989	4.310786	1.63e-05	0.0094680
ENSMUSG00000068699.13 Flnc	846.31933	0.8820455	0.2049507	4.303697	1.68e-05	0.0094680
ENSMUSG00000035105.6 Egln3	52.15483	1.2357813	0.2974319	4.154838	3.26e-05	0.0169370
ENSMUSG00000025612.6 Bach1	97.48204	0.8551240	0.2067572	4.135885	3.54e-05	0.0170834
ENSMUSG00000029167.14 Ppargc1a	224.72043	2.1684242	0.5271932	4.113149	3.90e-05	0.0172041
ENSMUSG00000113302.2 Gm32219	50.70732	1.3894107	0.3385921	4.103494	4.07e-05	0.0172041
ENSMUSG00000030852.19 Tacc2	817.34324	0.8484866	0.2084220	4.071003	4.68e-05	0.0186253
ENSMUSG00000026672.12 Optn	200.73678	0.8233917	0.2040656	4.034936	5.46e-05	0.0195031
ENSMUSG00000021903.12 Galnt15	41.28021	1.3592692	0.3369348	4.034220	5.48e-05	0.0195031
ENSMUSG00000043639.16 Rbm20	167.36912	0.8115622	0.2027473	4.002826	6.26e-05	0.0209901

dge8 <- dge
d8up <- rownames(subset(dge8,padj <= 0.05 & log2FoldChange > 0))
d8dn <- rownames(subset(dge8,padj <= 0.05 & log2FoldChange < 0))
write.table(dge8,file="dge8kd.tsv",quote=FALSE,sep="\t")

Conclusion

Some observations:

• Read counts were a bit low which translated to a smaller number of genes detected. Typically this is in the range of ~15000 but we detected ~9500 in the 4 contrasts considered here. This may impact the number of genes that are detected as differentially expressed.

• Acly (which encodes ATP-CL) was not measured at a sufficiently high level, and did not exhibit a decrease in KD mice.

• Acly knock-down did not result in any major gene expression changes in Sed or Ex mice.

• Acss2 was detected at a relatively low level, and it did exhibit a decrease in kd samples.

• Animal 6 might be an outlier in the Acss2 study (5176172).

• DESeq2 identified Acss2 as significantly downregulated in the ks mice but there were no other significant differences.

• In the list of top significant genes for Acss2 kd in Ex mice, there are some mitochondrial and metabolism genes. These are likely themes for pathway analysis.

• As the success of the knockdown was not consistent, we may need to omit animals or take a different approach to the expression analysis. For example by contructing the expression analysis by looking for genes that correlate with Acly and Acss2 instead of a case-control design.

Session information

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] beeswarm_0.4.0              vioplot_0.3.7              
##  [3] sm_2.2-5.7                  kableExtra_1.3.4           
##  [5] topconfects_1.12.0          limma_3.52.1               
##  [7] eulerr_6.1.1                mitch_1.8.0                
##  [9] MASS_7.3-58                 fgsea_1.22.0               
## [11] gplots_3.1.3                DESeq2_1.36.0              
## [13] SummarizedExperiment_1.26.1 Biobase_2.56.0             
## [15] MatrixGenerics_1.8.0        matrixStats_0.62.0         
## [17] GenomicRanges_1.48.0        GenomeInfoDb_1.32.2        
## [19] IRanges_2.30.0              S4Vectors_0.34.0           
## [21] BiocGenerics_0.42.0         reshape2_1.4.4             
## [23] forcats_0.5.1               stringr_1.4.0              
## [25] dplyr_1.0.9                 purrr_0.3.4                
## [27] readr_2.1.2                 tidyr_1.2.0                
## [29] tibble_3.1.7                ggplot2_3.3.6              
## [31] tidyverse_1.3.1             zoo_1.8-10                 
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3       ellipsis_0.3.2         XVector_0.36.0        
##   [4] fs_1.5.2               rstudioapi_0.13        bit64_4.0.5           
##   [7] AnnotationDbi_1.58.0   fansi_1.0.3            lubridate_1.8.0       
##  [10] xml2_1.3.3             codetools_0.2-18       splines_4.2.1         
##  [13] cachem_1.0.6           knitr_1.39             geneplotter_1.74.0    
##  [16] jsonlite_1.8.0         broom_0.8.0            annotate_1.74.0       
##  [19] dbplyr_2.2.1           png_0.1-7              shiny_1.7.1           
##  [22] compiler_4.2.1         httr_1.4.3             backports_1.4.1       
##  [25] assertthat_0.2.1       Matrix_1.4-1           fastmap_1.1.0         
##  [28] cli_3.3.0              later_1.3.0            htmltools_0.5.2       
##  [31] tools_4.2.1            gtable_0.3.0           glue_1.6.2            
##  [34] GenomeInfoDbData_1.2.8 fastmatch_1.1-3        Rcpp_1.0.8.3          
##  [37] jquerylib_0.1.4        cellranger_1.1.0       vctrs_0.4.1           
##  [40] Biostrings_2.64.0      svglite_2.1.0          xfun_0.31             
##  [43] rvest_1.0.2            mime_0.12              lifecycle_1.0.1       
##  [46] gtools_3.9.2.2         XML_3.99-0.10          zlibbioc_1.42.0       
##  [49] scales_1.2.0           hms_1.1.1              promises_1.2.0.1      
##  [52] parallel_4.2.1         RColorBrewer_1.1-3     yaml_2.3.5            
##  [55] memoise_2.0.1          gridExtra_2.3          sass_0.4.1            
##  [58] reshape_0.8.9          stringi_1.7.6          RSQLite_2.2.14        
##  [61] highr_0.9              genefilter_1.78.0      caTools_1.18.2        
##  [64] BiocParallel_1.30.3    echarts4r_0.4.4        systemfonts_1.0.4     
##  [67] rlang_1.0.3            pkgconfig_2.0.3        bitops_1.0-7          
##  [70] evaluate_0.15          lattice_0.20-45        htmlwidgets_1.5.4     
##  [73] bit_4.0.4              tidyselect_1.1.2       GGally_2.1.2          
##  [76] plyr_1.8.7             magrittr_2.0.3         R6_2.5.1              
##  [79] generics_0.1.2         DelayedArray_0.22.0    DBI_1.1.3             
##  [82] pillar_1.7.0           haven_2.5.0            withr_2.5.0           
##  [85] survival_3.4-0         KEGGREST_1.36.2        RCurl_1.98-1.7        
##  [88] modelr_0.1.8           crayon_1.5.1           KernSmooth_2.23-20    
##  [91] utf8_1.2.2             rmarkdown_2.14         tzdb_0.3.0            
##  [94] locfit_1.5-9.5         grid_4.2.1             readxl_1.4.0          
##  [97] data.table_1.14.2      blob_1.2.3             webshot_0.5.3         
## [100] reprex_2.0.1           digest_0.6.29          xtable_1.8-4          
## [103] httpuv_1.6.5           munsell_0.5.0          viridisLite_0.4.0     
## [106] bslib_0.3.1