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.

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 <- ll/colSums(ll) *1e6
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’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 <- xx_atp/colSums(xx_atp)*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 <- xx_acs/colSums(xx_acs)*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.

Effect of ATP-CL KD in Sed animals. (ss_atp_sed)
Effect of ATP-CL KD in Ex animals. (ss_atp_ex)
Effect of ACSS2 KD in Sed animals. (ss_acs_sed)
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")

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")

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")

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")

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.
The

Session information

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 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/liblapack.so.3
## 
## 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] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   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.8.0           limma_3.48.3               
##  [7] eulerr_6.1.1                mitch_1.5.1                
##  [9] MASS_7.3-57                 fgsea_1.18.0               
## [11] gplots_3.1.1                DESeq2_1.32.0              
## [13] SummarizedExperiment_1.22.0 Biobase_2.52.0             
## [15] MatrixGenerics_1.4.3        matrixStats_0.61.0         
## [17] GenomicRanges_1.44.0        GenomeInfoDb_1.28.4        
## [19] IRanges_2.26.0              S4Vectors_0.30.0           
## [21] BiocGenerics_0.38.0         reshape2_1.4.4             
## [23] forcats_0.5.1               stringr_1.4.0              
## [25] dplyr_1.0.8                 purrr_0.3.4                
## [27] readr_2.1.2                 tidyr_1.2.0                
## [29] tibble_3.1.6                ggplot2_3.3.5              
## [31] tidyverse_1.3.1             zoo_1.8-9                  
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3       ellipsis_0.3.2         XVector_0.32.0        
##   [4] fs_1.5.2               rstudioapi_0.13        bit64_4.0.5           
##   [7] AnnotationDbi_1.54.1   fansi_1.0.2            lubridate_1.8.0       
##  [10] xml2_1.3.3             splines_4.2.0          cachem_1.0.6          
##  [13] knitr_1.37             geneplotter_1.70.0     jsonlite_1.7.3        
##  [16] broom_0.7.12           annotate_1.70.0        dbplyr_2.1.1          
##  [19] png_0.1-7              shiny_1.7.1            compiler_4.2.0        
##  [22] httr_1.4.2             backports_1.4.1        assertthat_0.2.1      
##  [25] Matrix_1.4-1           fastmap_1.1.0          cli_3.2.0             
##  [28] later_1.3.0            htmltools_0.5.2        tools_4.2.0           
##  [31] gtable_0.3.0           glue_1.6.1             GenomeInfoDbData_1.2.6
##  [34] fastmatch_1.1-3        Rcpp_1.0.8             jquerylib_0.1.4       
##  [37] cellranger_1.1.0       vctrs_0.3.8            Biostrings_2.60.2     
##  [40] svglite_2.1.0          xfun_0.29              rvest_1.0.2           
##  [43] mime_0.12              lifecycle_1.0.1        gtools_3.9.2          
##  [46] XML_3.99-0.8           zlibbioc_1.38.0        scales_1.1.1          
##  [49] hms_1.1.1              promises_1.2.0.1       RColorBrewer_1.1-2    
##  [52] yaml_2.3.5             memoise_2.0.1          gridExtra_2.3         
##  [55] sass_0.4.0             reshape_0.8.8          stringi_1.7.6         
##  [58] RSQLite_2.2.10         highr_0.9              genefilter_1.74.0     
##  [61] caTools_1.18.2         BiocParallel_1.26.2    echarts4r_0.4.3       
##  [64] systemfonts_1.0.4      rlang_1.0.1            pkgconfig_2.0.3       
##  [67] bitops_1.0-7           evaluate_0.15          lattice_0.20-45       
##  [70] htmlwidgets_1.5.4      bit_4.0.4              tidyselect_1.1.2      
##  [73] GGally_2.1.2           plyr_1.8.6             magrittr_2.0.2        
##  [76] R6_2.5.1               generics_0.1.2         DelayedArray_0.18.0   
##  [79] DBI_1.1.2              pillar_1.7.0           haven_2.4.3           
##  [82] withr_2.4.3            survival_3.2-13        KEGGREST_1.32.0       
##  [85] RCurl_1.98-1.6         modelr_0.1.8           crayon_1.5.0          
##  [88] KernSmooth_2.23-20     utf8_1.2.2             rmarkdown_2.11        
##  [91] tzdb_0.2.0             locfit_1.5-9.4         grid_4.2.0            
##  [94] readxl_1.3.1           data.table_1.14.2      blob_1.2.2            
##  [97] webshot_0.5.2          reprex_2.0.1           digest_0.6.29         
## [100] xtable_1.8-4           httpuv_1.6.5           munsell_0.5.0         
## [103] viridisLite_0.4.0      bslib_0.3.1

Histone Acetyl RNA-seq DGE

Mark Ziemann

2022-05-30