Emily Crispr RNA-seq

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

Introduction

Here we analyse the effect of various CRISPR Cas13 constructs.

There are two different Cas13 sequences used, the normal RfxCas13d (WT) and the high-fidelity HfCas13d.

Sample	Cas13 type	Sample type
S1	Rfx	DsRed gRNA
S2	Rfx	DsRed gRNA
S3	Rfx	DsRed gRNA
S4	Rfx	Non-targeting gRNA
S5	Rfx	Non-targeting gRNA
S6	Rfx	Non-targeting gRNA
S7	Rfx	Transfection control
S8	Rfx	Transfection control
S9	Rfx	Transfection control
S10	Rfx	Non-transfected control
S11	Rfx	Non-transfected control
S12	Rfx	Non-transfected control
S13	Hf	DsRed gRNA
S14	Hf	DsRed gRNA
S15	Hf	DsRed gRNA
S16	Hf	Non-targeting gRNA
S17	Hf	Non-targeting gRNA
S18	Hf	Non-targeting gRNA
S19	Hf	Transfection control
S20	Hf	Transfection control
S21	Hf	Transfection control
S22	Hf	Non-transfected control
S23	Hf	Non-transfected control
S24	Hf	Non-transfected control

They are similar except for a few bp. The divergent bp were replaced with N for mapping with kallisto.

The goal is to understand the performance of each Cas, in terms of degree of knock-down and the degree of off-target degradation.

Methods

Reads were mapped to chicken transcripts from ensembl version 113.

Genes with mean counts > 10 are classified as detected.

Differential expression is conducted with DESeq2.

Pathway enrichment analysis is conducted with mitch.

suppressPackageStartupMessages({
    library("zoo")
    library("dplyr")
    library("reshape2")
    library("DESeq2")
    library("gplots")
    library("MASS")
    library("mitch")
    library("eulerr")
    library("kableExtra")
    library("beeswarm")
    library("UpSetR")
})

Sample sheet

ss <- read.table("samplesheet.tsv")
ss <- ss[order(rownames(ss)),]
ss

##        SampleGroup      SampleLabel
## S1  Rfx_DsRed_gRNA Rfx_DsRed_gRNA_1
## S10   Rfx_NTfx_ctl   Rfx_NTfx_ctl_1
## S11   Rfx_NTfx_ctl   Rfx_NTfx_ctl_2
## S12   Rfx_NTfx_ctl   Rfx_NTfx_ctl_3
## S13  Hf_DsRed_gRNA  Hf_DsRed_gRNA_1
## S14  Hf_DsRed_gRNA  Hf_DsRed_gRNA_2
## S15  Hf_DsRed_gRNA  Hf_DsRed_gRNA_3
## S16    Hf_NTC_gRNA    Hf_NTC_gRNA_1
## S17    Hf_NTC_gRNA    Hf_NTC_gRNA_2
## S18    Hf_NTC_gRNA    Hf_NTC_gRNA_3
## S19     Hf_Tfx_ctl     Hf_Tfx_ctl_1
## S2  Rfx_DsRed_gRNA Rfx_DsRed_gRNA_2
## S20     Hf_Tfx_ctl     Hf_Tfx_ctl_2
## S21     Hf_Tfx_ctl     Hf_Tfx_ctl_3
## S22    Hf_NTfx_ctl    Hf_NTfx_ctl_1
## S23    Hf_NTfx_ctl    Hf_NTfx_ctl_2
## S24    Hf_NTfx_ctl    Hf_NTfx_ctl_3
## S3  Rfx_DsRed_gRNA Rfx_DsRed_gRNA_3
## S4    Rfx_NTC_gRNA   Rfx_NTC_gRNA_1
## S5    Rfx_NTC_gRNA   Rfx_NTC_gRNA_2
## S6    Rfx_NTC_gRNA   Rfx_NTC_gRNA_3
## S7     Rfx_Tfx_ctl    Rfx_Tfx_ctl_1
## S8     Rfx_Tfx_ctl    Rfx_Tfx_ctl_2
## S9     Rfx_Tfx_ctl    Rfx_Tfx_ctl_3

Import gene ontology sets for pathway analysis later

go <- gmt_import("chicken_go_2024-11.gmt")

Import gene table and read counts

gt <- read.table("ref/Gallus_gallus.bGalGal1.mat.broiler.GRCg7b.113.tx2gene.tsv")
head(gt)

##                   V1                 V2    V3
## 1 ENSGALT00010001513 ENSGALG00010000711 SPRY2
## 2 ENSGALT00010001514 ENSGALG00010000711 SPRY2
## 3 ENSGALT00010001515 ENSGALG00010000711 SPRY2
## 4 ENSGALT00010001516 ENSGALG00010000711 SPRY2
## 5 ENSGALT00010001518 ENSGALG00010000711 SPRY2
## 6 ENSGALT00010001519 ENSGALG00010000711 SPRY2

rownames(gt) <- gt$V1
gt$gene <- paste(gt$V2,gt$V3)
gt$V1 = gt$V2 = gt$V3 = NULL
head(gt)

##                                        gene
## ENSGALT00010001513 ENSGALG00010000711 SPRY2
## ENSGALT00010001514 ENSGALG00010000711 SPRY2
## ENSGALT00010001515 ENSGALG00010000711 SPRY2
## ENSGALT00010001516 ENSGALG00010000711 SPRY2
## ENSGALT00010001518 ENSGALG00010000711 SPRY2
## ENSGALT00010001519 ENSGALG00010000711 SPRY2

# add cas13d and DsRed
add <- data.frame(c("cas13d cas13d","DsRed DsRed","Fbxw11-GFP Fbxw11-GFP"))
rownames(add) <- c("cas13d","DsRed","Fbxw11-GFP")
colnames(add) <- "gene"
gt2 <- as.data.frame(rbind(add,gt))
head(gt2)

##                                        gene
## cas13d                        cas13d cas13d
## DsRed                           DsRed DsRed
## Fbxw11-GFP            Fbxw11-GFP Fbxw11-GFP
## ENSGALT00010001513 ENSGALG00010000711 SPRY2
## ENSGALT00010001514 ENSGALG00010000711 SPRY2
## ENSGALT00010001515 ENSGALG00010000711 SPRY2

tmp <- read.table("AGRF_CAGRF24100018-1_22VGYMLT3/3col.tsv.gz",header=F)
tmp$V1 <-  sapply(strsplit(tmp$V1,"_"),"[[",1)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
rownames(x) <- sapply(strsplit( rownames(x) , "\\."),"[[",1)

xm <- merge(x,gt2,by=0)
rownames(xm) <- xm$Row.names
xm$Row.names = NULL
xmg <- aggregate(. ~ gene,xm,sum)
rownames(xmg) <- xmg$gene

rownames(xmg) <- xmg$gene
xmg$gene = NULL
xx <- round(xmg)
head(xx)

##                                          S1   S10   S11   S12   S13   S14   S15
## cas13d cas13d                         44039 37370 36900 40161 53106 48224 49161
## DsRed DsRed                            2095     0     0     0  3951  3910  3807
## ENSGALG00010000002 ENSGALG00010000002    17     7     9    16    15    17    17
## ENSGALG00010000003 ENSGALG00010000003 22387 19716 22466 19877 26183 25305 28846
## ENSGALG00010000004 ENSGALG00010000004     1     0     0     0     1     0     0
## ENSGALG00010000005 ENSGALG00010000005 18511 15071 19554 14291 16848 18041 20056
##                                         S16   S17   S18   S19    S2   S20   S21
## cas13d cas13d                         48976 44541 39904 46972 39665 44028 43130
## DsRed DsRed                           12568 10316  9196 10529  1855 10205 11838
## ENSGALG00010000002 ENSGALG00010000002    19    12    14    11    14    11    14
## ENSGALG00010000003 ENSGALG00010000003 28620 25906 21255 27990 17621 22773 23329
## ENSGALG00010000004 ENSGALG00010000004     1     1     0     0     0     2     0
## ENSGALG00010000005 ENSGALG00010000005 26851 23440 13559 19179 13201 17624 15785
##                                         S22   S23   S24    S3    S4    S5    S6
## cas13d cas13d                         38244 44665 39134 47778 43790 44436 39163
## DsRed DsRed                               0     0     0  2588 17544 15247 13850
## ENSGALG00010000002 ENSGALG00010000002    13    12    11    23    15     8    22
## ENSGALG00010000003 ENSGALG00010000003 24799 26109 26392 26643 25925 21893 23132
## ENSGALG00010000004 ENSGALG00010000004     0     0     1     1     1     0     0
## ENSGALG00010000005 ENSGALG00010000005 16355 17065 19702 20134 23215 19749 19385
##                                          S7    S8    S9
## cas13d cas13d                         47522 41347 46306
## DsRed DsRed                           19799 13067 18064
## ENSGALG00010000002 ENSGALG00010000002    20     9    13
## ENSGALG00010000003 ENSGALG00010000003 24793 21402 25391
## ENSGALG00010000004 ENSGALG00010000004     0     0     0
## ENSGALG00010000005 ENSGALG00010000005 16816 18267 15378

tail(xx)

##                             S1  S10  S11  S12  S13  S14  S15   S16   S17   S18
## ENSGALG00010030105 RAB40B   49   47   44   56   59   39   51    39    36    34
## ENSGALG00010030106 GPRC5C    0    1    4    1    5    2    4     2     1     4
## ENSGALG00010030107 RPTOR  1583 1434 1343 1352 1629 1415 1489  1619  1492  1381
## ENSGALG00010030108 LLGL2   996 1058  992 1030 1278 1291 1374  1383  1217  1038
## ENSGALG00010030109 FBF1    830  734  607  729  786  688  650   703   733   649
## Fbxw11-GFP Fbxw11-GFP     6811  307  308  303 6969 6301 6405 19976 15379 14478
##                             S19   S2   S20   S21  S22  S23  S24   S3    S4
## ENSGALG00010030105 RAB40B    49   36    43    45   42   52   33   52    53
## ENSGALG00010030106 GPRC5C     1    2     0     1    0    0    0    2     0
## ENSGALG00010030107 RPTOR   1462 1366  1432  1387 1207 1453 1260 1768  1756
## ENSGALG00010030108 LLGL2   1308 1036  1181  1210 1279 1456 1278 1253  1162
## ENSGALG00010030109 FBF1     773  741   682   626  617  726  575  853   879
## Fbxw11-GFP Fbxw11-GFP     10927 5319 10100 11594  311  306  309 7476 26922
##                              S5    S6    S7    S8    S9
## ENSGALG00010030105 RAB40B    46    42    55    51    40
## ENSGALG00010030106 GPRC5C     4     2     0     1     5
## ENSGALG00010030107 RPTOR   1609  1270  1697  1355  1586
## ENSGALG00010030108 LLGL2   1084   948  1100  1045  1242
## ENSGALG00010030109 FBF1     756   683   825   729   804
## Fbxw11-GFP Fbxw11-GFP     22776 20700 18592 12719 16902

# reordering the columns
rownames(ss) == colnames(xx)

##  [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

colnames(xx) <- ss$SampleLabel

xx <- xx[,order(colnames(xx))]
ss <- ss[order(ss$SampleLabel),]

# write counts
write.table(xx,"emily_counts.tsv",sep="\t",quote=FALSE)

# add colour and shape info to ss
ss$col <- c(rep("orange",12),rep("gray",12))
ss$shape <- rep(c(rep(17,3),rep(15,3),rep(18,3),rep(19,3)),2)

QC analysis

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

par(mar=c(5,10,3,1))
barplot(rev(colSums(xx)),horiz=TRUE,las=1,xlab="num reads")

colSums(xx)

##  Hf_DsRed_gRNA_1  Hf_DsRed_gRNA_2  Hf_DsRed_gRNA_3    Hf_NTC_gRNA_1 
##         36069938         33437002         35100639         35955443 
##    Hf_NTC_gRNA_2    Hf_NTC_gRNA_3    Hf_NTfx_ctl_1    Hf_NTfx_ctl_2 
##         33725699         29406778         29817493         34580661 
##    Hf_NTfx_ctl_3     Hf_Tfx_ctl_1     Hf_Tfx_ctl_2     Hf_Tfx_ctl_3 
##         30138577         32824476         31168254         30781130 
## Rfx_DsRed_gRNA_1 Rfx_DsRed_gRNA_2 Rfx_DsRed_gRNA_3   Rfx_NTC_gRNA_1 
##         33290779         30357975         36657286         35649182 
##   Rfx_NTC_gRNA_2   Rfx_NTC_gRNA_3   Rfx_NTfx_ctl_1   Rfx_NTfx_ctl_2 
##         35294526         30300316         31195104         30097571 
##   Rfx_NTfx_ctl_3    Rfx_Tfx_ctl_1    Rfx_Tfx_ctl_2    Rfx_Tfx_ctl_3 
##         31087737         37246761         31051428         36167680

par(mar = c(5.1, 4.1, 4.1, 2.1) )

MDS plot

Now check MDS clustering.

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

XMIN=min(mymds[,1])*1.3
XMAX=max(mymds[,1])*1.5

plot(mymds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=ss$shape, cex=4, col=ss$col,
  xlim=c(XMIN,XMAX), main = "MDS plot")

text(cmdscale(dist(t(xx))), labels=colnames(xx) )

legend("topleft", inset=.02, title="Cas13 type",
   c("Hf","Rfx"), fill=c("orange","gray"), horiz=TRUE, cex=1)

legend("topright", inset=.02, title="Sample group",
  c("DsRed_gRNA","NTC_gRNA","NTfx_ctl", "Tfx_ctl"),
  pch=c(17,15,18,19) , col="black", horiz=FALSE, cex=1)

It looks like there is clustering in terms of the Cas construct used but not the sample groups.

Let’s take a closer look at the different Cas types.

xrf <- xx[,grep("Rfx",colnames(xx))]
xhf <- xx[,grep("Hf",colnames(xx))]

srf <- ss[grep("Rfx",ss$SampleGroup),]
shf <- ss[grep("Hf",ss$SampleGroup),]

# Rfx Cas13d
mymds <- cmdscale(dist(t(xrf)))

XMIN=min(mymds[,1])*1.3
XMAX=max(mymds[,1])*1.5

plot(mymds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=srf$shape, cex=4, col=srf$col,
  xlim=c(XMIN,XMAX), main = "MDS plot: Rfx")

text(cmdscale(dist(t(xrf))), labels=colnames(xrf) )

legend("topleft", inset=.02, title="Sample group",
  c("DsRed_gRNA","NTC_gRNA","NTfx_ctl", "Tfx_ctl"),
  pch=c(17,15,18,19) , col="gray", horiz=FALSE, cex=1)

# Hf Cas13d
mymds <- cmdscale(dist(t(xhf)))

XMIN=min(mymds[,1])*1.3
XMAX=max(mymds[,1])*1.5

plot(mymds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=shf$shape, cex=4, col=shf$col,
  xlim=c(XMIN,XMAX), main = "MDS plot: Hf")

text(cmdscale(dist(t(xhf))), labels=colnames(xhf) )

legend("right", inset=.02, title="Sample group",
  c("DsRed_gRNA","NTC_gRNA","NTfx_ctl", "Tfx_ctl"),
  pch=c(17,15,18,19) , col="orange", horiz=FALSE, cex=1)

In the Hf experiment, non-transfected controls were clustered separately from the other samples. This was observed to a lesser extent in the Rfx samples. No other notable clustering was observed.

Correlation heatmap

mycor <- cor(xx)
mycor[mycor == 1] <- NA

heatmap.2(mycor,trace="n",main="Pearson correlation heatmap",mar=c(9,9))

mycor %>% kbl(caption = "Pearson correlation coefficients") %>% kable_paper("hover", full_width = F)

Pearson correlation coefficients

	Hf_DsRed_gRNA_1	Hf_DsRed_gRNA_2	Hf_DsRed_gRNA_3	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3	Hf_NTfx_ctl_1	Hf_NTfx_ctl_2	Hf_NTfx_ctl_3	Hf_Tfx_ctl_1	Hf_Tfx_ctl_2	Hf_Tfx_ctl_3	Rfx_DsRed_gRNA_1	Rfx_DsRed_gRNA_2	Rfx_DsRed_gRNA_3	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3	Rfx_NTfx_ctl_1	Rfx_NTfx_ctl_2	Rfx_NTfx_ctl_3	Rfx_Tfx_ctl_1	Rfx_Tfx_ctl_2	Rfx_Tfx_ctl_3
Hf_DsRed_gRNA_1	NA	0.9959796	0.9974034	0.9974165	0.9971398	0.9992745	0.9910852	0.9926766	0.9915194	0.9972227	0.9992909	0.9988950	0.9621010	0.9702130	0.9721555	0.9677693	0.9780049	0.9804914	0.9741073	0.9842743	0.9919543	0.9718543	0.9655399	0.9739115
Hf_DsRed_gRNA_2	0.9959796	NA	0.9989044	0.9930099	0.9982336	0.9965715	0.9948011	0.9953406	0.9957425	0.9989539	0.9964215	0.9957073	0.9434608	0.9542662	0.9549634	0.9530961	0.9626486	0.9693555	0.9597702	0.9769123	0.9825133	0.9561091	0.9498399	0.9595844
Hf_DsRed_gRNA_3	0.9974034	0.9989044	NA	0.9953891	0.9981515	0.9977535	0.9946609	0.9951077	0.9954056	0.9984339	0.9979725	0.9980730	0.9516428	0.9612926	0.9626762	0.9598735	0.9694568	0.9748110	0.9663345	0.9803817	0.9872154	0.9634640	0.9571040	0.9664875
Hf_NTC_gRNA_1	0.9974165	0.9930099	0.9953891	NA	0.9970908	0.9964691	0.9918128	0.9919882	0.9925493	0.9955035	0.9980350	0.9982397	0.9717448	0.9774876	0.9795271	0.9776647	0.9855782	0.9874559	0.9822617	0.9902062	0.9948006	0.9797484	0.9750126	0.9814711
Hf_NTC_gRNA_2	0.9971398	0.9982336	0.9981515	0.9970908	NA	0.9977274	0.9936839	0.9938413	0.9950531	0.9984874	0.9979693	0.9976965	0.9535992	0.9615955	0.9634180	0.9618815	0.9715831	0.9756109	0.9676774	0.9815698	0.9870500	0.9640808	0.9582332	0.9668464
Hf_NTC_gRNA_3	0.9992745	0.9965715	0.9977535	0.9964691	0.9977274	NA	0.9904103	0.9921692	0.9913427	0.9972773	0.9992909	0.9990666	0.9551461	0.9637210	0.9662539	0.9615645	0.9731976	0.9755989	0.9688494	0.9802242	0.9893228	0.9661077	0.9591316	0.9682763
Hf_NTfx_ctl_1	0.9910852	0.9948011	0.9946609	0.9918128	0.9936839	0.9904103	NA	0.9994242	0.9994985	0.9966429	0.9932834	0.9919778	0.9540719	0.9656453	0.9631534	0.9653816	0.9685619	0.9783961	0.9714295	0.9870410	0.9862911	0.9663478	0.9633857	0.9699769
Hf_NTfx_ctl_2	0.9926766	0.9953406	0.9951077	0.9919882	0.9938413	0.9921692	0.9994242	NA	0.9988137	0.9972577	0.9946409	0.9929909	0.9533450	0.9655920	0.9632437	0.9642140	0.9683798	0.9777484	0.9711220	0.9864766	0.9869968	0.9661606	0.9627263	0.9695830
Hf_NTfx_ctl_3	0.9915194	0.9957425	0.9954056	0.9925493	0.9950531	0.9913427	0.9994985	0.9988137	NA	0.9971487	0.9937722	0.9928073	0.9521592	0.9631497	0.9614631	0.9634359	0.9674701	0.9768346	0.9696447	0.9856247	0.9856144	0.9643072	0.9610845	0.9679984
Hf_Tfx_ctl_1	0.9972227	0.9989539	0.9984339	0.9955035	0.9984874	0.9972773	0.9966429	0.9972577	0.9971487	NA	0.9980671	0.9969898	0.9519753	0.9624284	0.9627117	0.9614057	0.9692855	0.9759783	0.9676205	0.9830046	0.9868791	0.9639302	0.9587466	0.9671196
Hf_Tfx_ctl_2	0.9992909	0.9964215	0.9979725	0.9980350	0.9979693	0.9992909	0.9932834	0.9946409	0.9937722	0.9980671	NA	0.9994164	0.9611362	0.9697521	0.9715769	0.9680148	0.9775598	0.9807668	0.9746758	0.9854452	0.9921599	0.9717292	0.9659574	0.9738680
Hf_Tfx_ctl_3	0.9988950	0.9957073	0.9980730	0.9982397	0.9976965	0.9990666	0.9919778	0.9929909	0.9928073	0.9969898	0.9994164	NA	0.9620190	0.9695946	0.9718051	0.9682352	0.9783244	0.9808679	0.9747173	0.9844215	0.9923434	0.9722322	0.9660124	0.9743267
Rfx_DsRed_gRNA_1	0.9621010	0.9434608	0.9516428	0.9717448	0.9535992	0.9551461	0.9540719	0.9533450	0.9521592	0.9519753	0.9611362	0.9620190	NA	0.9980555	0.9977787	0.9981284	0.9967016	0.9945627	0.9964082	0.9877442	0.9855619	0.9982592	0.9980322	0.9975290
Rfx_DsRed_gRNA_2	0.9702130	0.9542662	0.9612926	0.9774876	0.9615955	0.9637210	0.9656453	0.9655920	0.9631497	0.9624284	0.9697521	0.9695946	0.9980555	NA	0.9984265	0.9987679	0.9969579	0.9976436	0.9985534	0.9933774	0.9906351	0.9993970	0.9992563	0.9992905
Rfx_DsRed_gRNA_3	0.9721555	0.9549634	0.9626762	0.9795271	0.9634180	0.9662539	0.9631534	0.9632437	0.9614631	0.9627117	0.9715769	0.9718051	0.9977787	0.9984265	NA	0.9973730	0.9982602	0.9966610	0.9974699	0.9915804	0.9912411	0.9986829	0.9975959	0.9981605
Rfx_NTC_gRNA_1	0.9677693	0.9530961	0.9598735	0.9776647	0.9618815	0.9615645	0.9653816	0.9642140	0.9634359	0.9614057	0.9680148	0.9682352	0.9981284	0.9987679	0.9973730	NA	0.9971273	0.9978486	0.9985110	0.9935351	0.9891824	0.9989989	0.9993494	0.9989266
Rfx_NTC_gRNA_2	0.9780049	0.9626486	0.9694568	0.9855782	0.9715831	0.9731976	0.9685619	0.9683798	0.9674701	0.9692855	0.9775598	0.9783244	0.9967016	0.9969579	0.9982602	0.9971273	NA	0.9973448	0.9976182	0.9928745	0.9941695	0.9981635	0.9960024	0.9976632
Rfx_NTC_gRNA_3	0.9804914	0.9693555	0.9748110	0.9874559	0.9756109	0.9755989	0.9783961	0.9777484	0.9768346	0.9759783	0.9807668	0.9808679	0.9945627	0.9976436	0.9966610	0.9978486	0.9973448	NA	0.9981605	0.9978130	0.9954610	0.9982267	0.9972683	0.9989806
Rfx_NTfx_ctl_1	0.9741073	0.9597702	0.9663345	0.9822617	0.9676774	0.9688494	0.9714295	0.9711220	0.9696447	0.9676205	0.9746758	0.9747173	0.9964082	0.9985534	0.9974699	0.9985110	0.9976182	0.9981605	NA	0.9959400	0.9937477	0.9988922	0.9985741	0.9987924
Rfx_NTfx_ctl_2	0.9842743	0.9769123	0.9803817	0.9902062	0.9815698	0.9802242	0.9870410	0.9864766	0.9856247	0.9830046	0.9854452	0.9844215	0.9877442	0.9933774	0.9915804	0.9935351	0.9928745	0.9978130	0.9959400	NA	0.9961729	0.9935238	0.9929670	0.9949066
Rfx_NTfx_ctl_3	0.9919543	0.9825133	0.9872154	0.9948006	0.9870500	0.9893228	0.9862911	0.9869968	0.9856144	0.9868791	0.9921599	0.9923434	0.9855619	0.9906351	0.9912411	0.9891824	0.9941695	0.9954610	0.9937477	0.9961729	NA	0.9919964	0.9884013	0.9928792
Rfx_Tfx_ctl_1	0.9718543	0.9561091	0.9634640	0.9797484	0.9640808	0.9661077	0.9663478	0.9661606	0.9643072	0.9639302	0.9717292	0.9722322	0.9982592	0.9993970	0.9986829	0.9989989	0.9981635	0.9982267	0.9988922	0.9935238	0.9919964	NA	0.9991751	0.9997176
Rfx_Tfx_ctl_2	0.9655399	0.9498399	0.9571040	0.9750126	0.9582332	0.9591316	0.9633857	0.9627263	0.9610845	0.9587466	0.9659574	0.9660124	0.9980322	0.9992563	0.9975959	0.9993494	0.9960024	0.9972683	0.9985741	0.9929670	0.9884013	0.9991751	NA	0.9990639
Rfx_Tfx_ctl_3	0.9739115	0.9595844	0.9664875	0.9814711	0.9668464	0.9682763	0.9699769	0.9695830	0.9679984	0.9671196	0.9738680	0.9743267	0.9975290	0.9992905	0.9981605	0.9989266	0.9976632	0.9989806	0.9987924	0.9949066	0.9928792	0.9997176	0.9990639	NA

mycor <- cor(xx,method="spearman")
mycor[mycor == 1] <- NA
heatmap.2(mycor,trace="n",main="Spearman correlation heatmap",mar=c(9,9))

mycor %>% kbl(caption = "Spearman correlation coefficients") %>% kable_paper("hover", full_width = F)

Spearman correlation coefficients

	Hf_DsRed_gRNA_1	Hf_DsRed_gRNA_2	Hf_DsRed_gRNA_3	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3	Hf_NTfx_ctl_1	Hf_NTfx_ctl_2	Hf_NTfx_ctl_3	Hf_Tfx_ctl_1	Hf_Tfx_ctl_2	Hf_Tfx_ctl_3	Rfx_DsRed_gRNA_1	Rfx_DsRed_gRNA_2	Rfx_DsRed_gRNA_3	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3	Rfx_NTfx_ctl_1	Rfx_NTfx_ctl_2	Rfx_NTfx_ctl_3	Rfx_Tfx_ctl_1	Rfx_Tfx_ctl_2	Rfx_Tfx_ctl_3
Hf_DsRed_gRNA_1	NA	0.9599488	0.9593786	0.9602960	0.9592729	0.9584685	0.9579201	0.9601236	0.9577979	0.9601253	0.9596824	0.9588043	0.9524232	0.9537345	0.9536157	0.9549469	0.9536715	0.9525762	0.9520901	0.9510956	0.9521125	0.9537846	0.9531629	0.9542380
Hf_DsRed_gRNA_2	0.9599488	NA	0.9599999	0.9603373	0.9594227	0.9583982	0.9576127	0.9593367	0.9579286	0.9600871	0.9591133	0.9585310	0.9516034	0.9535027	0.9534788	0.9538586	0.9528460	0.9513969	0.9517635	0.9505957	0.9517369	0.9531870	0.9522520	0.9525369
Hf_DsRed_gRNA_3	0.9593786	0.9599999	NA	0.9603922	0.9596324	0.9591342	0.9582642	0.9591225	0.9582413	0.9599239	0.9597260	0.9595662	0.9526525	0.9540027	0.9531987	0.9545339	0.9535056	0.9520840	0.9528068	0.9520807	0.9518048	0.9531060	0.9514795	0.9535028
Hf_NTC_gRNA_1	0.9602960	0.9603373	0.9603922	NA	0.9614023	0.9589670	0.9574613	0.9596119	0.9594718	0.9608933	0.9603629	0.9603872	0.9534564	0.9529810	0.9533007	0.9562377	0.9554936	0.9523943	0.9526207	0.9513095	0.9520189	0.9542246	0.9537407	0.9541197
Hf_NTC_gRNA_2	0.9592729	0.9594227	0.9596324	0.9614023	NA	0.9577475	0.9576146	0.9586665	0.9579804	0.9588384	0.9589719	0.9591492	0.9527730	0.9529903	0.9533203	0.9550464	0.9544127	0.9520973	0.9527873	0.9511442	0.9518258	0.9541108	0.9522946	0.9531277
Hf_NTC_gRNA_3	0.9584685	0.9583982	0.9591342	0.9589670	0.9577475	NA	0.9564750	0.9585778	0.9572894	0.9585848	0.9576951	0.9586439	0.9506883	0.9513823	0.9517152	0.9529882	0.9528932	0.9513653	0.9511396	0.9499968	0.9505315	0.9515986	0.9519701	0.9518965
Hf_NTfx_ctl_1	0.9579201	0.9576127	0.9582642	0.9574613	0.9576146	0.9564750	NA	0.9580466	0.9578096	0.9579749	0.9572939	0.9575349	0.9494938	0.9506774	0.9507288	0.9525334	0.9505435	0.9495149	0.9513824	0.9496961	0.9503651	0.9517540	0.9503461	0.9518573
Hf_NTfx_ctl_2	0.9601236	0.9593367	0.9591225	0.9596119	0.9586665	0.9585778	0.9580466	NA	0.9591172	0.9602015	0.9588828	0.9585599	0.9507619	0.9523281	0.9511893	0.9535314	0.9525084	0.9504720	0.9519595	0.9508424	0.9514829	0.9517321	0.9516028	0.9520844
Hf_NTfx_ctl_3	0.9577979	0.9579286	0.9582413	0.9594718	0.9579804	0.9572894	0.9578096	0.9591172	NA	0.9581197	0.9575512	0.9578390	0.9502895	0.9516688	0.9512851	0.9535278	0.9515358	0.9499636	0.9517999	0.9509422	0.9509854	0.9510846	0.9510679	0.9519187
Hf_Tfx_ctl_1	0.9601253	0.9600871	0.9599239	0.9608933	0.9588384	0.9585848	0.9579749	0.9602015	0.9581197	NA	0.9599586	0.9593630	0.9524224	0.9532343	0.9534931	0.9547496	0.9537808	0.9520568	0.9520132	0.9512615	0.9522524	0.9535123	0.9527909	0.9533878
Hf_Tfx_ctl_2	0.9596824	0.9591133	0.9597260	0.9603629	0.9589719	0.9576951	0.9572939	0.9588828	0.9575512	0.9599586	NA	0.9593885	0.9527672	0.9530480	0.9530098	0.9542506	0.9530847	0.9516721	0.9524944	0.9511770	0.9519473	0.9534200	0.9519186	0.9529571
Hf_Tfx_ctl_3	0.9588043	0.9585310	0.9595662	0.9603872	0.9591492	0.9586439	0.9575349	0.9585599	0.9578390	0.9593630	0.9593885	NA	0.9526415	0.9529181	0.9527533	0.9543007	0.9530015	0.9513276	0.9525610	0.9502190	0.9513122	0.9533165	0.9515892	0.9531111
Rfx_DsRed_gRNA_1	0.9524232	0.9516034	0.9526525	0.9534564	0.9527730	0.9506883	0.9494938	0.9507619	0.9502895	0.9524224	0.9527672	0.9526415	NA	0.9568038	0.9573636	0.9597848	0.9588997	0.9555633	0.9562825	0.9544710	0.9549874	0.9585131	0.9568284	0.9585649
Rfx_DsRed_gRNA_2	0.9537345	0.9535027	0.9540027	0.9529810	0.9529903	0.9513823	0.9506774	0.9523281	0.9516688	0.9532343	0.9530480	0.9529181	0.9568038	NA	0.9577065	0.9593260	0.9580147	0.9557793	0.9562904	0.9551574	0.9563274	0.9577477	0.9561335	0.9577431
Rfx_DsRed_gRNA_3	0.9536157	0.9534788	0.9531987	0.9533007	0.9533203	0.9517152	0.9507288	0.9511893	0.9512851	0.9534931	0.9530098	0.9527533	0.9573636	0.9577065	NA	0.9591033	0.9575030	0.9562844	0.9566999	0.9555239	0.9555075	0.9583868	0.9563563	0.9579214
Rfx_NTC_gRNA_1	0.9549469	0.9538586	0.9545339	0.9562377	0.9550464	0.9529882	0.9525334	0.9535314	0.9535278	0.9547496	0.9542506	0.9543007	0.9597848	0.9593260	0.9591033	NA	0.9604467	0.9583005	0.9588380	0.9571838	0.9580201	0.9599535	0.9586112	0.9607349
Rfx_NTC_gRNA_2	0.9536715	0.9528460	0.9535056	0.9554936	0.9544127	0.9528932	0.9505435	0.9525084	0.9515358	0.9537808	0.9530847	0.9530015	0.9588997	0.9580147	0.9575030	0.9604467	NA	0.9561806	0.9576308	0.9557753	0.9568854	0.9589981	0.9578323	0.9590635
Rfx_NTC_gRNA_3	0.9525762	0.9513969	0.9520840	0.9523943	0.9520973	0.9513653	0.9495149	0.9504720	0.9499636	0.9520568	0.9516721	0.9513276	0.9555633	0.9557793	0.9562844	0.9583005	0.9561806	NA	0.9558707	0.9542115	0.9543574	0.9573377	0.9558973	0.9580704
Rfx_NTfx_ctl_1	0.9520901	0.9517635	0.9528068	0.9526207	0.9527873	0.9511396	0.9513824	0.9519595	0.9517999	0.9520132	0.9524944	0.9525610	0.9562825	0.9562904	0.9566999	0.9588380	0.9576308	0.9558707	NA	0.9554050	0.9567137	0.9570112	0.9564397	0.9572513
Rfx_NTfx_ctl_2	0.9510956	0.9505957	0.9520807	0.9513095	0.9511442	0.9499968	0.9496961	0.9508424	0.9509422	0.9512615	0.9511770	0.9502190	0.9544710	0.9551574	0.9555239	0.9571838	0.9557753	0.9542115	0.9554050	NA	0.9551307	0.9562743	0.9550974	0.9552916
Rfx_NTfx_ctl_3	0.9521125	0.9517369	0.9518048	0.9520189	0.9518258	0.9505315	0.9503651	0.9514829	0.9509854	0.9522524	0.9519473	0.9513122	0.9549874	0.9563274	0.9555075	0.9580201	0.9568854	0.9543574	0.9567137	0.9551307	NA	0.9583450	0.9554528	0.9567224
Rfx_Tfx_ctl_1	0.9537846	0.9531870	0.9531060	0.9542246	0.9541108	0.9515986	0.9517540	0.9517321	0.9510846	0.9535123	0.9534200	0.9533165	0.9585131	0.9577477	0.9583868	0.9599535	0.9589981	0.9573377	0.9570112	0.9562743	0.9583450	NA	0.9584630	0.9592594
Rfx_Tfx_ctl_2	0.9531629	0.9522520	0.9514795	0.9537407	0.9522946	0.9519701	0.9503461	0.9516028	0.9510679	0.9527909	0.9519186	0.9515892	0.9568284	0.9561335	0.9563563	0.9586112	0.9578323	0.9558973	0.9564397	0.9550974	0.9554528	0.9584630	NA	0.9577830
Rfx_Tfx_ctl_3	0.9542380	0.9525369	0.9535028	0.9541197	0.9531277	0.9518965	0.9518573	0.9520844	0.9519187	0.9533878	0.9529571	0.9531111	0.9585649	0.9577431	0.9579214	0.9607349	0.9590635	0.9580704	0.9572513	0.9552916	0.9567224	0.9592594	0.9577830	NA

# Rfx
mycor <- cor(xrf)
mycor[mycor == 1] <- NA
heatmap.2(mycor,trace="n",main="Rfx Pearson Cor",mar=c(11,11))

mycor <- cor(xrf,method="spearman")
mycor[mycor == 1] <- NA
heatmap.2(mycor,trace="n",main="Rfx Spearman Cor",mar=c(11,11))

# Hf
mycor <- cor(xhf)
mycor[mycor == 1] <- NA
heatmap.2(mycor,trace="n",main="Hf Pearson Cor",mar=c(11,11))

mycor <- cor(xhf,method="spearman")
mycor[mycor == 1] <- NA
heatmap.2(mycor,trace="n",main="Hf Spearman Cor",mar=c(11,11))

Now look specifically at the DsRed and Cas13 expression

rpm <- apply(xx,2, function(x) {x / sum(x) * 1e6 } )

cas13 <- unlist(rpm[1,,drop=TRUE])
dsred <- unlist(rpm[2,,drop=TRUE])
fbox <- rpm[nrow(rpm),,drop=TRUE]
par(mar=c(5,10,3,1))
barplot(cas13,horiz=TRUE,las=1,xlab="Reads per million",main="Cas13d")

barplot(dsred,horiz=TRUE,las=1,xlab="Reads per million",main="DsRed")

barplot(fbox,horiz=TRUE,las=1,xlab="Reads per million",main="Fbxw11-GFP")

par(mar = c(5.1, 4.1, 4.1, 2.1) )

There is a >1200 RPM for Cas13 in all these samples. DsRed was absent in NTfx samples. There was much less DsRed in targeting gDNA samples as compared to non-targeting controls and transfection controls. Fbxw11-GFP expression tracked DsRed.

## HF
ds <- dsred[grep("Hf",names(dsred))]

targeting_gRNA <- ds[grep("DsRed_gRNA",names(ds))]

nontargeting_gRNA <- ds[grep("NTC_gRNA",names(ds))]

dat <- list("Non-targeting gRNA"=nontargeting_gRNA,
"Targeting gRNA"=targeting_gRNA)

MAX=max(unlist(lapply(dat,max)))

boxplot(dat,ylab="Reads per million",ylim=c(0,MAX),
  main="DsRed expresssion (Hf)", col="white",frame=FALSE)
beeswarm(dat,add=TRUE,col="red",cex=2,pch=19)

dmean <- unlist(lapply(dat,mean))
reduction <- as.character(signif(unname(1-dmean[2]/dmean[1]),3)*100)
mtext(paste("Reduction:",reduction,"%"))

## RFX
ds <- dsred[grep("Rfx",names(dsred))]

targeting_gRNA <- ds[grep("DsRed_gRNA",names(ds))]

nontargeting_gRNA <- ds[grep("NTC_gRNA",names(ds))]

dat <- list("Non-targeting gRNA"=nontargeting_gRNA,
"Targeting gRNA"=targeting_gRNA)

MAX=max(unlist(lapply(dat,max)))

boxplot(dat,ylab="Reads per million",ylim=c(0,MAX),
  main="DsRed expresssion (Rfx)", col="white",frame=FALSE)
beeswarm(dat,add=TRUE,col="red",cex=2,pch=19)

dmean <- unlist(lapply(dat,mean))
reduction <- as.character(signif(unname(1-dmean[2]/dmean[1]),3)*100)
mtext(paste("Reduction:",reduction,"%"))

## HF
ds <- fbox[grep("Hf",names(fbox))]

targeting_gRNA <- ds[grep("DsRed_gRNA",names(ds))]

nontargeting_gRNA <- ds[grep("NTC_gRNA",names(ds))]

dat <- list("Non-targeting gRNA"=nontargeting_gRNA,
"Targeting gRNA"=targeting_gRNA)

MAX=max(unlist(lapply(dat,max)))

boxplot(dat,ylab="Reads per million",ylim=c(0,MAX),
  main="Fbxw11-GFP expresssion (Hf)", col="white",frame=FALSE)
beeswarm(dat,add=TRUE,col="red",cex=2,pch=19)

dmean <- unlist(lapply(dat,mean))
reduction <- as.character(signif(unname(1-dmean[2]/dmean[1]),3)*100)
mtext(paste("Reduction:",reduction,"%"))

## RFX
ds <- fbox[grep("Rfx",names(fbox))]

targeting_gRNA <- ds[grep("DsRed_gRNA",names(ds))]

nontargeting_gRNA <- ds[grep("NTC_gRNA",names(ds))]

dat <- list("Non-targeting gRNA"=nontargeting_gRNA,
"Targeting gRNA"=targeting_gRNA)

MAX=max(unlist(lapply(dat,max)))

boxplot(dat,ylab="Reads per million",ylim=c(0,MAX),
  main="Fbxw11-GFP expresssion (Rfx)", col="white",frame=FALSE)
beeswarm(dat,add=TRUE,col="red",cex=2,pch=19)

dmean <- unlist(lapply(dat,mean))
reduction <- as.character(signif(unname(1-dmean[2]/dmean[1]),3)*100)
mtext(paste("Reduction:",reduction,"%"))

Analysis of differential gene expression

I will be looking at the following comparisons (control vs case).

1. Rfx_NTC_gRNA versus Rfx_DsRed_gRNA

2. Rfx_Tfx_ctl versus Rfx_NTC_gRNA

3. Rfx_NTfc_ctl versus Rfx_Tfx_ctl

4. Hf_NTC_gRNA versus Hf_DsRed_gRNA

5. Hf_Tfx_ctl versus Hf_NTC_gRNA

6. Hf_NTfx_ctl versus Hf_Tfx_ctl

7. Rfx_DsRed_gRNA versus Hf_DsRed_gRNA

8. Rfx_NTC_gRNA versus Hf_NTC_gRNA

9. Rfx_Tfx_ctl versus Hf_Tfx_ctl

10. Rfx_NTfx_ctl versus Hf_NTfx_ctl

Plot functions

maplot <- function(de,contrast_name) {
  sig <-subset(de, padj < 0.05 )
  up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
  dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
  GENESUP <- length(up)
  GENESDN <- length(dn)
  DET=nrow(de)
  SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
  ns <-subset(de, padj > 0.05 )
  plot(log10(de$baseMean),de$log2FoldChange,
       xlab="log10 basemean", ylab="log2 foldchange",
       pch=19, cex=0.5, col="dark gray",
       main=contrast_name, cex.main=1)
  points(log10(sig$baseMean),sig$log2FoldChange,
         pch=19, cex=0.5, col="darkorange")
  mtext(SUBHEADER,cex = 1)
  gfprow <- grep("Fbxw11-GFP",rownames(de))
  points(log10(de$baseMean[gfprow]),de$log2FoldChange[gfprow],pch=19, cex=0.9, col="darkgreen")
  redrow <- grep("DsRed",rownames(de))
  points(log10(de$baseMean[redrow]),de$log2FoldChange[redrow],pch=19, cex=0.9, col="red")
}

make_volcano <- function(de,name) {
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$pval),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval", xlim=c(-6,6))
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$pval),cex=0.5,pch=19,col="darkorange")
    gfprow <- grep("Fbxw11-GFP",rownames(de))
    points(de$log2FoldChange[gfprow], -log10(de$pval)[gfprow],pch=19, cex=0.9, col="darkgreen")
    redrow <- grep("DsRed",rownames(de))
    points(de$log2FoldChange[redrow], -log10(de$pval)[redrow], pch=19, cex=0.9, col="red")
}

make_heatmap <- function(de,name,myss,mx,n=30){
  colfunc <- colorRampPalette(c("blue", "white", "red"))
  values <- myss$quickdash
  f <- colorRamp(c("yellow", "orange"))
  rr <- range(values)
  svals <- (values-rr[1])/diff(rr)
  colcols <- rgb(f(svals)/255)
  mxn <- mx/rowSums(mx)*1000000
  x <- mxn[which(rownames(mxn) %in% rownames(head(de,n))),]
  heatmap.2(as.matrix(x),trace="none",col=colfunc(25),scale="row",
    margins = c(7,20), cexRow=0.8, cexCol=0.8,
    main=paste("Top ranked",n,"genes in",name) )
}

1. Rfx_NTC_gRNA versus Rfx_DsRed_gRNA

ss1 <- subset(ss,SampleGroup=="Rfx_NTC_gRNA" | SampleGroup=="Rfx_DsRed_gRNA")
ss1$trt <- factor(grepl("DsRed",ss1$SampleGroup))
rownames(ss1) <- ss1$SampleLabel

xx1 <- xx[,colnames(xx) %in% ss1$SampleLabel]
xx1f <- xx1[rowMeans(xx1)>=10,]
dim(xx1)

## [1] 30111     6

dim(xx1f)

## [1] 14562     6

dds <- DESeqDataSetFromMatrix(countData = xx1f , colData = ss1, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by DsRed targeting in Rfx") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by DsRed targeting in Rfx

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Rfx_DsRed_gRNA_1	Rfx_DsRed_gRNA_2	Rfx_DsRed_gRNA_3	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3
DsRed DsRed	8779.4963	-2.8242664	0.0588480	-47.992526	0.0e+00	0.0000000	11.596088	11.553157	11.717575	14.095807	13.917811	13.997669
Fbxw11-GFP Fbxw11-GFP	14863.0158	-1.8337024	0.1030112	-17.801005	0.0e+00	0.0000000	12.948984	12.739591	12.948702	14.683145	14.464181	14.546489
ENSGALG00010023206 HSPB9	320.1663	-1.5530514	0.1390853	-11.166177	0.0e+00	0.0000000	9.855633	9.858837	9.856153	10.503028	10.297007	10.346312
ENSGALG00010005131 IL8L2	532.7586	0.7476336	0.1009870	7.403265	0.0e+00	0.0000000	10.611339	10.533044	10.669148	10.286743	10.235363	10.311760
ENSGALG00010016470 TNIP1	6411.7124	0.3313912	0.0467236	7.092582	0.0e+00	0.0000000	13.011005	13.006664	12.979349	12.767606	12.633823	12.733295
ENSGALG00010006134 K123	894.9555	0.4299283	0.0725414	5.926665	0.0e+00	0.0000055	10.953930	10.908059	10.939287	10.693424	10.714320	10.698324
ENSGALG00010025320 JAM3	5037.3056	-0.2966828	0.0535941	-5.535741	0.0e+00	0.0000470	12.348382	12.524727	12.434045	12.694746	12.716556	12.643131
ENSGALG00010019621 ENSGALG00010019621	13358.2914	0.2232910	0.0409066	5.458555	0.0e+00	0.0000636	13.928204	13.853888	13.961122	13.733939	13.676732	13.713168
ENSGALG00010001230 PTX3	1844.2211	0.3064948	0.0592945	5.169028	2.0e-07	0.0002773	11.530068	11.547733	11.602225	11.370614	11.304686	11.383077
ENSGALG00010015706 ENSGALG00010015706	2265.0717	0.2696937	0.0525580	5.131352	3.0e-07	0.0002984	11.740315	11.761993	11.778766	11.561940	11.541583	11.600211
ENSGALG00010025318 NUAK2	2131.7627	0.3066354	0.0599179	5.117593	3.0e-07	0.0002984	11.787872	11.676725	11.659855	11.496139	11.516273	11.466256
ENSGALG00010025057 SCP2	10658.6467	-0.2096236	0.0415755	-5.042001	5.0e-07	0.0004072	13.348541	13.453850	13.445804	13.618029	13.589614	13.616202
ENSGALG00010001754 APCDD1	13637.6828	-0.2158251	0.0433538	-4.978229	6.0e-07	0.0005236	13.653383	13.786873	13.774085	13.946430	13.951739	13.921061
ENSGALG00010029945 CCL4	18385.8762	0.2550680	0.0518998	4.914623	9.0e-07	0.0006740	14.313709	14.305542	14.463645	14.101595	14.049862	14.209415
ENSGALG00010022322 ENSGALG00010022322	3591.8592	0.3367394	0.0704462	4.780092	1.8e-06	0.0012390	12.171063	12.276865	12.427262	12.002925	11.963923	12.120157
ENSGALG00010019277 IGFBP5	16022.8418	0.1989457	0.0425197	4.678907	2.9e-06	0.0017365	14.128220	14.091270	14.228639	13.978997	13.943354	13.969404
ENSGALG00010020013 ADAM8	5047.5451	0.2162526	0.0462572	4.675008	2.9e-06	0.0017365	12.668944	12.597713	12.700948	12.447013	12.478143	12.502689
ENSGALG00010002340 RAB23	9745.1303	-0.1905774	0.0407697	-4.674489	2.9e-06	0.0017365	13.286394	13.359703	13.281169	13.478945	13.513220	13.452068
ENSGALG00010009432 BOP1	6378.4297	0.2138124	0.0462606	4.621912	3.8e-06	0.0021226	12.957643	12.929301	12.950216	12.685267	12.777846	12.817868
ENSGALG00010024658 CSF3	115.3912	0.8991719	0.1971158	4.561642	5.1e-06	0.0026918	9.824107	9.774722	9.875324	9.611872	9.604169	9.667846

dge1 <- dge
write.table(dge1,file="dge1.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge1,"Rfx: NTC gRNA vs DsRed gRNA")

make_volcano(dge1,"Rfx: NTC gRNA vs DsRed gRNA")

make_heatmap(dge1,"Rfx: NTC gRNA vs DsRed gRNA",ss1,xx1,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge1$gene <- sapply(strsplit(rownames(dge1)," "),"[[",2)
m1 <- aggregate(stat ~ gene,dge1,mean)
rownames(m1) <- m1$gene ; m1$gene <- NULL
mres1 <- mitch_calc(m1, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres1$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx: NTC gRNA vs DsRed gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx: NTC gRNA vs DsRed gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
779	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.6387612	0.0000000
937	GO:0031663 lipopolysaccharide-mediated signaling pathway	10	0.0013138	-0.5867493	0.0603783
615	GO:0009378 four-way junction helicase activity	13	0.0002707	0.5833426	0.0266915
435	GO:0006695 cholesterol biosynthetic process	16	0.0000586	0.5801984	0.0074325
1050	GO:0035145 exon-exon junction complex	10	0.0019511	0.5656887	0.0706760
1534	GO:0070182 DNA polymerase binding	13	0.0004891	0.5584911	0.0319696
1518	GO:0061749 forked DNA-dependent helicase activity	10	0.0024969	0.5522039	0.0777536
56	GO:0000793 condensed chromosome	13	0.0009537	0.5292204	0.0470477
1524	GO:0070034 telomerase RNA binding	12	0.0018880	0.5180580	0.0698152
135	GO:0001965 G-protein alpha-subunit binding	10	0.0046959	-0.5163223	0.1068619
166	GO:0003688 DNA replication origin binding	12	0.0020299	0.5144763	0.0720613
1724	GO:1902774 late endosome to lysosome transport	10	0.0059156	-0.5026722	0.1220952
805	GO:0030140 trans-Golgi network transport vesicle	11	0.0040195	-0.5009173	0.0967671
1086	GO:0036121 double-stranded DNA helicase activity	11	0.0040342	0.5007169	0.0967671
943	GO:0031929 TOR signaling	13	0.0018861	-0.4977983	0.0698152
781	GO:0022627 cytosolic small ribosomal subunit	25	0.0000172	-0.4966481	0.0038906
190	GO:0003899 DNA-directed 5’-3’ RNA polymerase activity	10	0.0072028	0.4907851	0.1404947
374	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0003277	-0.4760708	0.0276946
377	GO:0006268 DNA unwinding involved in DNA replication	14	0.0021458	0.4737925	0.0732445
1597	GO:0071806 protein transmembrane transport	12	0.0047896	-0.4703127	0.1076143

write.table(mres1$enrichment_result,file="mitch1.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres1$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 36

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx: NTC gRNA vs DsRed gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch1.html") ) {
  mitch_report(res=mres1,outfile="mitch1.html")
}

o <- dge1[order(dge1$baseMean),]
NBINS=20
brks <- quantile(o$baseMean, probs = c(seq(0,1,1/NBINS) ) )
o$quant <- cut(o$baseMean, breaks = brks, labels = 1:NBINS, include.lowest = TRUE)
z <- lapply(unique(o$quant) , function(i) { y <- subset(o,quant==i) ; y$log2FoldChange } )
names(z) <- 1:NBINS
XLAB=paste("Expression bin,",round(nrow(o)/NBINS),"genes each")
par(mar = c(5.1, 4.1, 4.1, 2.1) )
boxplot(z,ylim=c(-0.5,0.5),xlab=XLAB,ylab="Fold change",
  main="Rfx: NTC gRNA vs DsRed gRNA")
abline(h=0,col="red",lty=2,lwd=2)

There are some off target effects caused by the DsRed targeting gRNA.

2. Rfx_Tfx_ctl versus Rfx_NTC_gRNA

Now examining the effect of adding the NTC gRNA.

ss2 <- subset(ss,SampleGroup=="Rfx_Tfx_ctl" | SampleGroup=="Rfx_NTC_gRNA")
ss2$trt <- factor(grepl("NTC",ss2$SampleGroup))
rownames(ss2) <- ss2$SampleLabel
xx2 <- xx[,colnames(xx) %in% ss2$SampleLabel]
xx2f <- xx2[rowMeans(xx2)>=10,]
dim(xx2)

## [1] 30111     6

dim(xx2f)

## [1] 14612     6

dds <- DESeqDataSetFromMatrix(countData = xx2f , colData = ss2, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by NTC gRNA compared to transfection control") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by NTC gRNA compared to transfection control

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3	Rfx_Tfx_ctl_1	Rfx_Tfx_ctl_2	Rfx_Tfx_ctl_3
ENSGALG00010006303 ENSGALG00010006303	22.70974	-4.1955810	0.6651838	-6.307401	0.0000000	0.0000041	10.70987	10.79478	10.74731	10.92072	10.93430	10.96586
Fbxw11-GFP Fbxw11-GFP	19697.81377	0.5986278	0.1074391	5.571785	0.0000000	0.0001842	14.82634	14.62508	14.69901	14.30901	14.06570	14.22592
ENSGALG00010015565 ENSGALG00010015565	321.62945	0.5782863	0.1101357	5.250671	0.0000002	0.0007380	11.40690	11.38742	11.39188	11.27770	11.25566	11.28333
ENSGALG00010029439 ENSGALG00010029439	53.41059	-1.7167580	0.3318039	-5.174014	0.0000002	0.0008369	10.86412	10.93960	10.83401	11.01158	11.02694	11.04646
ENSGALG00010003609 ENSGALG00010003609	36.41332	1.6392185	0.3366935	4.868578	0.0000011	0.0032845	10.97456	10.96195	10.97688	10.85318	10.86208	10.85911
ENSGALG00010003545 ENSGALG00010003545	95.83669	1.0013060	0.2110708	4.743935	0.0000021	0.0046463	11.12821	11.11434	11.07764	11.00780	10.96274	11.00003
ENSGALG00010002373 POSTN	802.41475	0.4451845	0.0940849	4.731730	0.0000022	0.0046463	11.81707	11.69428	11.76187	11.61190	11.56352	11.66250
ENSGALG00010029479 ENSGALG00010029479	261.08343	0.5720726	0.1476259	3.875150	0.0001066	0.1670489	11.31782	11.26440	11.39484	11.19380	11.23448	11.22435
ENSGALG00010001609 ENSGALG00010001609	26.22936	-1.7170029	0.4438484	-3.868444	0.0001095	0.1670489	10.85621	10.84849	10.77472	10.92072	10.96005	10.91804
ENSGALG00010029468 ENSGALG00010029468	181.93398	1.2692646	0.3289992	3.857957	0.0001143	0.1670489	11.29650	11.34034	11.16198	11.10131	11.03328	11.08093
ENSGALG00010006565 ENSGALG00010006565	41.09502	1.3857372	0.3641313	3.805597	0.0001415	0.1878878	11.03473	10.96896	10.92799	10.84071	10.89806	10.89100
ENSGALG00010020995 ENSGALG00010020995	343.96907	0.4933009	0.1336160	3.691931	0.0002226	0.2679354	11.44676	11.41175	11.36764	11.32465	11.34235	11.23320
ENSGALG00010029501 ENSGALG00010029501	518.42604	0.4036486	0.1098543	3.674400	0.0002384	0.2679354	11.63096	11.51184	11.51375	11.47446	11.43979	11.42509
ENSGALG00010029462 ENSGALG00010029462	652.53692	-0.3228801	0.0932792	-3.461436	0.0005373	0.5607120	11.57527	11.53428	11.52933	11.57856	11.66660	11.68065
ENSGALG00010027753 ENSGALG00010027753	98.05600	0.7126180	0.2075726	3.433103	0.0005967	0.5812003	11.11393	11.11140	11.06614	11.04025	10.99310	11.00402
ENSGALG00010020993 ENSGALG00010020993	986.64075	-0.2416690	0.0719765	-3.357611	0.0007862	0.7178901	11.71375	11.74960	11.77821	11.79901	11.83880	11.86210
ENSGALG00010001531 ENSGALG00010001531	32.12185	-1.2296378	0.3697602	-3.325501	0.0008826	0.7585151	10.83425	10.86092	10.89318	10.93381	10.96805	10.93426
ENSGALG00010005536 TAF4B	155.89248	-0.5028528	0.1544522	-3.255718	0.0011311	0.9180434	11.11393	11.11140	11.10137	11.19148	11.17789	11.18360
ENSGALG00010019735 ENSGALG00010019735	426.47265	0.3593377	0.1116545	3.218299	0.0012895	0.9641510	11.51595	11.49308	11.40268	11.38217	11.39342	11.38111
ENSGALG00010008274 NFAT5	2116.12718	0.1561503	0.0486949	3.206710	0.0013426	0.9641510	12.31037	12.30360	12.26699	12.23259	12.21689	12.20334

dge2 <- dge
write.table(dge2,file="dge2.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge2,"Rfx: Tfx ctl vs NTC gRNA")

make_volcano(dge2,"Rfx: Tfx ctl vs NTC gRNA")

make_heatmap(dge2,"Rfx: Tfx ctl vs NTC gRNA",ss2,xx2,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge2$gene <- sapply(strsplit(rownames(dge2)," "),"[[",2)
m2 <- aggregate(stat ~ gene,dge2,mean)
rownames(m2) <- m2$gene ; m2$gene <- NULL
mres2 <- mitch_calc(m2, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres2$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx: Tfx ctl vs NTC gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx: Tfx ctl vs NTC gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
1691	GO:0140662 ATP-dependent protein folding chaperone	23	0.0000003	-0.6192405	0.0002419
1427	GO:0051315 attachment of mitotic spindle microtubules to kinetochore	10	0.0008472	-0.6093741	0.0328092
1625	GO:0090201 negative regulation of release of cytochrome c from mitochondria	11	0.0017417	-0.5452112	0.0514869
135	GO:0001965 G-protein alpha-subunit binding	10	0.0042201	0.5225364	0.0849583
36	GO:0000400 four-way junction DNA binding	14	0.0010045	-0.5077371	0.0356802
1110	GO:0042026 protein refolding	12	0.0023389	-0.5074125	0.0592357
958	GO:0032212 positive regulation of telomere maintenance via telomerase	17	0.0003508	-0.5007694	0.0228228
531	GO:0007339 binding of sperm to zona pellucida	11	0.0047070	-0.4921794	0.0905972
163	GO:0003678 DNA helicase activity	20	0.0001481	-0.4900838	0.0209802
118	GO:0001782 B cell homeostasis	12	0.0049993	0.4680165	0.0939853
473	GO:0007017 microtubule-based process	11	0.0077957	-0.4633425	0.1281966
1572	GO:0071233 cellular response to L-leucine	10	0.0116846	0.4605133	0.1704084
458	GO:0006913 nucleocytoplasmic transport	18	0.0007285	-0.4600079	0.0303438
154	GO:0003007 heart morphogenesis	16	0.0016006	0.4557102	0.0481822
833	GO:0030282 bone mineralization	24	0.0001259	0.4521741	0.0203305
1745	GO:1990166 protein localization to site of double-strand break	10	0.0146151	-0.4459511	0.1911528
1610	GO:0072542 protein phosphatase activator activity	10	0.0154031	-0.4424787	0.1977730
744	GO:0017116 single-stranded DNA helicase activity	16	0.0026153	-0.4346166	0.0615232
56	GO:0000793 condensed chromosome	13	0.0069061	-0.4327365	0.1179343
882	GO:0030863 cortical cytoskeleton	11	0.0129742	-0.4326464	0.1807799

write.table(mres2$enrichment_result,file="mitch2.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres2$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 59

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx: Tfx ctl vs NTC gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch2.html") ) {
  mitch_report(res=mres2,outfile="mitch2.html")
}

Could FBXW11 be a gene that is activated downstream of DsRed.

3. Rfx_NTfc_ctl versus Rfx_Tfx_ctl

ss3 <- subset(ss,SampleGroup=="Rfx_NTfx_ctl" | SampleGroup=="Rfx_Tfx_ctl")
ss3$trt <- factor(grepl("_Tfx",ss3$SampleGroup))
rownames(ss3) <- ss3$SampleLabel

xx3 <- xx[,colnames(xx) %in% ss3$SampleLabel]
xx3f <- xx3[rowMeans(xx3)>=10,]
dim(xx3)

## [1] 30111     6

dim(xx3f)

## [1] 14410     6

dds <- DESeqDataSetFromMatrix(countData = xx3f , colData = ss3, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by transfection in Rfx") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by transfection in Rfx

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Rfx_NTfx_ctl_1	Rfx_NTfx_ctl_2	Rfx_NTfx_ctl_3	Rfx_Tfx_ctl_1	Rfx_Tfx_ctl_2	Rfx_Tfx_ctl_3
Fbxw11-GFP Fbxw11-GFP	7663.1003	5.5328329	0.1064623	51.96987	0	0	11.07823	11.09253	11.07122	14.20839	13.95698	14.12272
ENSGALG00010006156 NFKBIZ	1377.9635	1.2274515	0.0586090	20.94304	0	0	11.48858	11.50677	11.45761	12.00869	12.02327	12.03853
ENSGALG00010022322 ENSGALG00010022322	2457.3842	0.9616168	0.0470333	20.44547	0	0	11.92477	11.89402	11.93224	12.42722	12.43647	12.47989
ENSGALG00010004483 STEAP4	1843.8349	0.9888966	0.0516806	19.13477	0	0	11.74661	11.68852	11.71688	12.21621	12.19929	12.21261
ENSGALG00010007787 C1orf198	19973.8427	-0.5029082	0.0302729	-16.61248	0	0	14.68887	14.65572	14.65494	14.19804	14.26734	14.19567
ENSGALG00010006816 ENSGALG00010006816	1017.1971	1.0143213	0.0655501	15.47398	0	0	11.38974	11.36635	11.38988	11.77225	11.75968	11.80539
ENSGALG00010025320 JAM3	6764.6010	-0.6325673	0.0439978	-14.37726	0	0	13.44786	13.36628	13.34556	12.91204	12.97579	12.86082
ENSGALG00010029178 NOTUM	1491.2883	-0.9612881	0.0693497	-13.86145	0	0	12.08913	12.03072	11.99314	11.58920	11.66044	11.54114
DsRed DsRed	7914.1464	16.3103829	1.1867882	13.74330	0	0	10.37505	10.37505	10.37505	14.28754	13.98998	14.20570
ENSGALG00010001230 PTX3	1296.6557	0.9377445	0.0687286	13.64416	0	0	11.48523	11.53166	11.55784	11.95844	11.86688	11.96865
ENSGALG00010005131 IL8L2	238.8007	2.1761892	0.1619716	13.43562	0	0	10.74986	10.68356	10.78549	11.13902	11.10426	11.19960
ENSGALG00010020187 MMP23B	2271.1355	0.7138761	0.0540994	13.19565	0	0	11.99605	11.90524	11.91430	12.35210	12.33732	12.28780
ENSGALG00010005086 ENSGALG00010005086	536.0812	1.1612824	0.0889803	13.05100	0	0	11.11110	11.07753	11.07571	11.45295	11.42726	11.40629
ENSGALG00010013207 C1QTNF3	4866.5343	-0.4797374	0.0399997	-11.99353	0	0	13.01473	12.99205	12.94659	12.61706	12.68564	12.66041
ENSGALG00010019277 IGFBP5	13152.1325	0.4175479	0.0353896	11.79859	0	0	13.76146	13.68912	13.80118	14.09885	14.08051	14.12961
ENSGALG00010018756 EML4	2064.6832	0.5676725	0.0490660	11.56957	0	0	11.88215	11.92095	11.94262	12.20283	12.21874	12.22004
ENSGALG00010005429 ZDHHC20	6047.6937	0.3952383	0.0348611	11.33751	0	0	12.89217	12.91233	12.89941	13.17498	13.16679	13.22402
ENSGALG00010015706 ENSGALG00010015706	1713.2009	0.6759317	0.0597946	11.30422	0	0	11.72360	11.76989	11.78342	12.05291	12.06468	12.14897
ENSGALG00010024266 ENSGALG00010024266	2123.2612	-0.6090502	0.0551053	-11.05249	0	0	12.25182	12.28288	12.20225	11.95629	11.93927	11.87221
ENSGALG00010025057 SCP2	12327.1468	-0.3086666	0.0282262	-10.93545	0	0	13.99725	13.97687	13.97376	13.74223	13.71670	13.71535

dge3 <- dge
write.table(dge3,file="dge3.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge3,"Rfx: NTfc ctl vs Tfx ctl")

make_volcano(dge3,"Rfx: NTfc ctl vs Tfx ctl")

make_heatmap(dge3,"Rfx: NTfc ctl vs Tfx ctl",ss3,xx3,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge3$gene <- sapply(strsplit(rownames(dge3)," "),"[[",2)
m3 <- aggregate(stat ~ gene,dge3,mean)
rownames(m3) <- m3$gene ; m3$gene <- NULL
mres3 <- mitch_calc(m3, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres3$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx: NTfc ctl vs Tfx ctl") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx: NTfc ctl vs Tfx ctl

	set	setSize	pANOVA	s.dist	p.adjustANOVA
1282	GO:0045943 positive regulation of transcription by RNA polymerase I	11	0.0000008	0.8611405	0.0000535
250	GO:0005049 nuclear export signal receptor activity	10	0.0000112	0.8021294	0.0005213
1686	GO:0140662 ATP-dependent protein folding chaperone	23	0.0000000	0.7504626	0.0000001
779	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.7359248	0.0000000
777	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.7211890	0.0000000
1047	GO:0035145 exon-exon junction complex	10	0.0000785	0.7210995	0.0025755
1520	GO:0070034 telomerase RNA binding	12	0.0000243	0.7037050	0.0008963
24	GO:0000228 nuclear chromosome	13	0.0000162	0.6904873	0.0006852
163	GO:0003678 DNA helicase activity	20	0.0000001	0.6889345	0.0000074
1083	GO:0036121 double-stranded DNA helicase activity	11	0.0000829	0.6852821	0.0026232
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.6802217	0.0000001
1746	GO:1990518 single-stranded 3’-5’ DNA helicase activity	11	0.0001110	0.6729702	0.0033925
26	GO:0000245 spliceosomal complex assembly	10	0.0005515	0.6308281	0.0120652
1197	GO:0043627 response to estrogen	10	0.0006530	0.6224635	0.0137127
148	GO:0002230 positive regulation of defense response to virus by host	11	0.0004561	0.6103480	0.0102730
1514	GO:0061749 forked DNA-dependent helicase activity	10	0.0009414	0.6040084	0.0179427
954	GO:0032212 positive regulation of telomere maintenance via telomerase	17	0.0000172	0.6020740	0.0006947
56	GO:0000793 condensed chromosome	13	0.0001952	0.5966719	0.0054055
741	GO:0017056 structural constituent of nuclear pore	19	0.0000102	0.5847965	0.0004875
515	GO:0007249 canonical NF-kappaB signal transduction	18	0.0000214	0.5785019	0.0008103

write.table(mres3$enrichment_result,file="mitch3.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres3$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 144

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx: NTfc ctl vs Tfx ctl",xlab="ES")
grid()

if ( ! file.exists("mitch3.html") ) {
  mitch_report(res=mres3,outfile="mitch3.html")
}

4. Hf_NTC_gRNA versus Hf_DsRed_gRNA

ss4 <- subset(ss,SampleGroup=="Hf_NTC_gRNA" | SampleGroup=="Hf_DsRed_gRNA")
ss4$trt <- factor(grepl("DsRed",ss4$SampleGroup))
rownames(ss4) <- ss4$SampleLabel

xx4 <- xx[,colnames(xx) %in% ss4$SampleLabel]
xx4f <- xx4[rowMeans(xx4)>=10,]
dim(xx4)

## [1] 30111     6

dim(xx4f)

## [1] 14412     6

dds <- DESeqDataSetFromMatrix(countData = xx4f , colData = ss4, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by DsRed targeting in Hf") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by DsRed targeting in Hf

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_DsRed_gRNA_1	Hf_DsRed_gRNA_2	Hf_DsRed_gRNA_3	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3
DsRed DsRed	7352.08543	-1.5278017	0.0545477	-28.008533	0	0.00e+00	12.304430	12.370374	12.292786	13.691844	13.520683	13.544964
Fbxw11-GFP Fbxw11-GFP	11668.72067	-1.4106970	0.0593248	-23.779219	0	0.00e+00	12.954770	12.918903	12.885270	14.303645	14.040457	14.138172
ENSGALG00010022322 ENSGALG00010022322	4495.62483	0.4402307	0.0439428	10.018268	0	0.00e+00	12.672143	12.679444	12.670138	12.352458	12.303991	12.335896
ENSGALG00010029945 CCL4	15140.89019	0.3971388	0.0551046	7.206992	0	0.00e+00	14.208537	14.089540	14.250945	13.847704	13.717597	13.892109
ENSGALG00010023206 HSPB9	260.75145	-0.8762077	0.1257649	-6.967030	0	0.00e+00	10.221969	10.234786	10.176926	10.458962	10.477272	10.441877
ENSGALG00010003545 ENSGALG00010003545	163.69670	1.2874066	0.1893633	6.798608	0	0.00e+00	10.229069	10.358763	10.301900	9.933674	10.073530	10.017911
ENSGALG00010012473 ENSGALG00010012473	320.82157	-0.9071409	0.1340674	-6.766306	0	0.00e+00	10.271866	10.285532	10.293674	10.668527	10.537800	10.476878
ENSGALG00010014984 LPP	1844.74856	-0.5102287	0.0774070	-6.591509	0	1.00e-07	11.517626	11.378005	11.412693	11.839862	11.772716	11.643400
ENSGALG00010013146 TNRC6B	1947.42852	-0.4049586	0.0627990	-6.448490	0	2.00e-07	11.559842	11.453182	11.545751	11.813786	11.786905	11.728374
ENSGALG00010022979 LAMA5	22283.56060	-0.3796638	0.0612877	-6.194777	0	7.00e-07	14.447199	14.258862	14.316458	14.818472	14.701789	14.568998
ENSGALG00010017057 NDUFAB1	2228.60631	0.3320152	0.0537368	6.178542	0	7.00e-07	11.853570	11.905614	11.898420	11.683380	11.640990	11.676247
ENSGALG00010025430 IRX5	1405.05185	0.5244179	0.0851349	6.159846	0	7.00e-07	11.590502	11.451509	11.485386	11.250507	11.094237	11.283739
ENSGALG00010024763 HOXD9	51.85914	1.9788799	0.3287829	6.018804	0	1.60e-06	9.995812	9.949732	9.983386	9.797615	9.650956	9.733606
ENSGALG00010006816 ENSGALG00010006816	1410.23510	0.3839373	0.0639937	5.999609	0	1.70e-06	11.522058	11.440866	11.470757	11.283642	11.235488	11.256630
ENSGALG00010024658 CSF3	137.23432	1.1157308	0.1863643	5.986827	0	1.70e-06	10.268536	10.188088	10.196254	10.006081	9.925697	10.027054
ENSGALG00010018888 ENSGALG00010018888	357.89695	-0.6430967	0.1083589	-5.934879	0	2.20e-06	10.403943	10.347576	10.372141	10.591928	10.607549	10.560079
ENSGALG00010015483 SGCD	1516.44549	-0.3459928	0.0601453	-5.752621	0	6.30e-06	11.359137	11.323602	11.302905	11.535146	11.545843	11.513301
ENSGALG00010018696 ENSGALG00010018696	879.22865	-0.4662691	0.0841336	-5.542007	0	2.02e-05	10.949169	10.842171	10.892515	11.170983	11.132963	11.051819
ENSGALG00010005913 NFKBIA	4135.89957	0.2560690	0.0465375	5.502424	0	2.40e-05	12.524904	12.514048	12.504439	12.347418	12.273541	12.331970
ENSGALG00010001704 ENSGALG00010001704	17388.94702	-0.3176630	0.0581646	-5.461451	0	2.87e-05	14.171830	13.927080	14.021872	14.398770	14.339768	14.272391

dge4 <- dge
write.table(dge4,file="dge4.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge4,"Hf: NTC gRNA vs DsRed gRNA")

make_volcano(dge4,"Hf: NTC gRNA vs DsRed gRNA")

make_heatmap(dge4,"Hf: NTC gRNA vs DsRed gRNA",ss4,xx4,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge4$gene <- sapply(strsplit(rownames(dge4)," "),"[[",2)
m4 <- aggregate(stat ~ gene,dge4,mean)
rownames(m4) <- m4$gene ; m4$gene <- NULL
mres4 <- mitch_calc(m4, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres4$enrichment_result,20) %>%
  kbl(caption = "Top GO in Hf: NTC gRNA vs DsRed gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Hf: NTC gRNA vs DsRed gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
145	GO:0002181 cytoplasmic translation	27	0.0000000	0.9176542	0.0000000
778	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.8789953	0.0000000
780	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.8238239	0.0000000
110	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0000010	0.8139527	0.0000658
347	GO:0005852 eukaryotic translation initiation factor 3 complex	10	0.0000159	0.7880601	0.0007383
180	GO:0003735 structural constituent of ribosome	118	0.0000000	0.7673033	0.0000000
409	GO:0006412 translation	96	0.0000000	0.6903554	0.0000000
700	GO:0016282 eukaryotic 43S preinitiation complex	14	0.0000127	0.6737194	0.0006060
345	GO:0005840 ribosome	43	0.0000000	0.6682555	0.0000000
1002	GO:0033290 eukaryotic 48S preinitiation complex	13	0.0000316	0.6664936	0.0012991
315	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	0.6650277	0.0000000
438	GO:0006695 cholesterol biosynthetic process	16	0.0000060	0.6536179	0.0002932
184	GO:0003774 cytoskeletal motor activity	12	0.0001145	-0.6430759	0.0037858
266	GO:0005201 extracellular matrix structural constituent	12	0.0001158	-0.6426235	0.0037858
994	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	0.6261081	0.0000000
1680	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0003707	-0.6198911	0.0093467
758	GO:0019843 rRNA binding	22	0.0000006	0.6166561	0.0000374
1650	GO:0098761 cellular response to interleukin-7	10	0.0007681	0.6143334	0.0152325
597	GO:0008320 protein transmembrane transporter activity	11	0.0004901	0.6070116	0.0112339
1232	GO:0045271 respiratory chain complex I	31	0.0000000	0.6032018	0.0000006

write.table(mres4$enrichment_result,file="mitch4.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres4$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 134

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Hf: NTC gRNA vs DsRed gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch4.html") ) {
  mitch_report(res=mres4,outfile="mitch4.html")
}

o <- dge4[order(dge4$baseMean),]
NBINS=20
brks <- quantile(o$baseMean, probs = c(seq(0,1,1/NBINS) ) )
o$quant <- cut(o$baseMean, breaks = brks, labels = 1:NBINS, include.lowest = TRUE)
z <- lapply(unique(o$quant) , function(i) { y <- subset(o,quant==i) ; y$log2FoldChange } )
names(z) <- 1:NBINS
XLAB=paste("Expression bin,",round(nrow(o)/NBINS),"genes each")
par(mar = c(5.1, 4.1, 4.1, 2.1) )
boxplot(z,ylim=c(-0.5,0.5),xlab=XLAB,ylab="Fold change",
  main="Hf: NTC gRNA vs DsRed gRNA")
abline(h=0,col="red",lty=2,lwd=2)

5. Hf_Tfx_ctl versus Hf_NTC_gRNA

Now examining the effect of adding the NTC gRNA.

ss5 <- subset(ss,SampleGroup=="Hf_Tfx_ctl" | SampleGroup=="Hf_NTC_gRNA")
ss5$trt <- factor(grepl("NTC",ss5$SampleGroup))
rownames(ss5) <- ss5$SampleLabel
xx5 <- xx[,colnames(xx) %in% ss5$SampleLabel]
xx5f <- xx5[rowMeans(xx5)>=10,]
dim(xx5)

## [1] 30111     6

dim(xx5f)

## [1] 14357     6

dds <- DESeqDataSetFromMatrix(countData = xx5f , colData = ss5, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by NTC gRNA compared to transfection control") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by NTC gRNA compared to transfection control

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3	Hf_Tfx_ctl_1	Hf_Tfx_ctl_2	Hf_Tfx_ctl_3
ENSGALG00010006303 ENSGALG00010006303	17.99008	-3.7098356	0.6683447	-5.550782	0.0000000	0.0004082	10.47085	10.41023	10.37142	10.60294	10.58012	10.61780
Fbxw11-GFP Fbxw11-GFP	13631.83580	0.5497622	0.1094397	5.023427	0.0000005	0.0036424	14.32532	14.07783	14.17050	13.68715	13.65521	13.84013
ENSGALG00010003545 ENSGALG00010003545	126.52418	-0.8367108	0.2120322	-3.946150	0.0000794	0.3799351	10.69192	10.79330	10.75310	10.90411	10.90173	10.80841
ENSGALG00010009114 SETD7	5938.92628	0.1358263	0.0358463	3.789127	0.0001512	0.5317574	13.08626	13.07643	13.08018	12.97351	12.95992	13.01302
ENSGALG00010019257 SLC1A4	5563.70273	0.1558961	0.0417022	3.738318	0.0001853	0.5317574	13.07137	12.99554	12.98867	12.88727	12.89179	12.94400
ENSGALG00010010914 ENSGALG00010010914	20.67785	1.6274111	0.4571709	3.559743	0.0003712	0.8879518	10.62045	10.58024	10.57466	10.50708	10.50469	10.47869
ENSGALG00010016508 CREBBP	2655.82671	0.1798783	0.0524290	3.430897	0.0006016	0.9998869	12.35346	12.33876	12.27839	12.17325	12.22361	12.26377
ENSGALG00010020691 ENSGALG00010020691	831.78306	-0.2484451	0.0738304	-3.365079	0.0007652	0.9998869	11.47163	11.42157	11.42498	11.55611	11.52101	11.51650
ENSGALG00010010496 MRPL17	1724.33797	-0.1797901	0.0538396	-3.339363	0.0008397	0.9998869	11.92661	11.87809	11.87526	11.99322	11.98201	11.97362
ENSGALG00010019068 SNRPA1	3875.88248	-0.1401798	0.0422012	-3.321706	0.0008947	0.9998869	12.56024	12.58803	12.53698	12.68566	12.62671	12.64628
ENSGALG00010019231 ENSGALG00010019231	790.13465	0.2950281	0.0889440	3.317008	0.0009099	0.9998869	11.46636	11.58294	11.47954	11.39378	11.37366	11.44658
ENSGALG00010005913 NFKBIA	3786.02649	-0.1415096	0.0437258	-3.236293	0.0012109	0.9998869	12.56474	12.50213	12.55230	12.66476	12.62376	12.60428
ENSGALG00010009392 MAPRE3	1810.30339	-0.1750275	0.0542114	-3.228610	0.0012439	0.9998869	11.91448	11.92532	11.94839	12.04316	11.97790	12.02994
ENSGALG00010017157 ENSGALG00010017157	12.90598	-2.2013829	0.6833555	-3.221431	0.0012755	0.9998869	10.47085	10.37142	10.48124	10.58220	10.57225	10.49304
ENSGALG00010015892 UBE2C	5334.69523	-0.1159327	0.0360174	-3.218795	0.0012873	0.9998869	12.88370	12.87372	12.88197	12.96990	12.97743	12.93787
ENSGALG00010021087 THBS1	68518.08634	0.0866757	0.0270778	3.200995	0.0013695	0.9998869	16.14315	16.15807	16.17732	16.04513	16.08494	16.09761
ENSGALG00010028900 TNC	20571.05631	0.1380925	0.0433031	3.188976	0.0014278	0.9998869	14.59575	14.55086	14.52822	14.32595	14.47739	14.49736
ENSGALG00010017057 NDUFAB1	2014.91869	-0.1872812	0.0587989	-3.185113	0.0014470	0.9998869	12.01684	11.98261	12.01161	12.17225	12.06066	12.06902
ENSGALG00010026886 JUN	3840.33281	-0.1291057	0.0411985	-3.133748	0.0017259	0.9998869	12.57349	12.55922	12.53838	12.67289	12.62229	12.62695
ENSGALG00010024763 HOXD9	34.62589	-1.3015227	0.4205421	-3.094869	0.0019690	0.9998869	10.59358	10.48782	10.54746	10.71438	10.59858	10.62110

dge5 <- dge
write.table(dge5,file="dge2.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge5,"Hf: Tfx ctl vs NTC gRNA")

make_volcano(dge5,"Hf: Tfx ctl vs NTC gRNA")

make_heatmap(dge5,"Hf: Tfx ctl vs NTC gRNA",ss5,xx5,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge5$gene <- sapply(strsplit(rownames(dge5)," "),"[[",2)
m5 <- aggregate(stat ~ gene,dge5,mean)
rownames(m5) <- m5$gene ; m5$gene <- NULL
mres5 <- mitch_calc(m5, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres5$enrichment_result,20) %>%
  kbl(caption = "Top GO in Hf: Tfx ctl vs NTC gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Hf: Tfx ctl vs NTC gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
1122	GO:0042474 middle ear morphogenesis	10	0.0008817	-0.6073479	0.0851343
26	GO:0000245 spliceosomal complex assembly	10	0.0021049	-0.5615702	0.1234197
775	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.5543046	0.0000008
777	GO:0022627 cytosolic small ribosomal subunit	25	0.0000101	-0.5102279	0.0025269
781	GO:0030016 myofibril	13	0.0023614	0.4870648	0.1298008
1731	GO:1990391 DNA repair complex	10	0.0082282	0.4826151	0.2378757
1727	GO:1990023 mitotic spindle midzone	10	0.0083479	-0.4817210	0.2378757
683	GO:0015813 L-glutamate transmembrane transport	11	0.0060594	0.4779268	0.1937918
99	GO:0001671 ATPase activator activity	18	0.0004729	-0.4759398	0.0669946
1646	GO:0098761 cellular response to interleukin-7	10	0.0095869	-0.4730879	0.2430366
144	GO:0002181 cytoplasmic translation	27	0.0000215	-0.4724372	0.0047235
436	GO:0006730 one-carbon metabolic process	13	0.0036063	0.4662676	0.1475211
805	GO:0030155 regulation of cell adhesion	19	0.0004951	0.4616335	0.0669946
158	GO:0003333 amino acid transmembrane transport	10	0.0124130	0.4566180	0.2726603
1031	GO:0035025 positive regulation of Rho protein signal transduction	10	0.0126236	0.4555284	0.2726603
134	GO:0001965 G-protein alpha-subunit binding	10	0.0131981	0.4526367	0.2726603
313	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000015	-0.4516761	0.0006038
1532	GO:0070412 R-SMAD binding	13	0.0048403	-0.4513387	0.1742071
1237	GO:0045505 dynein intermediate chain binding	18	0.0009780	0.4488706	0.0851343
1414	GO:0051382 kinetochore assembly	11	0.0106381	-0.4448170	0.2563354

write.table(mres5$enrichment_result,file="mitch5.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres5$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 9

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Hf: Tfx ctl vs NTC gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch5.html") ) {
  mitch_report(res=mres5,outfile="mitch5.html")
}

6. Hf_NTfx_ctl versus Hf_Tfx_ctl

ss6 <- subset(ss,SampleGroup=="Hf_NTfx_ctl" | SampleGroup=="Hf_Tfx_ctl")
ss6$trt <- factor(grepl("_Tfx",ss6$SampleGroup))
rownames(ss6) <- ss6$SampleLabel

xx6 <- xx[,colnames(xx) %in% ss6$SampleLabel]
xx6f <- xx6[rowMeans(xx6)>=10,]
dim(xx6)

## [1] 30111     6

dim(xx6f)

## [1] 14256     6

dds <- DESeqDataSetFromMatrix(countData = xx6f , colData = ss6, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by transfection in Hf") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by transfection in Hf

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_NTfx_ctl_1	Hf_NTfx_ctl_2	Hf_NTfx_ctl_3	Hf_Tfx_ctl_1	Hf_Tfx_ctl_2	Hf_Tfx_ctl_3
Fbxw11-GFP Fbxw11-GFP	5562.3847	5.1247670	0.1016093	50.43599	0	0	10.491681	10.414559	10.482699	13.53642	13.49822	13.70475
ENSGALG00010006156 NFKBIZ	1367.0066	1.4866924	0.0612919	24.25593	0	0	10.904194	10.900830	10.931792	11.72612	11.67492	11.72378
ENSGALG00010022322 ENSGALG00010022322	2739.2450	1.1082722	0.0478461	23.16328	0	0	11.558490	11.584260	11.582534	12.30698	12.30108	12.36787
ENSGALG00010004483 STEAP4	1517.0361	1.1008179	0.0605650	18.17580	0	0	11.059705	11.111586	11.156249	11.71324	11.73911	11.74660
ENSGALG00010025320 JAM3	5686.3159	-0.6266981	0.0370779	-16.90218	0	0	13.031027	13.028330	13.015446	12.53639	12.52881	12.49503
ENSGALG00010020187 MMP23B	1771.8426	0.8953471	0.0533494	16.78270	0	0	11.286628	11.290545	11.312623	11.85778	11.83214	11.81102
ENSGALG00010001230 PTX3	1050.4939	1.0418355	0.0689284	15.11476	0	0	10.870595	10.857983	10.921752	11.44021	11.37314	11.39107
ENSGALG00010007787 C1orf198	19044.0873	-0.4707627	0.0318848	-14.76450	0	0	14.513356	14.515700	14.550119	14.07114	14.07581	14.12087
ENSGALG00010006816 ENSGALG00010006816	914.0890	1.1176593	0.0797470	14.01506	0	0	10.820634	10.706765	10.806531	11.32138	11.25897	11.32980
ENSGALG00010019277 IGFBP5	11189.2227	0.4678286	0.0345033	13.55895	0	0	13.396896	13.379511	13.434935	13.80395	13.80327	13.84813
DsRed DsRed	5403.0502	15.8377798	1.1845114	13.37073	0	0	9.550636	9.550636	9.550636	13.48935	13.51133	13.73163
ENSGALG00010021087 THBS1	56029.8205	0.4628356	0.0358267	12.91875	0	0	15.585215	15.500877	15.618115	15.98807	16.02434	16.04573
ENSGALG00010017751 ST3GAL5	3061.3187	-0.5440060	0.0447750	-12.14978	0	0	12.330073	12.276981	12.272544	11.91086	11.91790	11.89650
ENSGALG00010029178 NOTUM	2075.3883	-0.6040073	0.0498738	-12.11072	0	0	11.892080	11.939610	11.908106	11.53241	11.53767	11.51985
ENSGALG00010005495 TXN	17201.8015	-0.3827643	0.0324677	-11.78909	0	0	14.391651	14.332222	14.338159	14.01775	13.98735	13.99845
ENSGALG00010013207 C1QTNF3	5734.9218	-0.4273741	0.0364043	-11.73967	0	0	12.956719	12.963672	12.981135	12.61916	12.63506	12.61060
ENSGALG00010018756 EML4	1934.5497	0.6277787	0.0548149	11.45270	0	0	11.483140	11.417190	11.484963	11.84222	11.82384	11.88839
ENSGALG00010005131 IL8L2	158.2816	2.0421669	0.1813879	11.25856	0	NA	9.973319	9.907713	10.006519	10.35711	10.35211	10.42679
ENSGALG00010016647 EGR1	2991.8776	-0.5471527	0.0500784	-10.92593	0	0	12.209319	12.306032	12.292123	11.92380	11.87731	11.85340
ENSGALG00010005086 ENSGALG00010005086	363.2340	1.1336860	0.1061348	10.68157	0	0	10.353130	10.333592	10.319582	10.72597	10.69374	10.67046

dge6 <- dge
write.table(dge6,file="dge6.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge6,"Hf: NTfc ctl vs Tfx ctl")

make_volcano(dge6,"Hf: NTfc ctl vs Tfx ctl")

make_heatmap(dge6,"Hf: NTfc ctl vs Tfx ctl",ss6,xx6,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge6$gene <- sapply(strsplit(rownames(dge6)," "),"[[",2)
m6 <- aggregate(stat ~ gene,dge6,mean)
rownames(m6) <- m6$gene ; m6$gene <- NULL
mres6 <- mitch_calc(m6, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres6$enrichment_result,20) %>%
  kbl(caption = "Top GO in Hf: NTfc ctl vs Tfx ctl") %>%
  kable_paper("hover", full_width = F)

Top GO in Hf: NTfc ctl vs Tfx ctl

	set	setSize	pANOVA	s.dist	p.adjustANOVA
775	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.8945404	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.8378811	0.0000000
1215	GO:0045116 protein neddylation	10	0.0000125	-0.7977631	0.0004218
777	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.7877051	0.0000000
483	GO:0007076 mitotic chromosome condensation	10	0.0000224	0.7741277	0.0006562
1076	GO:0036121 double-stranded DNA helicase activity	11	0.0000237	0.7359895	0.0006707
1041	GO:0035145 exon-exon junction complex	10	0.0000588	0.7336522	0.0013075
179	GO:0003735 structural constituent of ribosome	118	0.0000000	-0.7117635	0.0000000
771	GO:0021915 neural tube development	15	0.0000032	0.6942087	0.0001455
1708	GO:1902774 late endosome to lysosome transport	10	0.0001443	-0.6940208	0.0027582
24	GO:0000228 nuclear chromosome	13	0.0000326	0.6654525	0.0008178
1510	GO:0070034 telomerase RNA binding	12	0.0000835	0.6558440	0.0017686
249	GO:0005049 nuclear export signal receptor activity	10	0.0003750	0.6495779	0.0054931
1733	GO:1990518 single-stranded 3’-5’ DNA helicase activity	11	0.0001981	0.6479775	0.0034473
56	GO:0000793 condensed chromosome	13	0.0000550	0.6459904	0.0012564
1749	GO:2000737 negative regulation of stem cell differentiation	12	0.0001122	0.6439309	0.0022409
1675	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0002321	0.6409682	0.0039239
1504	GO:0061749 forked DNA-dependent helicase activity	10	0.0004562	0.6401097	0.0062655
1726	GO:1990023 mitotic spindle midzone	10	0.0006895	0.6197524	0.0085363
740	GO:0017056 structural constituent of nuclear pore	19	0.0000052	0.6036202	0.0002137

write.table(mres6$enrichment_result,file="mitch6.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres6$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 295

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Hf: NTfc ctl vs Tfx ctl",xlab="ES")
grid()

if ( ! file.exists("mitch6.html") ) {
  mitch_report(res=mres6,outfile="mitch6.html")
}

7. Rfx_DsRed_gRNA versus Hf_DsRed_gRNA

ss7 <- subset(ss,SampleGroup=="Rfx_DsRed_gRNA" | SampleGroup=="Hf_DsRed_gRNA")
ss7$trt <- factor(grepl("Hf",ss7$SampleGroup))
rownames(ss7) <- ss7$SampleLabel

xx7 <- xx[,colnames(xx) %in% ss7$SampleLabel]
xx7f <- xx7[rowMeans(xx7)>=10,]
dim(xx7)

## [1] 30111     6

dim(xx7f)

## [1] 14466     6

dds <- DESeqDataSetFromMatrix(countData = xx7f , colData = ss7, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes in Hf as compared to Rfx (DsRed gDNA)") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes in Hf as compared to Rfx (DsRed gDNA)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_DsRed_gRNA_1	Hf_DsRed_gRNA_2	Hf_DsRed_gRNA_3	Rfx_DsRed_gRNA_1	Rfx_DsRed_gRNA_2	Rfx_DsRed_gRNA_3
ENSGALG00010005474 COL8A1	5181.22579	1.4072932	0.0528465	26.62982	0	0	13.114134	13.042797	13.023945	11.907948	11.844188	11.938588
ENSGALG00010009793 TNFAIP6	4990.91092	1.1189490	0.0597584	18.72456	0	0	12.897318	12.927321	12.989725	11.909130	12.048842	12.073183
ENSGALG00010024266 ENSGALG00010024266	1306.26597	-1.3025659	0.0739627	-17.61113	0	0	10.691750	10.676942	10.657510	11.516537	11.479617	11.371810
ENSGALG00010005417 LOXL3	724.79066	1.5672084	0.0936381	16.73687	0	0	11.036075	10.942040	10.909233	10.187861	10.173152	10.260810
ENSGALG00010002949 ENSGALG00010002949	549403.78936	-0.7207412	0.0462421	-15.58625	0	0	18.734264	18.624489	18.637435	19.443675	19.370993	19.339998
ENSGALG00010016647 EGR1	3974.38230	-0.9745205	0.0626304	-15.55987	0	0	11.883861	11.885484	11.733738	12.551294	12.665595	12.626883
ENSGALG00010023078 TAGLN	4192.37168	-0.7629777	0.0500159	-15.25470	0	0	12.036081	11.974738	11.984511	12.648647	12.591584	12.610870
ENSGALG00010009067 ENSGALG00010009067	215523.51739	-0.7116253	0.0469641	-15.15255	0	0	17.409784	17.279862	17.289064	18.076927	18.032847	17.995610
ENSGALG00010007873 ENSGALG00010007873	460982.64947	-0.7402256	0.0489905	-15.10958	0	0	18.470030	18.353147	18.379165	19.214851	19.128931	19.073059
ENSGALG00010018616 ENSGALG00010018616	554903.05730	-0.7118731	0.0477575	-14.90601	0	0	18.740584	18.641571	18.674272	19.476403	19.371140	19.338455
ENSGALG00010001707 ENSGALG00010001707	180.88805	-6.5241663	0.4402708	-14.81853	0	0	9.174015	9.164188	9.144721	10.264639	10.109445	10.209583
ENSGALG00010013345 DCN	5137.39991	0.8515898	0.0581524	14.64411	0	0	12.862560	12.896837	12.929911	12.079362	12.234455	12.234527
ENSGALG00010024254 ENSGALG00010024254	316.79494	-1.9081965	0.1317474	-14.48376	0	0	9.769946	9.763161	9.736435	10.481080	10.342392	10.354795
ENSGALG00010018583 PPL	185.06075	2.5111425	0.1801697	13.93765	0	0	10.128891	10.154916	10.098050	9.549723	9.528337	9.426412
ENSGALG00010027337 AQP1	22210.83490	0.5544383	0.0403725	13.73307	0	0	14.718231	14.736076	14.781359	14.212548	14.211652	14.223964
ENSGALG00010007529 CDH13	17412.86597	0.5791960	0.0435489	13.29989	0	0	14.466430	14.409019	14.379171	13.842298	13.886576	13.887645
ENSGALG00010001347 MYOM1	520.80020	1.2918543	0.0996893	12.95881	0	0	10.625047	10.629195	10.727835	10.107849	10.124398	10.088422
ENSGALG00010022843 ENSGALG00010022843	87.34403	-4.8555794	0.3879446	-12.51617	0	0	9.220588	9.147228	9.189085	9.884512	9.844146	9.803835
ENSGALG00010023113 MYZAP	89.62103	4.0473234	0.3234829	12.51171	0	0	9.892770	9.773686	9.863557	9.249155	9.236873	9.230026
ENSGALG00010005836 NEBL	93.25190	3.7843647	0.3033558	12.47500	0	0	9.793911	9.871053	9.900015	9.193890	9.294744	9.280911

dge7 <- dge
write.table(dge7,file="dge7.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge7,"Rfx DsRed gRNA vs Hf DsRed gRNA")

make_volcano(dge7,"Rfx DsRed gRNA vs Hf DsRed gRNA")

make_heatmap(dge7,"Rfx DsRed gRNA vs Hf DsRed gRNA",ss7,xx7,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge7$gene <- sapply(strsplit(rownames(dge7)," "),"[[",2)
m7 <- aggregate(stat ~ gene,dge7,mean)
rownames(m7) <- m7$gene ; m7$gene <- NULL
mres7 <- mitch_calc(m7, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres7$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx DsRed gRNA vs Hf DsRed gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx DsRed gRNA vs Hf DsRed gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
778	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.8398847	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	0.8156863	0.0000000
780	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.8008800	0.0000000
1047	GO:0035145 exon-exon junction complex	10	0.0000317	-0.7598170	0.0017526
1686	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0000348	-0.7206734	0.0018312
180	GO:0003735 structural constituent of ribosome	118	0.0000000	0.6850152	0.0000000
572	GO:0008137 NADH dehydrogenase (ubiquinone) activity	10	0.0002977	0.6605657	0.0105434
1237	GO:0045271 respiratory chain complex I	31	0.0000000	0.6574504	0.0000000
377	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000008	0.6540957	0.0000741
440	GO:0006749 glutathione metabolic process	15	0.0000138	0.6482604	0.0008411
209	GO:0004602 glutathione peroxidase activity	12	0.0001049	0.6466394	0.0045325
200	GO:0004364 glutathione transferase activity	10	0.0005088	0.6347892	0.0166852
409	GO:0006412 translation	96	0.0000000	0.5956258	0.0000000
26	GO:0000245 spliceosomal complex assembly	10	0.0014952	-0.5799362	0.0351385
1719	GO:1902774 late endosome to lysosome transport	10	0.0015372	0.5784664	0.0351385
1742	GO:1990391 DNA repair complex	10	0.0017168	-0.5725735	0.0380066
995	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	0.5611098	0.0000004
482	GO:0007042 lysosomal lumen acidification	10	0.0023029	0.5566556	0.0463468
1225	GO:0045116 protein neddylation	10	0.0024540	0.5531614	0.0475309
1660	GO:0098869 cellular oxidant detoxification	28	0.0000006	0.5461365	0.0000557

write.table(mres7$enrichment_result,file="mitch7.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres7$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 92

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx DsRed gRNA vs Hf DsRed gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch7.html") ) {
  mitch_report(res=mres7,outfile="mitch7.html")
}

8. Rfx_NTC_gRNA versus Hf_NTC_gRNA

ss8 <- subset(ss,SampleGroup=="Rfx_NTC_gRNA" | SampleGroup=="Hf_NTC_gRNA")
ss8$trt <- factor(grepl("Hf",ss8$SampleGroup))
rownames(ss8) <- ss8$SampleLabel

xx8 <- xx[,colnames(xx) %in% ss8$SampleLabel]
xx8f <- xx8[rowMeans(xx8)>=10,]
dim(xx8)

## [1] 30111     6

dim(xx8f)

## [1] 14531     6

dds <- DESeqDataSetFromMatrix(countData = xx8f , colData = ss8, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes in Hf as compared to Rfx (NTC gDNA)") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes in Hf as compared to Rfx (NTC gDNA)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3
ENSGALG00010005474 COL8A1	4887.28026	1.4000111	0.0407401	34.36442	0	0	13.033606	13.030229	13.074111	11.936452	11.955421	11.993068
ENSGALG00010009793 TNFAIP6	4332.14842	0.9624357	0.0473261	20.33624	0	0	12.810498	12.752565	12.806856	12.069158	11.983858	12.103671
ENSGALG00010005417 LOXL3	733.65009	1.6386027	0.0851521	19.24325	0	0	11.158177	11.184751	11.198123	10.520789	10.433491	10.486890
ENSGALG00010016647 EGR1	3703.09860	-0.9145679	0.0514414	-17.77882	0	0	11.976831	11.893629	11.875465	12.653501	12.535289	12.594580
ENSGALG00010027337 AQP1	20305.86910	0.5719634	0.0322798	17.71892	0	0	14.617402	14.668769	14.659928	14.085251	14.105541	14.151400
ENSGALG00010007529 CDH13	17521.62235	0.5279203	0.0317741	16.61480	0	0	14.439400	14.417375	14.438277	13.903476	13.972360	13.952394
ENSGALG00010024266 ENSGALG00010024266	1335.25142	-1.1921850	0.0729832	-16.33506	0	0	10.907831	10.994855	10.964448	11.623680	11.674189	11.524813
ENSGALG00010001707 ENSGALG00010001707	177.23524	-5.9380826	0.3765330	-15.77042	0	0	9.585784	9.575079	9.663097	10.417210	10.499499	10.491322
ENSGALG00010013345 DCN	5257.02603	0.6339014	0.0411284	15.41273	0	0	12.939730	12.895469	12.928771	12.437493	12.362882	12.440636
ENSGALG00010001966 ARHGAP6	374.36645	1.5853355	0.1097062	14.45074	0	0	10.717642	10.710802	10.737895	10.194791	10.216143	10.227198
ENSGALG00010018583 PPL	193.48107	2.2727200	0.1611930	14.09937	0	0	10.424759	10.437238	10.431573	9.941605	9.897293	9.920485
ENSGALG00010023078 TAGLN	3886.16447	-0.6328938	0.0461591	-13.71115	0	0	12.034835	12.147572	12.086509	12.551823	12.583169	12.568158
ENSGALG00010022071 ACTA2	44250.86967	-0.3876839	0.0282967	-13.70065	0	0	15.269904	15.286664	15.280297	15.622160	15.672283	15.668759
ENSGALG00010006546 SMOC2	8454.33630	0.5122392	0.0387624	13.21485	0	0	13.465919	13.465196	13.470103	13.021671	12.968935	13.084408
ENSGALG00010024570 TPM1	45795.87484	-0.3664513	0.0285959	-12.81483	0	0	15.336329	15.369464	15.308541	15.685148	15.702617	15.693540
ENSGALG00010004981 ENSGALG00010004981	12008.18379	0.3869920	0.0316129	12.24160	0	0	13.868992	13.882645	13.868176	13.520056	13.551422	13.505426
ENSGALG00010022843 ENSGALG00010022843	84.40964	-5.1420390	0.4299877	-11.95857	0	0	9.628519	9.575079	9.581480	10.158912	10.181125	10.162026
ENSGALG00010025025 ENSGALG00010025025	130.72224	-2.4785877	0.2073377	-11.95435	0	0	9.816418	9.802509	9.847407	10.234036	10.327206	10.261521
ENSGALG00010009749 NDRG1	3131.46592	0.5252286	0.0442918	11.85837	0	0	12.295811	12.288884	12.298751	11.917419	11.885895	11.946596
ENSGALG00010005533 OPN4-1	79.28905	4.7699728	0.4142018	11.51606	0	0	10.232911	10.109868	10.077548	9.610357	9.599753	9.621508

dge8 <- dge
write.table(dge8,file="dge8.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge8,"Rfx NTC gRNA vs Hf NTC gRNA")

make_volcano(dge8,"Rfx NTC gRNA vs Hf NTC gRNA")

make_heatmap(dge8,"Rfx NTC gRNA vs Hf NTC gRNA",ss8,xx8,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge8$gene <- sapply(strsplit(rownames(dge8)," "),"[[",2)
m8 <- aggregate(stat ~ gene,dge8,mean)
rownames(m8) <- m8$gene ; m8$gene <- NULL
mres8 <- mitch_calc(m8, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres8$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx NTC gRNA vs Hf NTC gRNA") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx NTC gRNA vs Hf NTC gRNA

	set	setSize	pANOVA	s.dist	p.adjustANOVA
26	GO:0000245 spliceosomal complex assembly	10	0.0000652	-0.7291836	0.0086225
779	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.6493016	0.0000000
1047	GO:0035145 exon-exon junction complex	10	0.0004119	-0.6450487	0.0331612
781	GO:0022627 cytosolic small ribosomal subunit	25	0.0000005	0.5808455	0.0002931
799	GO:0030127 COPII vesicle coat	10	0.0020106	0.5640605	0.0890191
1737	GO:1905564 positive regulation of vascular endothelial cell proliferation	10	0.0021784	0.5596991	0.0899115
1771	GO:2001244 positive regulation of intrinsic apoptotic signaling pathway	11	0.0016124	0.5491436	0.0793190
1723	GO:1903076 regulation of protein localization to plasma membrane	18	0.0000663	0.5430965	0.0086225
145	GO:0002181 cytoplasmic translation	27	0.0000017	0.5323057	0.0006312
1724	GO:1903077 negative regulation of protein localization to plasma membrane	13	0.0010190	0.5262392	0.0622315
200	GO:0004364 glutathione transferase activity	10	0.0057273	0.5046029	0.1410611
209	GO:0004602 glutathione peroxidase activity	12	0.0025299	0.5034624	0.0974022
1577	GO:0071353 cellular response to interleukin-4	10	0.0059226	0.5026016	0.1410611
217	GO:0004683 calcium/calmodulin-dependent protein kinase activity	11	0.0055974	0.4824393	0.1410611
1213	GO:0044322 endoplasmic reticulum quality control compartment	11	0.0058386	0.4800427	0.1410611
1443	GO:0051642 centrosome localization	14	0.0019056	0.4792375	0.0865340
474	GO:0007017 microtubule-based process	11	0.0059738	0.4787377	0.1410611
788	GO:0030020 extracellular matrix structural constituent conferring tensile strength	16	0.0009995	0.4751812	0.0622315
798	GO:0030100 regulation of endocytosis	12	0.0047672	0.4705639	0.1321895
609	GO:0008593 regulation of Notch signaling pathway	10	0.0117287	0.4602719	0.1790655

write.table(mres8$enrichment_result,file="mitch8.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres8$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 26

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx NTC gRNA vs Hf NTC gRNA",xlab="ES")
grid()

if ( ! file.exists("mitch8.html") ) {
  mitch_report(res=mres8,outfile="mitch8.html")
}

9. Rfx_Tfx_ctl versus Hf_Tfx_ctl

ss9 <- subset(ss,SampleGroup=="Rfx_Tfx_ctl" | SampleGroup=="Hf_Tfx_ctl")
ss9$trt <- factor(grepl("Hf",ss9$SampleGroup))
rownames(ss9) <- ss9$SampleLabel

xx9 <- xx[,colnames(xx) %in% ss9$SampleLabel]
xx9f <- xx9[rowMeans(xx9)>=10,]
dim(xx9)

## [1] 30111     6

dim(xx9f)

## [1] 14432     6

dds <- DESeqDataSetFromMatrix(countData = xx9f , colData = ss9, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes in Hf as compared to Rfx (Tfx ctl)") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes in Hf as compared to Rfx (Tfx ctl)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_Tfx_ctl_1	Hf_Tfx_ctl_2	Hf_Tfx_ctl_3	Rfx_Tfx_ctl_1	Rfx_Tfx_ctl_2	Rfx_Tfx_ctl_3
ENSGALG00010005474 COL8A1	5058.9658	1.3742414	0.0342368	40.13926	0	0	13.39912	13.38567	13.37769	12.54041	12.54280	12.54890
ENSGALG00010007873 ENSGALG00010007873	448243.5748	-0.6993259	0.0274983	-25.43159	0	0	18.38479	18.39765	18.41131	19.06767	19.15229	19.05193
ENSGALG00010002949 ENSGALG00010002949	534836.9775	-0.6803653	0.0267654	-25.41954	0	0	18.66023	18.66924	18.65606	19.31972	19.39478	19.29568
ENSGALG00010009067 ENSGALG00010009067	203805.0933	-0.6980345	0.0278322	-25.08013	0	0	17.30344	17.30927	17.22611	17.96055	17.99866	17.93573
ENSGALG00010018616 ENSGALG00010018616	544278.5844	-0.6720191	0.0288915	-23.26007	0	0	18.65787	18.71084	18.70660	19.34623	19.41905	19.31035
ENSGALG00010009793 TNFAIP6	4464.2580	0.9719472	0.0456664	21.28363	0	0	13.16411	13.15566	13.21884	12.59233	12.54394	12.64832
ENSGALG00010022071 ACTA2	43846.8244	-0.4723792	0.0230126	-20.52696	0	0	15.28501	15.32576	15.29375	15.73850	15.73571	15.73870
ENSGALG00010016647 EGR1	3657.9302	-0.9109528	0.0457658	-19.90467	0	0	12.49625	12.46253	12.44206	12.95232	13.03656	12.95284
ENSGALG00010013345 DCN	5179.6563	0.6776020	0.0342031	19.81110	0	0	13.24818	13.23261	13.27551	12.82257	12.82740	12.81751
ENSGALG00010027337 AQP1	20799.0852	0.5037333	0.0257757	19.54293	0	0	14.75603	14.78298	14.78199	14.33215	14.32762	14.36792
ENSGALG00010005417 LOXL3	771.1802	1.5398981	0.0788635	19.52611	0	0	11.97615	11.98615	11.95391	11.53834	11.54718	11.49336
ENSGALG00010024570 TPM1	45467.9084	-0.4386000	0.0225627	-19.43914	0	0	15.35302	15.37821	15.37045	15.76927	15.78830	15.76219
ENSGALG00010023078 TAGLN	3845.3403	-0.7620142	0.0399556	-19.07152	0	0	12.56167	12.57109	12.52818	13.00323	12.96238	13.01143
ENSGALG00010007529 CDH13	17234.2465	0.5250673	0.0276973	18.95732	0	0	14.52554	14.55656	14.55737	14.10517	14.13993	14.08305
ENSGALG00010024266 ENSGALG00010024266	1239.9974	-1.1639972	0.0685222	-16.98716	0	0	11.79915	11.79426	11.73573	12.23202	12.21638	12.15865
ENSGALG00010018583 PPL	190.1852	3.0313472	0.1817951	16.67453	0	0	11.46364	11.44805	11.52280	11.08591	11.07339	11.09124
ENSGALG00010006546 SMOC2	8575.7594	0.5070699	0.0359863	14.09065	0	0	13.70860	13.69746	13.77144	13.34728	13.32509	13.39715
ENSGALG00010017983 SLC4A11	1553.3044	-0.7385209	0.0526192	-14.03521	0	0	11.96701	11.99134	11.98150	12.28009	12.27571	12.29094
ENSGALG00010024254 ENSGALG00010024254	315.3533	-1.8665545	0.1333077	-14.00185	0	0	11.24508	11.29439	11.22861	11.64451	11.56554	11.59804
ENSGALG00010001966 ARHGAP6	375.6951	1.5312095	0.1105491	13.85094	0	0	11.65129	11.65890	11.63328	11.34965	11.34659	11.29544

dge9 <- dge
write.table(dge9,file="dge9.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge9,"Rfx Tfx ctl vs Hf Tfx ctl")

make_volcano(dge9,"Rfx Tfx ctl vs Hf Tfx ctl")

make_heatmap(dge9,"Rfx Tfx ctl vs Hf Tfx ctl",ss9,xx9,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge9$gene <- sapply(strsplit(rownames(dge9)," "),"[[",2)
m9 <- aggregate(stat ~ gene,dge9,mean)
rownames(m9) <- m9$gene ; m9$gene <- NULL
mres9 <- mitch_calc(m9, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres9$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx Tfx ctl vs Hf Tfx ctl") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx Tfx ctl vs Hf Tfx ctl

	set	setSize	pANOVA	s.dist	p.adjustANOVA
777	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.7836617	0.0000000
1046	GO:0035145 exon-exon junction complex	10	0.0000517	-0.7391745	0.0060881
145	GO:0002181 cytoplasmic translation	27	0.0000000	0.7375536	0.0000000
779	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.7212603	0.0000002
1733	GO:1905564 positive regulation of vascular endothelial cell proliferation	10	0.0004731	0.6383407	0.0298708
1082	GO:0036121 double-stranded DNA helicase activity	11	0.0003015	-0.6292759	0.0231733
217	GO:0004683 calcium/calmodulin-dependent protein kinase activity	11	0.0006656	0.5926115	0.0346096
209	GO:0004602 glutathione peroxidase activity	12	0.0006052	0.5717165	0.0328475
1739	GO:1990391 DNA repair complex	10	0.0035511	-0.5324486	0.0910422
1684	GO:0140662 ATP-dependent protein folding chaperone	23	0.0000117	-0.5278440	0.0020712
429	GO:0006623 protein targeting to vacuole	11	0.0027829	-0.5207864	0.0793580
1076	GO:0035925 mRNA 3’-UTR AU-rich region binding	13	0.0014093	0.5114323	0.0534179
180	GO:0003735 structural constituent of ribosome	118	0.0000000	0.5093871	0.0000000
1758	GO:2000648 positive regulation of stem cell proliferation	12	0.0024527	0.5050266	0.0737111
798	GO:0030127 COPII vesicle coat	10	0.0058750	0.5030850	0.1180342
615	GO:0009378 four-way junction helicase activity	13	0.0018124	-0.4996819	0.0582618
175	GO:0003724 RNA helicase activity	37	0.0000002	-0.4918166	0.0000583
608	GO:0008593 regulation of Notch signaling pathway	10	0.0072324	0.4905364	0.1291627
1742	GO:1990518 single-stranded 3’-5’ DNA helicase activity	11	0.0053756	-0.4847285	0.1138289
1696	GO:0140861 DNA repair-dependent chromatin remodeling	14	0.0017337	-0.4835379	0.0568600

write.table(mres9$enrichment_result,file="mitch9.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres9$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 40

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx Tfx ctl vs Hf Tfx ctl",xlab="ES")
grid()

if ( ! file.exists("mitch9.html") ) {
  mitch_report(res=mres9,outfile="mitch9.html")
}

10. Rfx_NTfx_ctl versus Hf_NTfx_ctl

ss10 <- subset(ss,SampleGroup=="Rfx_NTfx_ctl" | SampleGroup=="Hf_NTfx_ctl")
ss10$trt <- factor(grepl("Hf",ss10$SampleGroup))
rownames(ss10) <- ss10$SampleLabel

xx10 <- xx[,colnames(xx) %in% ss10$SampleLabel]
xx10f <- xx10[rowMeans(xx10)>=10,]
dim(xx10)

## [1] 30111     6

dim(xx10f)

## [1] 14237     6

dds <- DESeqDataSetFromMatrix(countData = xx10f , colData = ss10, design = ~ trt )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes in Hf as compared to Rfx (NTfx ctl)") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes in Hf as compared to Rfx (NTfx ctl)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_NTfx_ctl_1	Hf_NTfx_ctl_2	Hf_NTfx_ctl_3	Rfx_NTfx_ctl_1	Rfx_NTfx_ctl_2	Rfx_NTfx_ctl_3
ENSGALG00010005474 COL8A1	4913.2772	1.2209302	0.0402450	30.33741	0	0	13.10232	13.12808	13.08417	12.20203	12.25901	12.21044
ENSGALG00010009793 TNFAIP6	3492.1439	1.1160675	0.0515550	21.64810	0	0	12.72016	12.69834	12.69551	11.92403	11.92121	12.04246
ENSGALG00010013345 DCN	6085.8214	0.7182705	0.0379330	18.93525	0	0	13.22778	13.17361	13.24523	12.67281	12.65105	12.66953
ENSGALG00010022071 ACTA2	38383.0388	-0.5404750	0.0291876	-18.51727	0	0	15.00890	15.04445	15.01582	15.50711	15.56130	15.54005
ENSGALG00010024266 ENSGALG00010024266	1777.1059	-1.1510012	0.0630063	-18.26804	0	0	11.46499	11.47697	11.37332	12.06177	12.09515	12.00912
ENSGALG00010023078 TAGLN	3380.0513	-0.8167866	0.0447485	-18.25284	0	0	12.06826	12.05156	12.04397	12.55389	12.63734	12.59883
ENSGALG00010027337 AQP1	21189.7609	0.5811939	0.0333915	17.40543	0	0	14.74510	14.72609	14.75695	14.15978	14.25249	14.22204
ENSGALG00010005417 LOXL3	709.4946	1.3398635	0.0797463	16.80157	0	0	11.37818	11.38317	11.39803	10.85610	10.91211	10.90847
ENSGALG00010000588 ENSGALG00010000588	20314.0353	0.7127700	0.0429404	16.59906	0	0	14.63852	14.78564	14.77269	14.06399	14.14292	14.04919
ENSGALG00010024254 ENSGALG00010024254	344.6289	-2.0923129	0.1273135	-16.43434	0	0	10.53561	10.53628	10.47439	10.97980	11.07570	11.07734
ENSGALG00010001707 ENSGALG00010001707	172.7742	-5.9309245	0.3714765	-15.96581	0	0	10.11689	10.07775	10.13334	10.82744	10.81121	10.83053
ENSGALG00010006546 SMOC2	10127.1906	0.6035196	0.0385533	15.65414	0	0	13.82491	13.73471	13.81814	13.25254	13.28246	13.32158
ENSGALG00010004783 ENSGALG00010004783	3764.9784	0.6258509	0.0410313	15.25301	0	0	12.64119	12.66616	12.66422	12.25132	12.22757	12.20279
ENSGALG00010018583 PPL	148.1320	2.9775210	0.2033634	14.64138	0	0	10.74340	10.72310	10.71310	10.30926	10.24258	10.23297
ENSGALG00010011313 RASSF9	2586.0599	0.7225685	0.0496568	14.55124	0	0	12.30650	12.26242	12.35030	11.84396	11.85453	11.88642
ENSGALG00010001966 ARHGAP6	346.0891	1.5194017	0.1107348	13.72109	0	0	11.00521	11.00006	11.00677	10.58600	10.58652	10.63268
ENSGALG00010002566 TMSB4X	27168.3486	0.3973974	0.0298101	13.33099	0	0	15.02359	14.98211	15.01263	14.62974	14.61464	14.66332
ENSGALG00010018867 RPS2	52755.4226	0.4178890	0.0316454	13.20536	0	0	15.97908	15.89776	15.90823	15.50210	15.54748	15.52946
ENSGALG00010007529 CDH13	18511.8221	0.4107463	0.0317293	12.94534	0	0	14.49735	14.46503	14.53089	14.14369	14.13130	14.10502
ENSGALG00010025318 NUAK2	1241.2072	-0.8088486	0.0625057	-12.94040	0	0	11.31167	11.30671	11.31519	11.66072	11.74002	11.68401

dge10 <- dge
write.table(dge10,file="dge10.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge10,"Rfx NTfx ctl vs Hf NTfx ctl")

make_volcano(dge10,"Rfx NTfx ctl vs Hf NTfx ctl")

make_heatmap(dge10,"Rfx NTfx ctl vs Hf NTfx ctl",ss10,xx10,n=30)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge10$gene <- sapply(strsplit(rownames(dge10)," "),"[[",2)
m10 <- aggregate(stat ~ gene,dge10,mean)
rownames(m10) <- m10$gene ; m10$gene <- NULL
mres10 <- mitch_calc(m10, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres10$enrichment_result,20) %>%
  kbl(caption = "Top GO in Rfx NTfx ctl vs Hf NTfx ctl") %>%
  kable_paper("hover", full_width = F)

Top GO in Rfx NTfx ctl vs Hf NTfx ctl

	set	setSize	pANOVA	s.dist	p.adjustANOVA
778	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.8670008	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	0.7528653	0.0000000
1222	GO:0045116 protein neddylation	10	0.0000560	0.7356906	0.0024108
780	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.7293922	0.0000001
1046	GO:0035145 exon-exon junction complex	10	0.0001121	-0.7053744	0.0039546
1422	GO:0051393 alpha-actinin binding	10	0.0001289	-0.6991054	0.0044601
1320	GO:0048041 focal adhesion assembly	14	0.0000252	-0.6502703	0.0012205
1082	GO:0036121 double-stranded DNA helicase activity	11	0.0003045	-0.6288199	0.0082647
180	GO:0003735 structural constituent of ribosome	118	0.0000000	0.6260841	0.0000000
973	GO:0032585 multivesicular body membrane	12	0.0002988	0.6028883	0.0082364
390	GO:0006337 nucleosome disassembly	10	0.0009718	-0.6023808	0.0200929
893	GO:0030992 intraciliary transport particle B	15	0.0000999	0.5802048	0.0036723
32	GO:0000380 alternative mRNA splicing, via spliceosome	16	0.0000890	-0.5658604	0.0034112
1682	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0014852	-0.5533057	0.0256858
1487	GO:0060395 SMAD protein signal transduction	14	0.0003811	-0.5484092	0.0097434
1107	GO:0042073 intraciliary transport	16	0.0001646	0.5440684	0.0054792
155	GO:0003148 outflow tract septum morphogenesis	12	0.0013319	-0.5350123	0.0237320
375	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000727	0.5257519	0.0029804
56	GO:0000793 condensed chromosome	13	0.0010318	-0.5256826	0.0202529
1537	GO:0070411 I-SMAD binding	11	0.0026225	-0.5239376	0.0392035

write.table(mres10$enrichment_result,file="mitch10.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres10$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 131

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Rfx NTfx ctl vs Hf NTfx ctl",xlab="ES")
grid()

par(mar = c(5.1, 4.1, 4.1, 2.1) )

if ( ! file.exists("mitch10.html") ) {
  mitch_report(res=mres10,outfile="mitch10.html")
}

11. Interaction between cell line and target gRNA

ss11 <- subset(ss,SampleGroup=="Rfx_NTC_gRNA" | SampleGroup=="Rfx_DsRed_gRNA" | SampleGroup=="Hf_NTC_gRNA" | SampleGroup=="Hf_DsRed_gRNA")
ss11$tgt <- factor(grepl("DsRed",ss11$SampleGroup))
ss11$hf <- factor(grepl("Hf",ss11$SampleGroup))

rownames(ss11) <- ss11$SampleLabel

xx11 <- xx[,colnames(xx) %in% ss11$SampleLabel]
xx11f <- xx11[rowMeans(xx11)>=10,]
dim(xx11)

## [1] 30111    12

dim(xx11f)

## [1] 14471    12

dds <- DESeqDataSetFromMatrix(countData = xx11f , colData = ss11, design = ~ tgt * hf )

## 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),])

head(dge,20) %>%
  kbl(caption = "Top gene expression changes caused by interaction of DsRed targeting and  Rfx vs Hf") %>%
  kable_paper("hover", full_width = F)

Top gene expression changes caused by interaction of DsRed targeting and Rfx vs Hf

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	Hf_DsRed_gRNA_1	Hf_DsRed_gRNA_2	Hf_DsRed_gRNA_3	Hf_NTC_gRNA_1	Hf_NTC_gRNA_2	Hf_NTC_gRNA_3	Rfx_DsRed_gRNA_1	Rfx_DsRed_gRNA_2	Rfx_DsRed_gRNA_3	Rfx_NTC_gRNA_1	Rfx_NTC_gRNA_2	Rfx_NTC_gRNA_3
DsRed DsRed	8066.52333	1.3020176	0.0869722	14.970495	0.0000000	0.0000000	12.197207	12.270675	12.187288	13.632701	13.460708	13.487294	11.606560	11.567849	11.731227	14.109481	13.933426	14.015058
ENSGALG00010017057 NDUFAB1	1951.86547	0.4325047	0.0922153	4.690164	0.0000027	0.0116847	11.717552	11.776832	11.767789	11.530794	11.487702	11.527004	11.249755	11.376986	11.424503	11.360960	11.387851	11.497924
ENSGALG00010020691 ENSGALG00010020691	817.01552	0.5469877	0.1168280	4.681992	0.0000028	0.0116847	10.832841	10.873360	10.982117	10.770724	10.700098	10.705296	10.556225	10.669116	10.675794	10.719825	10.752426	10.780337
ENSGALG00010023749 RPP25	91.34135	1.5552370	0.3371352	4.613096	0.0000040	0.0116847	9.845658	9.801078	9.792541	9.703250	9.522536	9.632288	9.483165	9.564841	9.558902	9.574391	9.653987	9.739156
ENSGALG00010024763 HOXD9	32.18087	3.0220335	0.6556184	4.609440	0.0000040	0.0116847	9.658433	9.606746	9.644846	9.431427	9.264160	9.358889	9.232229	9.192004	9.325291	9.238688	9.363163	9.393097
ENSGALG00010017499 SNRPD1	3744.21816	0.3055852	0.0680723	4.489125	0.0000072	0.0167117	12.308570	12.385796	12.337484	12.206192	12.190173	12.202358	12.056307	12.132723	12.137909	12.138153	12.217252	12.268531
ENSGALG00010003570 ENSGALG00010003570	669.33628	-6.7238577	1.5065997	-4.462936	0.0000081	0.0167117	9.416812	9.406989	9.399598	9.933820	11.417810	9.846013	11.800665	11.280347	10.567840	10.617005	10.279986	9.482196
ENSGALG00010027900 REN	876.48513	-0.4530758	0.1084750	-4.176776	0.0000296	0.0534792	10.723937	10.643288	10.661910	10.917137	10.815309	10.818966	10.908402	10.853193	10.888229	10.765858	10.889729	10.794441
ENSGALG00010009059 MED29	667.25990	0.4635630	0.1125726	4.117903	0.0000382	0.0614711	10.699818	10.733586	10.792496	10.604665	10.570080	10.578123	10.471373	10.495802	10.567840	10.557231	10.562018	10.629159
ENSGALG00010001531 ENSGALG00010001531	38.84716	-2.1118212	0.5215305	-4.049277	0.0000514	0.0743414	9.411516	9.489638	9.387984	9.524114	9.559058	9.366532	9.561388	9.474553	9.614114	9.296499	9.343141	9.399824
ENSGALG00010019231 ENSGALG00010019231	860.66187	-0.4775430	0.1250308	-3.819402	0.0001338	0.1548681	10.737999	10.640221	10.634727	10.763085	10.933171	10.784660	10.810046	10.872829	10.962899	10.765858	10.781344	10.818239
ENSGALG00010023836 NDUFS3	2082.62673	0.2826528	0.0741048	3.814233	0.0001366	0.1548681	11.645296	11.736136	11.733682	11.551201	11.579437	11.573240	11.464507	11.525416	11.494658	11.563816	11.507895	11.596568
ENSGALG00010021073 NDUFB8	2906.00702	0.3462428	0.0913499	3.790293	0.0001505	0.1548681	11.972696	12.120472	12.111108	11.936143	11.854028	11.944728	11.706600	11.878073	11.864494	11.905521	11.884076	11.980338
ENSGALG00010020607 ENSGALG00010020607	2050.30660	0.3133224	0.0828182	3.783256	0.0001548	0.1548681	11.766191	11.738308	11.736176	11.568125	11.484085	11.617413	11.379977	11.471153	11.518040	11.467223	11.487623	11.514587
ENSGALG00010023206 HSPB9	290.56821	0.6822730	0.1807743	3.774171	0.0001605	0.1548681	9.915701	9.931631	9.865365	10.181759	10.203777	10.164554	9.863152	9.867996	9.864787	10.512639	10.306989	10.357217
ENSGALG00010020624 MINOS1	1460.06212	0.3039089	0.0815506	3.726631	0.0001941	0.1677978	11.268463	11.352899	11.331250	11.245348	11.238210	11.198956	11.112868	11.209838	11.143517	11.278166	11.240361	11.256464
ENSGALG00010007964 RECQL4	1055.76080	-0.4017280	0.1079143	-3.722658	0.0001971	0.1677978	10.862338	10.798819	10.857670	10.949696	10.950461	10.963268	11.168024	11.045682	11.000023	11.036413	10.917126	10.913564
ENSGALG00010025365 SIKE1	3887.99423	0.2847988	0.0780400	3.649393	0.0002629	0.2113110	12.250367	12.345370	12.425101	12.183367	12.142894	12.220531	12.189145	12.222402	12.257848	12.258585	12.250038	12.366356
ENSGALG00010003549 ENSGALG00010003549	1027.11817	1.9886030	0.5502423	3.614050	0.0003015	0.2295790	11.350463	11.548611	11.597743	10.958430	10.600865	10.180898	10.696008	10.611555	10.353940	10.824314	10.569321	11.006260
ENSGALG00010023617 ENSGALG00010023617	186.29156	0.9762562	0.2719689	3.589588	0.0003312	0.2396242	10.035607	10.012714	10.059341	9.935786	9.757533	10.017433	9.702894	9.759882	9.840549	9.819975	9.974632	9.941505

dge11 <- dge
write.table(dge11,file="dge11.tsv",quote=FALSE,sep='\t')

# plots
par(mar = c(5.1, 4.1, 4.1, 2.1) )
maplot(dge11,"gRNA target * Rfx/HF interaction")

make_volcano(dge11,"gRNA target * Rfx/HF interaction")

make_heatmap(dge11,"gRNA target * Rfx/HF interaction",ss11,xx11,n=20)

## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf

## Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

# pathway
dge11$gene <- sapply(strsplit(rownames(dge11)," "),"[[",2)
m11 <- aggregate(stat ~ gene,dge11,mean)
rownames(m11) <- m11$gene ; m11$gene <- NULL
mres11 <- mitch_calc(m11, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres11$enrichment_result,20) %>%
  kbl(caption = "Top GO in gRNA target * Rfx/Hf interaction") %>%
  kable_paper("hover", full_width = F)

Top GO in gRNA target * Rfx/Hf interaction

	set	setSize	pANOVA	s.dist	p.adjustANOVA
778	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.8982640	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	0.8508806	0.0000000
780	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.8118646	0.0000000
180	GO:0003735 structural constituent of ribosome	118	0.0000000	0.7504334	0.0000000
1237	GO:0045271 respiratory chain complex I	31	0.0000000	0.6687426	0.0000000
408	GO:0006412 translation	96	0.0000000	0.6627106	0.0000000
440	GO:0006749 glutathione metabolic process	15	0.0000114	0.6544757	0.0006710
994	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	0.6445550	0.0000000
376	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000015	0.6368706	0.0001181
315	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	0.6263827	0.0000000
1130	GO:0042474 middle ear morphogenesis	10	0.0008469	0.6093964	0.0220313
110	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0002570	0.6093631	0.0089141
1018	GO:0034314 Arp2/3 complex-mediated actin nucleation	11	0.0005318	0.6032004	0.0156783
1527	GO:0070131 positive regulation of mitochondrial translation	10	0.0013457	0.5854896	0.0313989
1685	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0007949	-0.5841114	0.0214273
1741	GO:1990391 DNA repair complex	10	0.0015341	-0.5785739	0.0343517
572	GO:0008137 NADH dehydrogenase (ubiquinone) activity	10	0.0017401	0.5718523	0.0370233
347	GO:0005852 eukaryotic translation initiation factor 3 complex	10	0.0020126	0.5640080	0.0395580
806	GO:0030150 protein import into mitochondrial matrix	12	0.0007243	0.5635410	0.0200197
200	GO:0004364 glutathione transferase activity	10	0.0020656	0.5625944	0.0398510

write.table(mres11$enrichment_result,file="mitch11.tsv",quote=FALSE,sep='\t')

par(mar=c(5,25,3,3))
top <- mres11$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)

## [1] 110

up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.7,
  main="Top GO in gRNA target * Rfx/Hf interaction",xlab="ES")
grid()

par(mar = c(5.1, 4.1, 4.1, 2.1) )

if ( ! file.exists("mitch11.html") ) {
  mitch_report(res=mres11,outfile="mitch11.html")
}

Cross comparison

Make a Euler diagram of the differentially expressed genes.

Specifically compare contrast 1 and 4.

up1 <- unique(sapply(strsplit(rownames(subset(dge1,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
dn1 <- unique(sapply(strsplit(rownames(subset(dge1,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))

up4 <- unique(sapply(strsplit(rownames(subset(dge4,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
dn4 <- unique(sapply(strsplit(rownames(subset(dge4,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))

v1 <- list("Rfx up"=up1,"Rfx  dn"=dn1,
  "Hf up"=up4,"Hf dn"=dn4)

plot(euler(v1), quantities = TRUE, main="NTC vs DsRed gRNA")

message("Common up genes")

## Common up genes

intersect(up1,up4)

##  [1] "IL8L2"              "TNIP1"              "K123"              
##  [4] "PTX3"               "ENSGALG00010015706" "CCL4"              
##  [7] "ENSGALG00010022322" "IGFBP5"             "CSF3"              
## [10] "ENSGALG00010022313" "TNFSF15"            "MYC"

message("Common down genes")

## Common down genes

intersect(dn1,dn4)

## [1] "DsRed"      "Fbxw11-GFP" "HSPB9"      "JAM3"       "WNT9A"

message("Up genes specific to Rfx")

## Up genes specific to Rfx

setdiff(up1,up4)

##  [1] "ENSGALG00010019621" "NUAK2"              "ADAM8"             
##  [4] "BOP1"               "SQSTM1"             "TMEM63B"           
##  [7] "PPRC1"              "NFKB2"              "PHF19"             
## [10] "PIEZO1"             "UTP4"               "ENSGALG00010029462"
## [13] "NOP56"              "ENSGALG00010012858" "ENSGALG00010000177"
## [16] "SLC8A1"             "ENSGALG00010027165" "BRD2"              
## [19] "TNFRSF10B"

message("Up genes specific to Hf")

## Up genes specific to Hf

setdiff(up4,up1)

##  [1] "ENSGALG00010003545" "NDUFAB1"            "IRX5"              
##  [4] "HOXD9"              "ENSGALG00010006816" "NFKBIA"            
##  [7] "ENSGALG00010025570" "RRS1"               "TNFAIP6"           
## [10] "ENSGALG00010020607" "AVD"                "EIF1"              
## [13] "AVPI1"              "TNFRSF11B"          "JUND"              
## [16] "MRPL17"             "CX3CL1"             "RPP25"             
## [19] "ENSGALG00010024292" "MED29"              "MRPS15"            
## [22] "ENSGALG00010020691" "SNRPD1"             "ENSGALG00010000278"
## [25] "MRPS28"             "FAM43A"             "ENSGALG00010009123"
## [28] "WBP1"               "DTX2"               "FABP7"             
## [31] "RPS28"              "GATAD1"             "EXOSC6"            
## [34] "ENSGALG00010028461" "SIKE1"              "FAM212A"           
## [37] "NDUFS3"             "SLC35B1"            "JAC1"              
## [40] "NDUFB8"             "ID4"                "DIO3"

message("Down genes specific to Rfx")

## Down genes specific to Rfx

setdiff(dn1,dn4)

##  [1] "SCP2"               "APCDD1"             "RAB23"             
##  [4] "C1orf198"           "TGFBR3"             "CPED1"             
##  [7] "C1QTNF3"            "SASH1"              "WBP1L"             
## [10] "ENSGALG00010028177" "MXRA8"              "TRANK1"            
## [13] "PTBP3"

message("Down genes specific to Hf")

## Down genes specific to Hf

setdiff(dn4,dn1)

##   [1] "ENSGALG00010012473" "LPP"                "TNRC6B"            
##   [4] "LAMA5"              "ENSGALG00010018888" "SGCD"              
##   [7] "ENSGALG00010018696" "ENSGALG00010001704" "TCF4"              
##  [10] "ENSGALG00010012687" "ENSGALG00010016594" "ENSGALG00010012677"
##  [13] "GLI2"               "ASH1L"              "ENSGALG00010018886"
##  [16] "SH3PXD2A"           "TTC28"              "MGA"               
##  [19] "ZBTB20"             "SPEN"               "GOLGB1"            
##  [22] "CDC42BPA"           "KMT2D"              "EP300"             
##  [25] "CEP350"             "AKAP13"             "ENSGALG00010001318"
##  [28] "FASN"               "PRRC2B"             "COL6A3"            
##  [31] "ENSGALG00010004063" "MYO9A"              "ENSGALG00010000697"
##  [34] "PCNT"               "FBN1"               "PRRC2C"            
##  [37] "ZNF469"             "PLEKHG4"            "ENSGALG00010009160"
##  [40] "FLNB"               "ADAMTS9"            "ENSGALG00010020265"
##  [43] "EP400"              "NAV2"               "SPTAN1"            
##  [46] "ARHGAP31"           "MYH11"              "NOTCH2"            
##  [49] "ANKRD17"            "CHD6"               "ENSGALG00010007190"
##  [52] "IGF1R"              "FRYL"               "ENSGALG00010008225"
##  [55] "REN"                "FBN2"               "ENSGALG00010013957"
##  [58] "NOTUM"              "ASPM"               "ENSGALG00010012478"
##  [61] "AGO2"               "KIF1B"              "ARHGEF12"          
##  [64] "MYO5A"              "NOTCH1"             "ENSGALG00010023327"
##  [67] "TNC"                "USP34"              "HOOK3"             
##  [70] "NCOR1"              "NBEAL1"             "ENSGALG00010004662"
##  [73] "ENSGALG00010025824" "ENSGALG00010010289" "CEP250"            
##  [76] "TNRC18"             "CHD7"               "PEAK1"             
##  [79] "ENSGALG00010028926" "ABCA5"              "AXIN2"             
##  [82] "ENSGALG00010019437" "BOD1L1"             "SLK"               
##  [85] "CHD9"               "ENSGALG00010013747" "ENSGALG00010018018"
##  [88] "COL5A1"             "NF1"                "ENSGALG00010012044"
##  [91] "RC3H1"              "ENSGALG00010008220" "POU2F1"            
##  [94] "TANC2"              "ZZEF1"              "ENSGALG00010026717"
##  [97] "TNS1"               "DOT1L"              "GAS7"              
## [100] "MYH10"              "ENSGALG00010024799" "ENSGALG00010009880"
## [103] "N4BP2"              "MED13"              "ENSGALG00010002289"
## [106] "CENPF"              "COL6A1"             "AKAP11"            
## [109] "ENSGALG00010007324" "DGKH"               "CCDC88A"           
## [112] "RALGAPB"            "CREBBP"             "SYNE2"             
## [115] "XRN1"               "CKAP5"              "NAV3"              
## [118] "VPS13A"             "PLCE1"              "SPTBN1"            
## [121] "ENSGALG00010013811" "SETD1A"             "DNAJC13"           
## [124] "BMPR2"              "ENSGALG00010019440" "ENSGALG00010010124"
## [127] "ENSGALG00010019735" "OIP5-AS1"           "NCOA3"             
## [130] "MDN1"               "MACF1"              "ENSGALG00010023657"
## [133] "EIF4EBP3"           "PHIP"               "BRWD3"             
## [136] "SYNM"               "KMT2C"              "APC"               
## [139] "DNAH10"             "UBR4"

up1 <- unique(sapply(strsplit(rownames(subset(dge1,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up2 <- unique(sapply(strsplit(rownames(subset(dge2,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up3 <- unique(sapply(strsplit(rownames(subset(dge3,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up4 <- unique(sapply(strsplit(rownames(subset(dge4,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up5 <- unique(sapply(strsplit(rownames(subset(dge5,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up6 <- unique(sapply(strsplit(rownames(subset(dge6,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up7 <- unique(sapply(strsplit(rownames(subset(dge7,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up8 <- unique(sapply(strsplit(rownames(subset(dge8,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up9 <- unique(sapply(strsplit(rownames(subset(dge9,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))
up10 <- unique(sapply(strsplit(rownames(subset(dge10,padj < 0.05 & log2FoldChange > 0))," "),"[[",2))

dn1 <- unique(sapply(strsplit(rownames(subset(dge1,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn2 <- unique(sapply(strsplit(rownames(subset(dge2,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn3 <- unique(sapply(strsplit(rownames(subset(dge3,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn4 <- unique(sapply(strsplit(rownames(subset(dge4,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn5 <- unique(sapply(strsplit(rownames(subset(dge5,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn6 <- unique(sapply(strsplit(rownames(subset(dge6,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn7 <- unique(sapply(strsplit(rownames(subset(dge7,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn8 <- unique(sapply(strsplit(rownames(subset(dge8,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn9 <- unique(sapply(strsplit(rownames(subset(dge9,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))
dn10 <- unique(sapply(strsplit(rownames(subset(dge10,padj < 0.05 & log2FoldChange < 0))," "),"[[",2))


v1 <- list("Up1"=up1,"Up2"=up2,"Up3"=up3,
  "Dn1"=dn1,"Dn2"=dn2,"Dn3"=dn3)

plot(euler(v1), quantities = TRUE, main="Rfx gene overlap")

v1 <- list("Up4"=up4,"Up5"=up5,"Up6"=up6,
  "Dn4"=dn4,"Dn5"=dn5,"Dn6"=dn6)

plot(euler(v1), quantities = TRUE, main="Hf gene overlap")

v1 <- list("Up7"=up7,"Up8"=up8,"Up9"=up9,"Up10"=up10,
  "Dn7"=dn7,"Dn8"=dn8,"Dn9"=dn9,"Dn10"=dn10)

plot(euler(v1), quantities = TRUE, main="Rfx-Hf comparison gene overlap")

v1 <- list("Up1"=up1,"Up2"=up2,"Up3"=up3, "Up4"=up4,"Up5"=up5,"Up6"=up6,"Up7"=up7,"Up8"=up8,"Up9"=up9,"Up10"=up10,
  "Dn1"=dn1,"Dn2"=dn2,"Dn3"=dn3,"Dn4"=dn4,"Dn5"=dn5,"Dn6"=dn6,"Dn7"=dn7,"Dn8"=dn8,"Dn9"=dn9,"Dn10"=dn10)

upset(fromList(v1), order.by = "freq")

Scatterplot of NTC vs DsRed gRNA

l <- list("Rfx"=m1, "Hf"=m4)
m <- mitch_import(l, DEtype="prescored")

## Note: Mean no. genes in input = 14456

## Note: no. genes in output = 14101

## Note: estimated proportion of input genes in output = 0.975

head(m)

##                  Rfx          Hf
## 5_8S_rRNA  1.8503649 -0.01023976
## A4GALT    -1.8915624  0.23351781
## AAAS       2.1906313  0.21120150
## AACS       1.1170591  0.57504626
## AADAC      0.7635603 -0.09350280
## AADAT     -0.2696192 -1.21123488

plot(m,ylim=c(-50,10),xlim=c(-50,10),main="Test stat")
abline(1,1,col="red",lty=2,lwd=2)
sig <- subset(m,Rfx < -10)
points(sig,pch=19)
text(sig$Rfx,sig$Hf-2,labels=rownames(sig))

sig

##                  Rfx        Hf
## DsRed      -47.99253 -28.00853
## Fbxw11-GFP -17.80101 -23.77922
## HSPB9      -11.16618  -6.96703

plot(m,ylim=c(-11,10),xlim=c(-11,10),main="Test stat")
abline(1,1,col="red",lty=2,lwd=2)
sig <- subset(m,Hf > 5 | Rfx > 5 | Hf < -5 | Rfx < -5)
points(sig,pch=19)
text(sig$Rfx,sig$Hf+0.5,labels=rownames(sig),cex=0.7)

plot(m,ylim=c(-2,11),xlim=c(-2,11),main="Test stat")
abline(1,1,col="red",lty=2,lwd=2)
sig <- subset(m,Hf > 5 | Rfx > 5)
points(sig,pch=19)
text(sig$Rfx,sig$Hf+0.5,labels=rownames(sig),cex=0.7)

sig

##                           Rfx         Hf
## CCL4                4.9146227  7.2069917
## CSF3                4.5616423  5.9868265
## ENSGALG00010003545 -0.7941936  6.7986079
## ENSGALG00010006816  2.8476270  5.9996088
## ENSGALG00010015706  5.1313522  3.6962549
## ENSGALG00010019621  5.4585553  0.7078594
## ENSGALG00010022322  4.7800917 10.0182679
## ENSGALG00010025570  1.1389740  5.2588473
## HOXD9              -1.6078076  6.0188043
## IL8L2               7.4032653  5.2112555
## IRX5                0.9078357  6.1598455
## K123                5.9266651  4.3901759
## NDUFAB1            -1.2928359  6.1785418
## NFKBIA              2.7841672  5.5024241
## NUAK2               5.1175934  1.1493253
## PTX3                5.1690275  4.3454939
## TNIP1               7.0925819  4.9665784

plot(m,ylim=c(-11,0),xlim=c(-11,0),main="Test stat")
abline(1,1,col="red",lty=2,lwd=2)
sig <- subset(m,Hf < -5 | Rfx < -5)
points(sig,pch=19)
text(sig$Rfx,sig$Hf-0.5,labels=rownames(sig),cex=0.7)

sig

##                             Rfx         Hf
## ASH1L               -0.31840850  -5.263332
## DsRed              -47.99252613 -28.008533
## ENSGALG00010001704  -0.50218477  -5.461451
## ENSGALG00010012473   0.07297691  -6.766306
## ENSGALG00010012677  -0.61319658  -5.375869
## ENSGALG00010012687  -0.61848790  -5.418866
## ENSGALG00010016594  -0.81309349  -5.397298
## ENSGALG00010018696  -0.57295985  -5.542007
## ENSGALG00010018886  -0.45556734  -5.262599
## ENSGALG00010018888  -0.57261510  -5.934879
## Fbxw11-GFP         -17.80100530 -23.779219
## GLI2                -1.61966555  -5.286938
## GOLGB1              -0.50227317  -5.061648
## HSPB9              -11.16617730  -6.967030
## JAM3                -5.53574121  -4.449473
## LAMA5               -0.64694229  -6.194777
## LPP                 -0.84694703  -6.591509
## MGA                 -0.50704251  -5.197213
## SCP2                -5.04200145  -1.220863
## SGCD                -2.63371796  -5.752621
## SH3PXD2A            -1.64903446  -5.241714
## SPEN                -0.01313174  -5.111272
## TCF4                -0.39153844  -5.419533
## TNRC6B              -0.13518674  -6.448490
## TTC28               -1.35069465  -5.217414
## ZBTB20              -1.15314382  -5.161435

Multi-dimensional pathway enrichment

This is a unique feature of mitch package that allows us to look at enrichment in two or more contrasts.

First look at NTC vs DsRed gRNA in both Rfx and HF.

m0 <- mitch_calc(m, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(m0$enrichment_result,20) %>%
  kbl(caption = "Top pathways for gRNA in Rfx & Hf") %>%
  kable_paper("hover", full_width = F)

Top pathways for gRNA in Rfx & Hf

	set	setSize	pMANOVA	s.Rfx	s.Hf	p.Rfx	p.Hf	s.dist	SD	p.adjustMANOVA
772	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.6379063	0.8780030	0.0000000	0.0000000	1.0852712	1.0719097	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.3703914	0.9165417	0.0008663	0.0000000	0.9885537	0.9099991	0.0000000
774	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.4967036	0.8227678	0.0000172	0.0000000	0.9610730	0.9330072	0.0000000
180	GO:0003735 structural constituent of ribosome	118	0.0000000	-0.4185603	0.7655773	0.0000000	0.0000000	0.8725259	0.8373118	0.0000000
434	GO:0006695 cholesterol biosynthetic process	16	0.0000000	0.5782570	0.6513578	0.0000621	0.0000064	0.8710041	0.0516901	0.0000022
110	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0000072	0.0428467	0.8112949	0.7972188	0.0000011	0.8124255	0.5433749	0.0003137
345	GO:0005852 eukaryotic translation initiation factor 3 complex	10	0.0000960	0.0647505	0.7850117	0.7229665	0.0000171	0.7876776	0.5093016	0.0029063
406	GO:0006412 translation	96	0.0000000	-0.3422498	0.6883464	0.0000000	0.0000000	0.7687364	0.7287416	0.0000000
1223	GO:0045271 respiratory chain complex I	31	0.0000000	-0.4369672	0.6004677	0.0000255	0.0000000	0.7426317	0.7335773	0.0000000
374	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000002	-0.4761883	0.5361978	0.0003266	0.0000520	0.7171216	0.7158651	0.0000131
988	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	-0.3500900	0.6231320	0.0001888	0.0000000	0.7147423	0.6881719	0.0000000
314	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	-0.2319974	0.6620995	0.0133698	0.0000000	0.7015686	0.6322220	0.0000000
437	GO:0006749 glutathione metabolic process	15	0.0000116	-0.4675375	0.5116002	0.0017189	0.0006025	0.6930556	0.6923549	0.0004717
343	GO:0005840 ribosome	43	0.0000000	-0.1888191	0.6661340	0.0322511	0.0000000	0.6923779	0.6045431	0.0000000
1093	GO:0042026 protein refolding	12	0.0003047	0.3522488	0.5876925	0.0346373	0.0004233	0.6851728	0.1664839	0.0074024
694	GO:0016282 eukaryotic 43S preinitiation complex	14	0.0000655	0.1250494	0.6712065	0.4179700	0.0000137	0.6827558	0.3861913	0.0022925
996	GO:0033290 eukaryotic 48S preinitiation complex	13	0.0001751	0.0892413	0.6638274	0.5775218	0.0000340	0.6697991	0.4062937	0.0046392
1011	GO:0034314 Arp2/3 complex-mediated actin nucleation	11	0.0005276	-0.4657333	0.4668688	0.0074842	0.0073399	0.6594498	0.6594493	0.0109864
184	GO:0003774 cytoskeletal motor activity	12	0.0005764	-0.0552559	-0.6434689	0.7403685	0.0001135	0.6458370	0.4159294	0.0115275
266	GO:0005201 extracellular matrix structural constituent	12	0.0005891	-0.0379137	-0.6429247	0.8201394	0.0001150	0.6440417	0.4278074	0.0115275

top <- head(subset(m0$enrichment_result,p.adjustMANOVA<0.05),30)
mycols <- grep("^s\\.",colnames(top))
mycols <- mycols[1:length(mycols)-1]
mx <- as.matrix(top[,mycols])
rownames(mx) <- top$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mx,trace="none",col=colfunc(25),scale="none",
    margins = c(7,25), cexRow=0.8, cexCol=0.8,
    main="Top pathways for gRNA in Rfx & Hf" )

par(mar = c(5.1, 4.1, 4.1, 2.1) )

plot(m0$enrichment_result$s.Rfx,m0$enrichment_result$s.Hf,
  col=alpha("darkgray", 0.6), pch=19,
  main="Enrichment score" , xlab="Rfx", ylab="Hf")

sig <- subset(m0$enrichment_result,p.adjustMANOVA<0.05)

points(sig$s.Rfx,sig$s.Hf, col="darkorange", pch=19 )

abline(v=0,h=0,col="blue",lwd=2,lty=2)

if ( ! file.exists("mitch0.html") ) {
  mitch_report(res=m0,outfile="mitch0.html")
}

# SD=discordant
m0 <- mitch_calc(m, go, priority="SD",cores=8)

## Note: Prioritisation by SD after selecting sets with 
##             p.adjustMANOVA<=0.05.

head(m0$enrichment_result,20) %>%
  kbl(caption = "Discordant pathways for gRNA in Rfx & Hf") %>%
  kable_paper("hover", full_width = F)

Discordant pathways for gRNA in Rfx & Hf

	set	setSize	pMANOVA	s.Rfx	s.Hf	p.Rfx	p.Hf	s.dist	SD	p.adjustMANOVA
772	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.6379063	0.8780030	0.0000000	0.0000000	1.0852712	1.0719097	0.0000000
774	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.4967036	0.8227678	0.0000172	0.0000000	0.9610730	0.9330072	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.3703914	0.9165417	0.0008663	0.0000000	0.9885537	0.9099991	0.0000000
180	GO:0003735 structural constituent of ribosome	118	0.0000000	-0.4185603	0.7655773	0.0000000	0.0000000	0.8725259	0.8373118	0.0000000
1223	GO:0045271 respiratory chain complex I	31	0.0000000	-0.4369672	0.6004677	0.0000255	0.0000000	0.7426317	0.7335773	0.0000000
406	GO:0006412 translation	96	0.0000000	-0.3422498	0.6883464	0.0000000	0.0000000	0.7687364	0.7287416	0.0000000
374	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000002	-0.4761883	0.5361978	0.0003266	0.0000520	0.7171216	0.7158651	0.0000131
437	GO:0006749 glutathione metabolic process	15	0.0000116	-0.4675375	0.5116002	0.0017189	0.0006025	0.6930556	0.6923549	0.0004717
988	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	-0.3500900	0.6231320	0.0001888	0.0000000	0.7147423	0.6881719	0.0000000
1011	GO:0034314 Arp2/3 complex-mediated actin nucleation	11	0.0005276	-0.4657333	0.4668688	0.0074842	0.0073399	0.6594498	0.6594493	0.0109864
569	GO:0008137 NADH dehydrogenase (ubiquinone) activity	10	0.0014909	-0.3658363	0.5291037	0.0451749	0.0037655	0.6432627	0.6328182	0.0230760
314	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	-0.2319974	0.6620995	0.0133698	0.0000000	0.7015686	0.6322220	0.0000000
343	GO:0005840 ribosome	43	0.0000000	-0.1888191	0.6661340	0.0322511	0.0000000	0.6923779	0.6045431	0.0000000
1721	GO:1990391 DNA repair complex	10	0.0033375	0.4307998	-0.4193883	0.0183354	0.0216598	0.6012279	0.6011737	0.0399814
1451	GO:0060047 heart contraction	10	0.0038236	-0.4609751	0.3751047	0.0116018	0.0399974	0.5943076	0.5911977	0.0429174
1422	GO:0051603 proteolysis involved in protein catabolic process	22	0.0000064	-0.4001511	0.4299246	0.0011596	0.0004822	0.5873296	0.5869522	0.0002850
1433	GO:0051881 regulation of mitochondrial membrane potential	14	0.0005492	-0.3412366	0.4736383	0.0270800	0.0021532	0.5837599	0.5762036	0.0113010
200	GO:0004364 glutathione transferase activity	10	0.0034686	-0.2750976	0.5353062	0.1320250	0.0033772	0.6018567	0.5730420	0.0405111
1665	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0010111	0.1510678	-0.6205045	0.3857149	0.0003658	0.6386293	0.5455840	0.0178627
110	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0000072	0.0428467	0.8112949	0.7972188	0.0000011	0.8124255	0.5433749	0.0003137

top <- head(subset(m0$enrichment_result,p.adjustMANOVA<0.05),30)
mycols <- grep("^s\\.",colnames(top))
mycols <- mycols[1:length(mycols)-1]
mx <- as.matrix(top[,mycols])
rownames(mx) <- top$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mx,trace="none",col=colfunc(25),scale="none",
    margins = c(7,25), cexRow=0.8, cexCol=0.8,
    main="Discordant pathways for gRNA in Rfx & Hf" )

Next look at Rfx.

l <- list("NTC v DsRed gRNA"=m1,
  "Tfx v NTC"=m2,
  "NTfc v Tfx"=m3)

mm <- mitch_import(l, DEtype="prescored")

## Note: Mean no. genes in input = 14497.3333333333

## Note: no. genes in output = 14182

## Note: estimated proportion of input genes in output = 0.978

head(mm)

##           NTC v DsRed gRNA  Tfx v NTC  NTfc v Tfx
## 5_8S_rRNA        1.8503649  0.5180755 -1.41779880
## A4GALT          -1.8915624 -0.4952824 -1.93247212
## AAAS             2.1906313 -1.7533504  0.03040747
## AACS             1.1170591  0.8714301  4.57047093
## AADAC            0.7635603 -0.4594811  1.06088899
## AADAT           -0.2696192 -0.0656957  0.66719319

mmres_rfx <- mitch_calc(mm, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mmres_rfx$enrichment_result,20) %>%
  kbl(caption = "Top pathways for Rfx contrasts") %>%
  kable_paper("hover", full_width = F)

Top pathways for Rfx contrasts

	set	setSize	pMANOVA	s.NTC.v.DsRed.gR	s.Tfx.v.NTC	s.NTfc.v.Tfx	p.NTC.v.DsRed.gR	p.Tfx.v.NTC	p.NTfc.v.Tfx	s.dist	SD	p.adjustMANOVA
1681	GO:0140662 ATP-dependent protein folding chaperone	23	0.0000000	0.4263351	-0.6169897	0.7485145	0.0004016	0.0000003	0.0000000	1.0595810	0.7137843	0.0000000
774	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.6377351	0.3625760	-0.7202517	0.0000000	0.0000728	0.0000000	1.0280710	0.6027641	0.0000000
776	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.4962944	0.3652158	-0.7353620	0.0000174	0.0015758	0.0000000	0.9593998	0.5788818	0.0000001
1043	GO:0035145 exon-exon junction complex	10	0.0003118	0.5649591	-0.2761784	0.7184025	0.0019777	0.1305121	0.0000835	0.9547541	0.5354512	0.0057320
1277	GO:0045943 positive regulation of transcription by RNA polymerase I	11	0.0000100	0.3806814	-0.1145425	0.8594826	0.0288162	0.5107384	0.0000008	0.9469681	0.4870356	0.0005219
1515	GO:0070034 telomerase RNA binding	12	0.0001324	0.5173489	-0.3461656	0.7019642	0.0019155	0.0378857	0.0000254	0.9382080	0.5595111	0.0031581
163	GO:0003678 DNA helicase activity	20	0.0000003	0.3801017	-0.4870004	0.6868804	0.0032571	0.0001631	0.0000001	0.9238243	0.6088206	0.0000220
250	GO:0005049 nuclear export signal receptor activity	10	0.0001316	0.4144087	-0.1476997	0.7999718	0.0232674	0.4187238	0.0000118	0.9129647	0.4765686	0.0031581
1079	GO:0036121 double-stranded DNA helicase activity	11	0.0005080	0.4988870	-0.3191473	0.6830146	0.0041709	0.0668625	0.0000875	0.9040200	0.5334497	0.0084584
56	GO:0000793 condensed chromosome	13	0.0002928	0.5272670	-0.4293935	0.5953137	0.0009962	0.0073526	0.0002019	0.9037631	0.5729826	0.0054402
1509	GO:0061749 forked DNA-dependent helicase activity	10	0.0021201	0.5500706	-0.3583686	0.6025120	0.0025952	0.0497469	0.0009694	0.8910816	0.5402628	0.0247893
1740	GO:1990518 single-stranded 3’-5’ DNA helicase activity	11	0.0011133	0.3863524	-0.3684798	0.6710696	0.0265166	0.0343521	0.0001161	0.8575430	0.5371993	0.0159749
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.3697607	0.3619320	-0.6798985	0.0008841	0.0011355	0.0000000	0.8543885	0.5349409	0.0000014
951	GO:0032212 positive regulation of telomere maintenance via telomerase	17	0.0000287	0.3447727	-0.4976101	0.5998505	0.0138651	0.0003823	0.0000185	0.8522351	0.5743247	0.0010771
612	GO:0009378 four-way junction helicase activity	13	0.0003258	0.5814155	-0.2706070	0.5597757	0.0002836	0.0912000	0.0004747	0.8512467	0.4857891	0.0058678
1100	GO:0042026 protein refolding	12	0.0013426	0.3536815	-0.5044225	0.5612679	0.0339073	0.0024824	0.0007611	0.8333993	0.5649676	0.0177244
1525	GO:0070182 DNA polymerase binding	13	0.0012004	0.5580493	-0.3493488	0.4892208	0.0004942	0.0292076	0.0022580	0.8202442	0.5051910	0.0162710
529	GO:0007339 binding of sperm to zona pellucida	11	0.0082636	0.4048152	-0.4887382	0.4557900	0.0200933	0.0050061	0.0088607	0.7813353	0.5312202	0.0672135
135	GO:0001965 G-protein alpha-subunit binding	10	0.0100828	-0.5155095	0.5269404	-0.2350974	0.0047619	0.0039101	0.1980415	0.7737487	0.5394495	0.0760518
1714	GO:1902774 late endosome to lysosome transport	10	0.0000241	-0.5016653	-0.3184166	-0.4826559	0.0060164	0.0812700	0.0082229	0.7655155	0.1007604	0.0009437

top <- head(subset(mmres_rfx$enrichment_result,p.adjustMANOVA<0.05),30)
mycols <- grep("^s\\.",colnames(top))
mycols <- mycols[1:length(mycols)-1]
mx <- as.matrix(top[,mycols])
rownames(mx) <- top$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mx,trace="none",col=colfunc(25),scale="none",
    margins = c(7,25), cexRow=0.8, cexCol=0.8,
    main="Top pathways for Rfx contrasts" )

Now Hf.

l <- list("NTC v DsRed gRNA"=m4,
  "Tfx v NTC"=m5,
  "NTfc v Tfx"=m6)

mm <- mitch_import(l, DEtype="prescored")

## Note: Mean no. genes in input = 14311.6666666667

## Note: no. genes in output = 14014

## Note: estimated proportion of input genes in output = 0.979

head(mm)

##           NTC v DsRed gRNA   Tfx v NTC NTfc v Tfx
## 5_8S_rRNA      -0.01023976  0.05721893  1.0783442
## A4GALT          0.23351781 -0.27502670 -2.3287576
## AAAS            0.21120150 -0.44096365  0.4658765
## AACS            0.57504626  0.38006625  2.9719878
## AADAT          -1.21123488  1.31421790 -1.1764787
## AAGAB           1.19842357 -0.60759817  1.0624617

mmres_hf <- mitch_calc(mm, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mmres_hf$enrichment_result,20) %>%
  kbl(caption = "Top pathways for Hf contrasts") %>%
  kable_paper("hover", full_width = F)

Top pathways for Hf contrasts

	set	setSize	pMANOVA	s.NTC.v.DsRed.gR	s.Tfx.v.NTC	s.NTfc.v.Tfx	p.NTC.v.DsRed.gR	p.Tfx.v.NTC	p.NTfc.v.Tfx	s.dist	SD	p.adjustMANOVA
771	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	0.8779054	-0.5519715	-0.8940282	0.0000000	0.0000000	0.0000000	1.3691884	0.9399734	0.0000000
144	GO:0002181 cytoplasmic translation	27	0.0000000	0.9164727	-0.4698569	-0.8367532	0.0000000	0.0000238	0.0000000	1.3269678	0.9246913	0.0000000
773	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	0.8227236	-0.5079791	-0.7868211	0.0000000	0.0000110	0.0000000	1.2465971	0.8601508	0.0000000
178	GO:0003735 structural constituent of ribosome	118	0.0000000	0.7655465	-0.4359150	-0.7111667	0.0000000	0.0000000	0.0000000	1.1321844	0.7852765	0.0000000
404	GO:0006412 translation	96	0.0000000	0.6883861	-0.4251404	-0.5986343	0.0000000	0.0000000	0.0000000	1.0064705	0.6983866	0.0000000
1209	GO:0045116 protein neddylation	10	0.0000725	0.4485433	-0.2016281	-0.7960154	0.0140516	0.2696391	0.0000130	0.9356738	0.6224877	0.0012660
1666	GO:0140658 ATP-dependent chromatin remodeler activity	11	0.0000609	-0.6199126	0.2469666	0.6390384	0.0003705	0.1561667	0.0002424	0.9239341	0.6442253	0.0011334
341	GO:0005840 ribosome	43	0.0000000	0.6661623	-0.3070097	-0.5223728	0.0000000	0.0004976	0.0000000	0.9005001	0.6332536	0.0000000
311	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	0.6620786	-0.4492438	-0.4088927	0.0000000	0.0000017	0.0000130	0.8985329	0.6302969	0.0000000
343	GO:0005852 eukaryotic translation initiation factor 3 complex	10	0.0001820	0.7851328	-0.3171094	-0.2812625	0.0000171	0.0825288	0.1235808	0.8922446	0.6262883	0.0026523
109	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0000036	0.8113722	-0.2041732	-0.2501428	0.0000011	0.2207870	0.1335780	0.8732600	0.6000361	0.0001030
1221	GO:0045271 respiratory chain complex I	31	0.0000000	0.6003673	-0.3177476	-0.5289487	0.0000000	0.0022053	0.0000003	0.8609246	0.6004020	0.0000000
481	GO:0007076 mitotic chromosome condensation	10	0.0000731	0.0334761	-0.3628820	0.7712796	0.8545824	0.0469401	0.0000240	0.8530394	0.5755833	0.0012660
1721	GO:1990391 DNA repair complex	10	0.0012493	-0.4190803	0.4868466	0.5346758	0.0217567	0.0076824	0.0034150	0.8357787	0.5373766	0.0114081
1070	GO:0036121 double-stranded DNA helicase activity	11	0.0000028	-0.0533717	0.3924808	0.7341089	0.7592622	0.0242120	0.0000248	0.8341497	0.3948882	0.0000831
985	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	0.6230945	-0.2928027	-0.4701292	0.0000000	0.0017944	0.0000005	0.8336676	0.5867212	0.0000000
372	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000009	0.5359546	-0.2071680	-0.5948403	0.0000524	0.1180639	0.0000072	0.8270434	0.5746329	0.0000298
198	GO:0004364 glutathione transferase activity	10	0.0009257	0.5352614	-0.1882034	-0.5965153	0.0033799	0.3028260	0.0010892	0.8232592	0.5731547	0.0093053
1506	GO:0070064 proline-rich region binding	10	0.0006947	-0.5391745	0.1453585	0.6003428	0.0031539	0.4261396	0.0010112	0.8199083	0.5735991	0.0072757
433	GO:0006695 cholesterol biosynthetic process	16	0.0000004	0.6513698	-0.3477193	0.3529076	0.0000064	0.0160521	0.0145391	0.8183735	0.5128575	0.0000164

top <- head(subset(mmres_hf$enrichment_result,p.adjustMANOVA<0.05),30)
mycols <- grep("^s\\.",colnames(top))
mycols <- mycols[1:length(mycols)-1]
mx <- as.matrix(top[,mycols])
rownames(mx) <- top$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mx,trace="none",col=colfunc(25),scale="none",
    margins = c(7,25), cexRow=0.8, cexCol=0.8,
    main="Top pathways for Hf contrasts" )

Now some Rfx vs Hf.

l <- list("Rfx C v gRNA"=m1,
  "Hf C v gRNA"=m4)

mm <- mitch_import(l, DEtype="prescored")

## Note: Mean no. genes in input = 14456

## Note: no. genes in output = 14101

## Note: estimated proportion of input genes in output = 0.975

head(mm)

##           Rfx C v gRNA Hf C v gRNA
## 5_8S_rRNA    1.8503649 -0.01023976
## A4GALT      -1.8915624  0.23351781
## AAAS         2.1906313  0.21120150
## AACS         1.1170591  0.57504626
## AADAC        0.7635603 -0.09350280
## AADAT       -0.2696192 -1.21123488

mmres_vs1 <- mitch_calc(mm, go, priority="effect",cores=8)

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mmres_vs1$enrichment_result,20) %>%
  kbl(caption = "Top pathways for NTC vs gRNA in Rfx and Hf") %>%
  kable_paper("hover", full_width = F)

Top pathways for NTC vs gRNA in Rfx and Hf

	set	setSize	pMANOVA	s.Rfx.C.v.gRNA	s.Hf.C.v.gRNA	p.Rfx.C.v.gRNA	p.Hf.C.v.gRNA	s.dist	SD	p.adjustMANOVA
772	GO:0022625 cytosolic large ribosomal subunit	40	0.0000000	-0.6379063	0.8780030	0.0000000	0.0000000	1.0852712	1.0719097	0.0000000
145	GO:0002181 cytoplasmic translation	27	0.0000000	-0.3703914	0.9165417	0.0008663	0.0000000	0.9885537	0.9099991	0.0000000
774	GO:0022627 cytosolic small ribosomal subunit	25	0.0000000	-0.4967036	0.8227678	0.0000172	0.0000000	0.9610730	0.9330072	0.0000000
180	GO:0003735 structural constituent of ribosome	118	0.0000000	-0.4185603	0.7655773	0.0000000	0.0000000	0.8725259	0.8373118	0.0000000
434	GO:0006695 cholesterol biosynthetic process	16	0.0000000	0.5782570	0.6513578	0.0000621	0.0000064	0.8710041	0.0516901	0.0000022
110	GO:0001732 formation of cytoplasmic translation initiation complex	12	0.0000072	0.0428467	0.8112949	0.7972188	0.0000011	0.8124255	0.5433749	0.0003137
345	GO:0005852 eukaryotic translation initiation factor 3 complex	10	0.0000960	0.0647505	0.7850117	0.7229665	0.0000171	0.7876776	0.5093016	0.0029063
406	GO:0006412 translation	96	0.0000000	-0.3422498	0.6883464	0.0000000	0.0000000	0.7687364	0.7287416	0.0000000
1223	GO:0045271 respiratory chain complex I	31	0.0000000	-0.4369672	0.6004677	0.0000255	0.0000000	0.7426317	0.7335773	0.0000000
374	GO:0006120 mitochondrial electron transport, NADH to ubiquinone	19	0.0000002	-0.4761883	0.5361978	0.0003266	0.0000520	0.7171216	0.7158651	0.0000131
988	GO:0032981 mitochondrial respiratory chain complex I assembly	38	0.0000000	-0.3500900	0.6231320	0.0001888	0.0000000	0.7147423	0.6881719	0.0000000
314	GO:0005762 mitochondrial large ribosomal subunit	38	0.0000000	-0.2319974	0.6620995	0.0133698	0.0000000	0.7015686	0.6322220	0.0000000
437	GO:0006749 glutathione metabolic process	15	0.0000116	-0.4675375	0.5116002	0.0017189	0.0006025	0.6930556	0.6923549	0.0004717
343	GO:0005840 ribosome	43	0.0000000	-0.1888191	0.6661340	0.0322511	0.0000000	0.6923779	0.6045431	0.0000000
1093	GO:0042026 protein refolding	12	0.0003047	0.3522488	0.5876925	0.0346373	0.0004233	0.6851728	0.1664839	0.0074024
694	GO:0016282 eukaryotic 43S preinitiation complex	14	0.0000655	0.1250494	0.6712065	0.4179700	0.0000137	0.6827558	0.3861913	0.0022925
996	GO:0033290 eukaryotic 48S preinitiation complex	13	0.0001751	0.0892413	0.6638274	0.5775218	0.0000340	0.6697991	0.4062937	0.0046392
1011	GO:0034314 Arp2/3 complex-mediated actin nucleation	11	0.0005276	-0.4657333	0.4668688	0.0074842	0.0073399	0.6594498	0.6594493	0.0109864
184	GO:0003774 cytoskeletal motor activity	12	0.0005764	-0.0552559	-0.6434689	0.7403685	0.0001135	0.6458370	0.4159294	0.0115275
266	GO:0005201 extracellular matrix structural constituent	12	0.0005891	-0.0379137	-0.6429247	0.8201394	0.0001150	0.6440417	0.4278074	0.0115275

top <- head(subset(mmres_vs1$enrichment_result,p.adjustMANOVA<0.05),30)
mycols <- grep("^s\\.",colnames(top))
mycols <- mycols[1:length(mycols)-1]
mx <- as.matrix(top[,mycols])
rownames(mx) <- top$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mx,trace="none",col=colfunc(25),scale="none",
    margins = c(7,25), cexRow=0.8, cexCol=0.8,
    main="Top pathways for NTC vs gRNA in Rfx and Hf" )

Session information

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## 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       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] pkgload_1.4.0               GGally_2.2.1               
##  [3] ggplot2_3.5.1               gtools_3.9.5               
##  [5] tibble_3.2.1                echarts4r_0.4.5            
##  [7] UpSetR_1.4.0                beeswarm_0.4.0             
##  [9] kableExtra_1.4.0            eulerr_7.0.2               
## [11] mitch_1.19.2                MASS_7.3-61                
## [13] gplots_3.2.0                DESeq2_1.44.0              
## [15] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [17] MatrixGenerics_1.16.0       matrixStats_1.4.1          
## [19] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [21] IRanges_2.38.1              S4Vectors_0.42.1           
## [23] BiocGenerics_0.50.0         reshape2_1.4.4             
## [25] dplyr_1.1.4                 zoo_1.8-12                 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        viridisLite_0.4.2       farver_2.1.2           
##  [4] bitops_1.0-9            fastmap_1.2.0           promises_1.3.2         
##  [7] digest_0.6.37           mime_0.12               lifecycle_1.0.4        
## [10] polylabelr_0.3.0        magrittr_2.0.3          compiler_4.4.2         
## [13] rlang_1.1.4             sass_0.4.9              tools_4.4.2            
## [16] yaml_2.3.10             knitr_1.49              labeling_0.4.3         
## [19] S4Arrays_1.4.1          htmlwidgets_1.6.4       DelayedArray_0.30.1    
## [22] xml2_1.3.6              plyr_1.8.9              RColorBrewer_1.1-3     
## [25] abind_1.4-8             BiocParallel_1.38.0     KernSmooth_2.23-24     
## [28] withr_3.0.2             purrr_1.0.2             polyclip_1.10-7        
## [31] grid_4.4.2              caTools_1.18.3          xtable_1.8-4           
## [34] colorspace_2.1-1        scales_1.3.0            cli_3.6.3              
## [37] rmarkdown_2.29          crayon_1.5.3            generics_0.1.3         
## [40] rstudioapi_0.17.1       httr_1.4.7              cachem_1.1.0           
## [43] stringr_1.5.1           zlibbioc_1.50.0         parallel_4.4.2         
## [46] XVector_0.44.0          vctrs_0.6.5             Matrix_1.7-1           
## [49] jsonlite_1.8.9          systemfonts_1.1.0       locfit_1.5-9.10        
## [52] jquerylib_0.1.4         tidyr_1.3.1             glue_1.8.0             
## [55] ggstats_0.7.0           codetools_0.2-20        stringi_1.8.4          
## [58] gtable_0.3.6            later_1.4.1             UCSC.utils_1.0.0       
## [61] munsell_0.5.1           pillar_1.10.0           htmltools_0.5.8.1      
## [64] GenomeInfoDbData_1.2.12 R6_2.5.1                evaluate_1.0.1         
## [67] shiny_1.10.0            lattice_0.22-6          httpuv_1.6.15          
## [70] bslib_0.8.0             Rcpp_1.0.13-1           svglite_2.1.3          
## [73] gridExtra_2.3           SparseArray_1.4.8       xfun_0.49              
## [76] pkgconfig_2.0.3