Shoulder instability transcriptome analysis
Mark Ziemann
2024-07-25
Source: https://github.com/markziemann/shoulder-instability-osteroarthritis
Introduction
Here we are analysing the gene expression of six patients in these tissues.
B: Bone
C: Capsule
F: Fat
M: Muscle
S: Synovium
They have a range of QuickDASH scores which I will examine.
suppressPackageStartupMessages({
library("zoo")
library("tidyverse")
library("reshape2")
library("DESeq2")
library("gplots")
library("fgsea")
library("MASS")
library("mitch")
library("eulerr")
library("limma")
library("topconfects")
library("kableExtra")
})Import read counts
tmp <- read.table("3col_oa.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
colnames <- gsub("T0R","T0",colnames(xx))
xx$geneid = NULL
xx <- round(xx)
colnames(xx) <- sapply(strsplit(colnames(xx),"-"),"[[",2)
head(xx)## 41B 41C 41F 41M 41S 42B 42C 42F 42M 42S 43B 43C
## ENSG00000000003.15 TSPAN6 327 231 613 27 346 432 301 579 46 336 279 683
## ENSG00000000005.6 TNMD 9 45 477 1 30 52 2 275 1 308 17 356
## ENSG00000000419.13 DPM1 304 497 485 132 502 349 333 422 126 339 851 815
## ENSG00000000457.14 SCYL3 245 245 280 83 231 202 185 183 93 168 287 281
## ENSG00000000460.17 C1orf112 101 91 79 15 82 86 71 53 9 49 229 105
## ENSG00000000938.13 FGR 2962 828 661 91 736 562 1012 713 95 479 1294 1404
## 43F 43M 43S 44B 44C 44F 44M 44S 45B 45C 45F 45M
## ENSG00000000003.15 TSPAN6 885 69 370 526 386 428 57 301 85 371 368 78
## ENSG00000000005.6 TNMD 772 12 387 3 22 480 2 10 60 157 246 4
## ENSG00000000419.13 DPM1 486 328 617 366 529 218 281 595 395 300 168 493
## ENSG00000000457.14 SCYL3 326 138 278 214 199 177 141 262 261 149 140 157
## ENSG00000000460.17 C1orf112 68 21 89 149 58 32 38 103 179 60 58 38
## ENSG00000000938.13 FGR 421 69 605 4440 775 406 74 803 8028 533 450 72
## 45S 46B 46C 46F 46M 46S
## ENSG00000000003.15 TSPAN6 231 722 342 385 55 363
## ENSG00000000005.6 TNMD 139 40 4 171 13 11
## ENSG00000000419.13 DPM1 414 525 434 273 238 516
## ENSG00000000457.14 SCYL3 184 300 242 147 113 237
## ENSG00000000460.17 C1orf112 70 134 60 47 19 116
## ENSG00000000938.13 FGR 1392 620 968 362 364 1021write.table(xx,"quickdash.tsv",sep="\t")QC analysis
Here I’ll look at a few different quality control measures.
par(mar=c(5,8,3,1))
barplot(colSums(xx),horiz=TRUE,las=1,xlab="num reads")
sums <- colSums(xx)
sums <- sums[order(sums)]
barplot(sums,horiz=TRUE,las=1,xlab="num reads")
abline(v=15000000,col="red")
Sample sheet
URL="https://raw.githubusercontent.com/markziemann/shoulder-instability-osteroarthritis/main/OA_samplesheet.tsv"
ss <- read.table(URL,header=TRUE)
ss$col <- as.numeric(factor(ss$TissueType))
ss$id <- sapply(strsplit(ss$Dataset,"-"),"[[",1)
ss$Dataset <- sapply(strsplit(ss$Dataset,"-"),"[[",2)
URL="https://raw.githubusercontent.com/markziemann/shoulder-instability-osteroarthritis/main/OA_patients.tsv"
pat <- read.table(URL,header=TRUE,sep="\t")
pat$id <- sapply(strsplit(pat$Record_number," "),"[[",2)MDS plot
Multidimensional scaling plot to show the variation between samples, very similar to PCA.
colvec <- c("orange2","pink","purple1","yellow","lightblue")
cols <- rep(colvec,6)
plot(cmdscale(dist(t(xx))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n",pch=19, cex=4, col=cols)
text(cmdscale(dist(t(xx))), labels=ss$ParticipantID )
legend("topright", inset=0.1, title="tissue type", c("Bone","Capsule","Fat","Muscle","Synovium"),
pch=19,col=colvec, cex=1.5)
Correlation heatmap
heatmap.2(cor(xx),trace="n",main="Pearson correlation heatmap")
There were 4 samples with fewer than 15M reads and this might impact data quality: 4002-42S, 4006-46B, 4005-42F and 4006-46C.
MDS plot shows samples clustered by tissue type as expected. For example M samples to the right, F samples in the centre, C samples in the lower right and S samples between C and F. However the B samples were somewhat not clustered, which is something to look into.
The result of the MDS is confirmed in the correlation heatmap.
Separate data by tissue
Also a good point to filter out any genes with low expression (average < 10 counts).
Below I show the dimensions of each dataset (separated by tissue type), in terms of number of genes detected and samples.
ssb <- subset(ss,TissueType=="B")
xxb <- xx[,which(colnames(xx) %in% ssb$Dataset)]
xxb <- xxb[which(rowMeans(xxb)>10),]
ssc <- subset(ss,TissueType=="C")
xxc <- xx[,which(colnames(xx) %in% ssc$Dataset)]
xxc <- xxc[which(rowMeans(xxc)>10),]
ssf <- subset(ss,TissueType=="F")
xxf <- xx[,which(colnames(xx) %in% ssf$Dataset)]
xxf <- xxf[which(rowMeans(xxf)>10),]
ssm <- subset(ss,TissueType=="M")
xxm <- xx[,which(colnames(xx) %in% ssm$Dataset)]
xxm <- xxm[which(rowMeans(xxm)>10),]
sss <- subset(ss,TissueType=="S")
xxs <- xx[,which(colnames(xx) %in% sss$Dataset)]
xxs <- xxs[which(rowMeans(xxs)>10),]
lapply(list(xx,xxb,xxc,xxf,xxm,xxs),dim)## [[1]]
## [1] 60651 30
##
## [[2]]
## [1] 20241 6
##
## [[3]]
## [1] 18726 6
##
## [[4]]
## [1] 19072 6
##
## [[5]]
## [1] 15920 6
##
## [[6]]
## [1] 19361 6Analysis of gene expression correlates with basline quickdash scores
ssb$quickdash <- scale(pat[match(ssb$id,pat$id),"QuickDASH_baseline"])
ssb %>% kbl(caption="samplesheet for bone") %>% kable_paper("hover", full_width = F)| Dataset | ParticipantID | TissueType | Age | Sex | col | id | quickdash | |
|---|---|---|---|---|---|---|---|---|
| 1 | 41B | 41 | B | 23 | M | 1 | 4001 | -0.51634096 |
| 6 | 42B | 42 | B | 35 | M | 1 | 4002 | 1.15023003 |
| 11 | 43B | 42 | B | 35 | M | 1 | 4003 | 0.07923613 |
| 16 | 44B | 44 | B | 27 | M | 1 | 4004 | 1.15023003 |
| 21 | 45B | 45 | B | 27 | M | 1 | 4005 | -1.34701427 |
| 26 | 46B | 46 | B | 20 | M | 1 | 4006 | -0.51634096 |
dds <- DESeqDataSetFromMatrix(countData = xxb , colData = ssb, design = ~ quickdash )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- 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 correlates with QuickDASH in bone") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | 41B | 42B | 43B | 44B | 45B | 46B | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000211638.2 IGLV8-61 | 440.91837 | -2.931254 | 0.5159649 | -5.681111 | 0.0000000 | 0.0002709 | 8.930659 | 4.205568 | 8.597901 | 3.501237 | 10.671973 | 7.232449 |
| ENSG00000274419.6 TBC1D3D | 46.71931 | -3.506396 | 0.7067664 | -4.961181 | 0.0000007 | 0.0048521 | 5.716401 | 1.861318 | 5.328218 | 1.861318 | 7.103516 | 6.290270 |
| ENSG00000145423.5 SFRP2 | 460.55116 | -1.786447 | 0.3604527 | -4.956120 | 0.0000007 | 0.0048521 | 9.225523 | 5.445464 | 6.876516 | 6.747061 | 10.415843 | 9.162634 |
| ENSG00000196436.8 NPIPB15 | 83.47206 | -2.540972 | 0.5571418 | -4.560727 | 0.0000051 | 0.0257955 | 7.177856 | 3.405667 | 3.060644 | 3.501237 | 7.693135 | 7.305542 |
| ENSG00000271130.1 IGHV3OR16-8 | 27.99582 | -3.564167 | 0.8218232 | -4.336903 | 0.0000145 | 0.0584983 | 5.040078 | 1.861318 | 3.840140 | 1.861318 | 7.123073 | 2.839618 |
| ENSG00000095752.7 IL11 | 298.67564 | 3.846244 | 0.9014247 | 4.266851 | 0.0000198 | 0.0668802 | 3.780852 | 3.079033 | 2.721216 | 10.789836 | 2.586595 | 4.535245 |
| ENSG00000073756.13 PTGS2 | 141.73199 | 1.285224 | 0.3072010 | 4.183659 | 0.0000287 | 0.0829461 | 5.815436 | 8.399398 | 6.955978 | 8.123029 | 5.620540 | 5.551745 |
| ENSG00000270149.5 AL591806.2 | 105.30108 | -1.121805 | 0.2727571 | -4.112835 | 0.0000391 | 0.0988847 | 7.021321 | 5.213977 | 6.373073 | 5.544574 | 8.085718 | 6.729005 |
| ENSG00000281383.1 FP671120.7 | 60.01779 | 1.433140 | 0.3545168 | 4.042517 | 0.0000529 | 0.1189281 | 4.910451 | 6.903751 | 6.396072 | 7.003666 | 4.151981 | 4.419797 |
| ENSG00000214940.8 NPIPA8 | 88.24011 | -2.685878 | 0.6790047 | -3.955611 | 0.0000763 | 0.1545185 | 4.767596 | 3.079033 | 3.514704 | 3.342451 | 8.967258 | 3.220658 |
| ENSG00000130720.13 FIBCD1 | 12.86486 | 3.555247 | 0.9047775 | 3.929416 | 0.0000852 | 0.1566884 | 2.583587 | 3.767650 | 1.861318 | 6.248497 | 1.861318 | 1.861318 |
| ENSG00000119508.18 NR4A3 | 210.56959 | 2.481732 | 0.6488879 | 3.824593 | 0.0001310 | 0.2106384 | 4.663511 | 10.190786 | 4.810847 | 4.675312 | 3.982721 | 5.350149 |
| ENSG00000153993.13 SEMA3D | 263.02297 | 1.762272 | 0.4617345 | 3.816634 | 0.0001353 | 0.2106384 | 7.681316 | 8.579345 | 7.144141 | 9.615365 | 3.787985 | 6.539845 |
| ENSG00000229183.9 PGA4 | 13.52739 | -3.644256 | 0.9753663 | -3.736295 | 0.0001868 | 0.2700033 | 3.780852 | 1.861318 | 3.060644 | 1.861318 | 4.070076 | 6.083706 |
| ENSG00000039537.14 C6 | 98.93901 | -1.159004 | 0.3170880 | -3.655151 | 0.0002570 | 0.3468372 | 6.489877 | 5.475657 | 6.147608 | 4.947782 | 7.775736 | 7.497997 |
| ENSG00000134184.13 GSTM1 | 42.95966 | -7.719138 | 2.1431726 | -3.601734 | 0.0003161 | 0.3899222 | 1.861318 | 1.861318 | 1.861318 | 1.861318 | 1.861318 | 8.049712 |
| ENSG00000223931.2 AC142381.1 | 57.50480 | -1.931489 | 0.5376410 | -3.592526 | 0.0003275 | 0.3899222 | 5.268256 | 3.538885 | 3.514704 | 3.984908 | 8.153817 | 4.478738 |
| ENSG00000136634.7 IL10 | 70.99991 | 2.008719 | 0.5632486 | 3.566310 | 0.0003620 | 0.4071179 | 4.663511 | 8.387616 | 5.465747 | 4.335672 | 3.270504 | 5.001029 |
| ENSG00000137077.8 CCL21 | 15.70882 | -1.701534 | 0.4914766 | -3.462085 | 0.0005360 | 0.5710177 | 3.975064 | 2.864993 | 4.098614 | 3.157099 | 6.005978 | 4.225850 |
| ENSG00000153266.13 FEZF2 | 26.86972 | 2.060832 | 0.5992079 | 3.439260 | 0.0005833 | 0.5903348 | 1.861318 | 6.677080 | 3.840140 | 5.810519 | 3.270504 | 3.499068 |
dgeb <- dge
write.table(dgeb,file="deseq2_bone.tsv",quote=FALSE,sep='\t')
ssc$quickdash <- scale(pat[match(ssc$id,pat$id),"QuickDASH_baseline"])
ssb %>% kbl(caption="samplesheet for capsule") %>% kable_paper("hover", full_width = F)| Dataset | ParticipantID | TissueType | Age | Sex | col | id | quickdash | |
|---|---|---|---|---|---|---|---|---|
| 1 | 41B | 41 | B | 23 | M | 1 | 4001 | -0.51634096 |
| 6 | 42B | 42 | B | 35 | M | 1 | 4002 | 1.15023003 |
| 11 | 43B | 42 | B | 35 | M | 1 | 4003 | 0.07923613 |
| 16 | 44B | 44 | B | 27 | M | 1 | 4004 | 1.15023003 |
| 21 | 45B | 45 | B | 27 | M | 1 | 4005 | -1.34701427 |
| 26 | 46B | 46 | B | 20 | M | 1 | 4006 | -0.51634096 |
dds <- DESeqDataSetFromMatrix(countData = xxc , colData = ssc, design = ~ quickdash )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- 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 correlates with QuickDASH in capsule") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | 41C | 42C | 43C | 44C | 45C | 46C | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000248144.6 ADH1C | 49.68522 | -2.0385840 | 0.3247349 | -6.277686 | 0.00e+00 | 0.0000034 | 6.248232 | 4.879783 | 5.703559 | 5.255560 | 7.706001 | 6.636002 |
| ENSG00000274419.6 TBC1D3D | 73.36200 | -4.1811328 | 0.6672278 | -6.266425 | 0.00e+00 | 0.0000034 | 6.658670 | 4.403795 | 5.678897 | 4.403795 | 7.624139 | 7.922841 |
| ENSG00000039537.14 C6 | 299.22535 | -2.4279663 | 0.4208040 | -5.769827 | 0.00e+00 | 0.0000486 | 10.016795 | 5.449903 | 7.719263 | 5.867846 | 9.028609 | 7.267074 |
| ENSG00000142178.9 SIK1 | 110.70007 | 2.3193982 | 0.4383154 | 5.291619 | 1.00e-07 | 0.0004610 | 5.965574 | 7.119941 | 5.355726 | 9.086798 | 5.067327 | 5.773043 |
| ENSG00000125730.17 C3 | 2208.32354 | -1.8826918 | 0.3562169 | -5.285240 | 1.00e-07 | 0.0004610 | 12.832108 | 8.257105 | 10.190717 | 8.351984 | 11.663371 | 10.183791 |
| ENSG00000186431.19 FCAR | 30.38540 | 2.5646086 | 0.4913965 | 5.219021 | 2.00e-07 | 0.0004781 | 5.078262 | 7.325994 | 5.242080 | 6.292669 | 4.403795 | 5.254083 |
| ENSG00000166482.12 MFAP4 | 3157.73936 | -0.9271698 | 0.1777349 | -5.216589 | 2.00e-07 | 0.0004781 | 12.438460 | 10.259131 | 11.298717 | 10.421206 | 12.450246 | 11.519873 |
| ENSG00000176945.17 MUC20 | 344.82717 | -1.1779858 | 0.2346290 | -5.020632 | 5.00e-07 | 0.0011821 | 9.574743 | 7.107541 | 8.080757 | 7.126881 | 9.084261 | 8.841379 |
| ENSG00000211669.3 IGLV3-10 | 19.79640 | -3.9906631 | 0.8137022 | -4.904329 | 9.00e-07 | 0.0019127 | 5.008137 | 4.403795 | 4.403795 | 4.403795 | 7.260605 | 4.729155 |
| ENSG00000159263.16 SIM2 | 117.11826 | -1.1230524 | 0.2331657 | -4.816542 | 1.50e-06 | 0.0026821 | 7.279374 | 5.919030 | 7.111892 | 6.451796 | 8.469127 | 7.102895 |
| ENSG00000124107.5 SLPI | 135.25645 | -1.2871914 | 0.2741228 | -4.695674 | 2.70e-06 | 0.0044357 | 8.295381 | 5.886992 | 7.401278 | 6.214895 | 7.919685 | 7.529725 |
| ENSG00000188257.12 PLA2G2A | 32861.23105 | -1.5584845 | 0.3479723 | -4.478760 | 7.50e-06 | 0.0113600 | 15.972762 | 11.078605 | 15.408404 | 12.933800 | 15.194721 | 15.362567 |
| ENSG00000197766.8 CFD | 10476.63157 | -1.2618786 | 0.2833101 | -4.454055 | 8.40e-06 | 0.0113600 | 14.658820 | 10.791169 | 12.998718 | 11.843394 | 13.934758 | 12.957715 |
| ENSG00000163735.7 CXCL5 | 57.91874 | 2.7383282 | 0.6156688 | 4.447730 | 8.70e-06 | 0.0113600 | 4.707618 | 8.423334 | 5.281803 | 5.754112 | 4.737747 | 5.590062 |
| ENSG00000128383.13 APOBEC3A | 68.97410 | 1.2533513 | 0.2827128 | 4.433303 | 9.30e-06 | 0.0113600 | 5.915336 | 7.408607 | 6.106610 | 7.773314 | 5.873228 | 5.838137 |
| ENSG00000136689.19 IL1RN | 78.20007 | 1.3717917 | 0.3116023 | 4.402381 | 1.07e-05 | 0.0122876 | 5.989888 | 8.106177 | 6.824687 | 6.787031 | 5.439799 | 6.402780 |
| ENSG00000145824.13 CXCL14 | 1562.57168 | -1.2371931 | 0.2835391 | -4.363396 | 1.28e-05 | 0.0138319 | 11.450332 | 8.940588 | 11.098540 | 8.177586 | 10.959583 | 10.863248 |
| ENSG00000162747.12 FCGR3B | 332.89923 | 1.4785624 | 0.3458081 | 4.275673 | 1.91e-05 | 0.0189012 | 7.158334 | 9.732960 | 6.338019 | 9.844083 | 6.942508 | 6.875852 |
| ENSG00000205221.12 VIT | 242.40661 | -1.5450981 | 0.3618600 | -4.269878 | 1.96e-05 | 0.0189012 | 9.380209 | 5.545030 | 7.375552 | 6.862148 | 8.900090 | 7.750909 |
| ENSG00000122121.12 XPNPEP2 | 92.06925 | -1.4711424 | 0.3566169 | -4.125274 | 3.70e-05 | 0.0339967 | 8.059610 | 5.343401 | 6.918292 | 5.894526 | 7.590016 | 6.568593 |
dgec <- dge
write.table(dgec,file="deseq2_capsule.tsv",quote=FALSE,sep='\t')
ssf$quickdash <- scale(pat[match(ssf$id,pat$id),"QuickDASH_baseline"])
ssb %>% kbl(caption="samplesheet for fat") %>% kable_paper("hover", full_width = F)| Dataset | ParticipantID | TissueType | Age | Sex | col | id | quickdash | |
|---|---|---|---|---|---|---|---|---|
| 1 | 41B | 41 | B | 23 | M | 1 | 4001 | -0.51634096 |
| 6 | 42B | 42 | B | 35 | M | 1 | 4002 | 1.15023003 |
| 11 | 43B | 42 | B | 35 | M | 1 | 4003 | 0.07923613 |
| 16 | 44B | 44 | B | 27 | M | 1 | 4004 | 1.15023003 |
| 21 | 45B | 45 | B | 27 | M | 1 | 4005 | -1.34701427 |
| 26 | 46B | 46 | B | 20 | M | 1 | 4006 | -0.51634096 |
dds <- DESeqDataSetFromMatrix(countData = xxf , colData = ssf, design = ~ quickdash )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- 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 correlates with QuickDASH in fat") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | 41F | 42F | 43F | 44F | 45F | 46F | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000073756.13 PTGS2 | 119.90779 | 3.7022183 | 0.6818779 | 5.429445 | 1.00e-07 | 0.0010781 | 4.734032 | 9.468892 | 4.825193 | 4.775833 | 3.498061 | 4.413838 |
| ENSG00000278299.5 TBC1D3C | 21.01376 | -3.8365212 | 0.7254307 | -5.288612 | 1.00e-07 | 0.0011753 | 5.291506 | 3.498061 | 3.498061 | 3.498061 | 6.857765 | 4.913062 |
| ENSG00000196436.8 NPIPB15 | 52.46229 | -2.5979241 | 0.4991609 | -5.204582 | 2.00e-07 | 0.0012361 | 6.745880 | 4.262497 | 4.251850 | 4.152412 | 7.328100 | 6.680265 |
| ENSG00000273294.1 C1QTNF3-AMACR | 25.76986 | -3.9396329 | 0.7653204 | -5.147691 | 3.00e-07 | 0.0012574 | 5.349143 | 3.498061 | 3.498061 | 3.498061 | 6.857765 | 5.889474 |
| ENSG00000019169.11 MARCO | 1154.06922 | -1.0710154 | 0.2205026 | -4.857155 | 1.20e-06 | 0.0045424 | 10.428500 | 8.474129 | 9.269474 | 9.092977 | 11.571676 | 10.126887 |
| ENSG00000205502.4 C2CD4B | 204.79684 | 2.4598164 | 0.5102746 | 4.820574 | 1.40e-06 | 0.0045501 | 5.134134 | 10.143152 | 5.379069 | 6.061604 | 4.891409 | 5.473738 |
| ENSG00000214940.8 NPIPA8 | 80.18742 | -2.6254709 | 0.5669613 | -4.630776 | 3.60e-06 | 0.0099255 | 5.064806 | 4.477605 | 4.338507 | 4.152412 | 8.850932 | 5.094199 |
| ENSG00000258430.1 AL583722.1 | 21.01845 | -3.2268797 | 0.7033708 | -4.587736 | 4.50e-06 | 0.0104211 | 5.134134 | 3.498061 | 4.613647 | 3.498061 | 6.261418 | 6.082966 |
| ENSG00000164694.17 FNDC1 | 1731.21142 | -1.3540428 | 0.2964017 | -4.568269 | 4.90e-06 | 0.0104211 | 10.629144 | 8.177412 | 9.233994 | 9.428736 | 12.584622 | 10.246329 |
| ENSG00000164825.4 DEFB1 | 158.38918 | -1.3044346 | 0.2959097 | -4.408218 | 1.04e-05 | 0.0181199 | 7.372966 | 5.552629 | 7.900770 | 5.893925 | 8.427900 | 7.802637 |
| ENSG00000215914.4 MMP23A | 17.82400 | -2.9069569 | 0.6595285 | -4.407629 | 1.05e-05 | 0.0181199 | 4.991015 | 3.498061 | 4.725032 | 3.498061 | 6.346294 | 5.473738 |
| ENSG00000169429.11 CXCL8 | 93.23866 | 2.4216015 | 0.5614329 | 4.313252 | 1.61e-05 | 0.0255677 | 4.826800 | 9.026023 | 5.287173 | 5.348505 | 4.593292 | 4.150954 |
| ENSG00000153266.13 FEZF2 | 25.54163 | 2.5841825 | 0.6045138 | 4.274811 | 1.91e-05 | 0.0280651 | 4.225928 | 7.068360 | 4.416223 | 5.667402 | 4.140369 | 3.498061 |
| ENSG00000261654.1 AL360270.2 | 356.86548 | 2.1074011 | 0.4975544 | 4.235519 | 2.28e-05 | 0.0310634 | 6.879139 | 10.852957 | 7.096235 | 5.833049 | 5.561161 | 5.889474 |
| ENSG00000274419.6 TBC1D3D | 48.89398 | -2.3660608 | 0.5609489 | -4.217961 | 2.47e-05 | 0.0313443 | 5.028511 | 4.477605 | 5.644417 | 3.498061 | 6.828880 | 7.440265 |
| ENSG00000106976.21 DNM1 | 1568.41911 | -0.9708034 | 0.2322853 | -4.179358 | 2.92e-05 | 0.0338734 | 11.424041 | 9.081865 | 10.265546 | 9.475019 | 11.629322 | 10.161399 |
| ENSG00000277027.1 RMRP | 17.32159 | -2.8477992 | 0.6825976 | -4.172003 | 3.02e-05 | 0.0338734 | 5.099987 | 3.498061 | 4.871869 | 3.498061 | 5.942595 | 5.828673 |
| ENSG00000100473.18 COCH | 128.93323 | -1.3349220 | 0.3223252 | -4.141538 | 3.45e-05 | 0.0353936 | 7.835377 | 6.094936 | 6.324899 | 5.348505 | 8.640548 | 6.382682 |
| ENSG00000215246.5 AC116351.1 | 47.55904 | 2.3707474 | 0.5731249 | 4.136528 | 3.53e-05 | 0.0353936 | 5.134134 | 8.023617 | 5.116596 | 4.775833 | 3.954101 | 4.413838 |
| ENSG00000104921.15 FCER2 | 130.89946 | -1.2035532 | 0.2934860 | -4.100888 | 4.12e-05 | 0.0375994 | 8.151345 | 5.402904 | 6.809562 | 6.187101 | 8.053827 | 7.085475 |
dgef <- dge
write.table(dgef,file="deseq2_fat.tsv",quote=FALSE,sep='\t')
ssm$quickdash <- scale(pat[match(ssm$id,pat$id),"QuickDASH_baseline"])
ssb %>% kbl(caption="samplesheet for muscle") %>% kable_paper("hover", full_width = F)| Dataset | ParticipantID | TissueType | Age | Sex | col | id | quickdash | |
|---|---|---|---|---|---|---|---|---|
| 1 | 41B | 41 | B | 23 | M | 1 | 4001 | -0.51634096 |
| 6 | 42B | 42 | B | 35 | M | 1 | 4002 | 1.15023003 |
| 11 | 43B | 42 | B | 35 | M | 1 | 4003 | 0.07923613 |
| 16 | 44B | 44 | B | 27 | M | 1 | 4004 | 1.15023003 |
| 21 | 45B | 45 | B | 27 | M | 1 | 4005 | -1.34701427 |
| 26 | 46B | 46 | B | 20 | M | 1 | 4006 | -0.51634096 |
dds <- DESeqDataSetFromMatrix(countData = xxm , colData = ssm, design = ~ quickdash )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- 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 correlates with QuickDASH in muscle") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | 41M | 42M | 43M | 44M | 45M | 46M | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000274419.6 TBC1D3D | 19.41891 | -3.6142012 | 0.7764588 | -4.654724 | 0.0000032 | 0.0516470 | 5.482788 | 4.094096 | 4.710388 | 4.094096 | 5.934661 | 6.578833 |
| ENSG00000160888.7 IER2 | 518.83277 | 1.4210313 | 0.3193194 | 4.450188 | 0.0000086 | 0.0612357 | 7.714334 | 11.051849 | 7.892072 | 8.332358 | 7.098796 | 8.135541 |
| ENSG00000081041.9 CXCL2 | 76.78205 | 1.7761791 | 0.4104752 | 4.327128 | 0.0000151 | 0.0612357 | 4.808246 | 8.520217 | 5.875837 | 6.418142 | 5.287051 | 5.776348 |
| ENSG00000171617.15 ENC1 | 190.11873 | 0.8976600 | 0.2076431 | 4.323090 | 0.0000154 | 0.0612357 | 6.962449 | 8.411764 | 7.872623 | 8.627597 | 6.851619 | 6.965453 |
| ENSG00000196436.8 NPIPB15 | 16.76047 | -2.9244111 | 0.7262633 | -4.026654 | 0.0000566 | 0.1610387 | 5.398903 | 4.460127 | 4.451708 | 4.094096 | 6.337441 | 5.969007 |
| ENSG00000134184.13 GSTM1 | 119.57051 | -8.5934024 | 2.1429395 | -4.010100 | 0.0000607 | 0.1610387 | 4.094096 | 4.094096 | 4.094096 | 4.094096 | 4.094096 | 9.553037 |
| ENSG00000154734.16 ADAMTS1 | 666.03357 | 0.9346697 | 0.2388152 | 3.913778 | 0.0000909 | 0.2066490 | 8.827198 | 10.886483 | 9.145379 | 9.225506 | 8.121563 | 8.668636 |
| ENSG00000204460.3 LINC01854 | 381.98127 | 1.0571224 | 0.2877602 | 3.673623 | 0.0002391 | 0.4709364 | 8.706825 | 8.866629 | 8.038808 | 10.070861 | 7.117544 | 6.996709 |
| ENSG00000224687.2 RASAL2-AS1 | 38.21144 | 2.0203818 | 0.5541188 | 3.646117 | 0.0002662 | 0.4709364 | 4.808246 | 4.904082 | 4.885817 | 7.868738 | 4.819832 | 5.117292 |
| ENSG00000081181.8 ARG2 | 139.55961 | -1.4436456 | 0.3990482 | -3.617723 | 0.0002972 | 0.4731530 | 5.941149 | 5.726977 | 5.957945 | 6.267292 | 9.412878 | 6.399891 |
| ENSG00000168621.15 GDNF | 610.01979 | 0.9797612 | 0.2741525 | 3.573782 | 0.0003519 | 0.5092399 | 8.931970 | 10.869537 | 8.249185 | 9.371930 | 8.121563 | 8.000207 |
| ENSG00000136244.12 IL6 | 13.43905 | 2.4760095 | 0.7132823 | 3.471290 | 0.0005180 | 0.6871673 | 4.453894 | 6.656734 | 4.803969 | 4.698627 | 4.421385 | 4.559315 |
| ENSG00000122877.16 EGR2 | 27.52304 | 1.7208326 | 0.5036746 | 3.416556 | 0.0006342 | 0.7766346 | 4.601627 | 7.254736 | 5.251880 | 5.278482 | 4.744556 | 5.410203 |
| ENSG00000162772.17 ATF3 | 180.42876 | 1.0884805 | 0.3249612 | 3.349570 | 0.0008094 | 0.9203697 | 5.734328 | 8.688426 | 6.844778 | 8.616434 | 6.432825 | 7.739545 |
| ENSG00000125740.14 FOSB | 727.74879 | 3.5458419 | 1.0869717 | 3.262129 | 0.0011058 | 0.9999305 | 4.453894 | 12.028100 | 7.467668 | 4.870848 | 4.819832 | 6.646910 |
| ENSG00000170345.10 FOS | 5286.85345 | 3.0332649 | 0.9371373 | 3.236735 | 0.0012091 | 0.9999305 | 6.094375 | 14.861025 | 8.847656 | 8.446097 | 6.432825 | 10.221335 |
| ENSG00000177606.8 JUN | 7279.99453 | 0.9818296 | 0.3081936 | 3.185756 | 0.0014438 | 0.9999305 | 11.723421 | 14.584456 | 11.885568 | 11.831670 | 11.273073 | 12.558502 |
| ENSG00000267247.1 AP005264.2 | 43.98595 | 1.1211783 | 0.3547236 | 3.160710 | 0.0015739 | 0.9999305 | 5.797137 | 7.470043 | 5.931180 | 5.926055 | 5.364608 | 5.410203 |
| ENSG00000159399.10 HK2 | 766.44623 | -0.7636722 | 0.2417055 | -3.159515 | 0.0015803 | 0.9999305 | 9.865551 | 7.884404 | 9.131868 | 9.139829 | 10.107042 | 10.454723 |
| ENSG00000128965.13 CHAC1 | 31.98198 | 1.4449630 | 0.4574674 | 3.158614 | 0.0015852 | 0.9999305 | 5.207028 | 7.296406 | 5.251880 | 5.409718 | 4.744556 | 5.776348 |
dgem <- dge
write.table(dgem,file="deseq2_muscle.tsv",quote=FALSE,sep='\t')
sss$quickdash <- scale(pat[match(sss$id,pat$id),"QuickDASH_baseline"])
ssb %>% kbl(caption="samplesheet for synovium") %>% kable_paper("hover", full_width = F)| Dataset | ParticipantID | TissueType | Age | Sex | col | id | quickdash | |
|---|---|---|---|---|---|---|---|---|
| 1 | 41B | 41 | B | 23 | M | 1 | 4001 | -0.51634096 |
| 6 | 42B | 42 | B | 35 | M | 1 | 4002 | 1.15023003 |
| 11 | 43B | 42 | B | 35 | M | 1 | 4003 | 0.07923613 |
| 16 | 44B | 44 | B | 27 | M | 1 | 4004 | 1.15023003 |
| 21 | 45B | 45 | B | 27 | M | 1 | 4005 | -1.34701427 |
| 26 | 46B | 46 | B | 20 | M | 1 | 4006 | -0.51634096 |
dds <- DESeqDataSetFromMatrix(countData = xxs , colData = sss, design = ~ quickdash )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz<- 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 correlates with QuickDASH in synovium") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | 41S | 42S | 43S | 44S | 45S | 46S | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000119508.18 NR4A3 | 481.25337 | 2.5912836 | 0.4224586 | 6.133817 | 0.0e+00 | 0.0000089 | 6.342950 | 11.312129 | 7.252746 | 7.529529 | 5.661256 | 6.794226 |
| ENSG00000211893.4 IGHG2 | 610.05143 | -2.7334236 | 0.4464040 | -6.123206 | 0.0e+00 | 0.0000089 | 6.594382 | 5.936293 | 5.962233 | 5.715193 | 11.567969 | 9.184508 |
| ENSG00000198435.4 NRARP | 381.16544 | 1.1179232 | 0.1899825 | 5.884347 | 0.0e+00 | 0.0000258 | 8.244739 | 9.925151 | 8.107532 | 9.270060 | 7.120048 | 7.842529 |
| ENSG00000137962.13 ARHGAP29 | 654.33270 | 0.7928283 | 0.1461255 | 5.425668 | 1.0e-07 | 0.0002795 | 9.006135 | 10.173490 | 9.329650 | 10.002710 | 8.044563 | 9.176952 |
| ENSG00000069122.19 ADGRF5 | 958.81660 | 0.6851472 | 0.1322419 | 5.181016 | 2.0e-07 | 0.0008545 | 9.642179 | 10.460410 | 9.784922 | 10.690402 | 8.798755 | 9.682777 |
| ENSG00000160200.18 CBS | 153.06842 | -1.2082613 | 0.2348533 | -5.144749 | 3.0e-07 | 0.0008644 | 7.266678 | 6.662419 | 7.209423 | 6.090989 | 8.967252 | 7.535559 |
| ENSG00000172201.12 ID4 | 1039.66904 | 0.8189104 | 0.1620821 | 5.052442 | 4.0e-07 | 0.0011039 | 9.646803 | 11.144889 | 9.912144 | 10.429460 | 8.931841 | 9.224151 |
| ENSG00000091879.14 ANGPT2 | 263.11042 | 1.0103110 | 0.2003035 | 5.043900 | 5.0e-07 | 0.0011039 | 7.618551 | 8.864783 | 7.501996 | 9.301480 | 6.891241 | 7.861650 |
| ENSG00000128917.8 DLL4 | 474.37653 | 0.8229100 | 0.1672282 | 4.920881 | 9.0e-07 | 0.0018534 | 8.392112 | 9.962565 | 9.172649 | 9.315911 | 7.897936 | 8.290204 |
| ENSG00000255422.3 AP002954.1 | 34.26502 | -1.4750875 | 0.3026374 | -4.874108 | 1.1e-06 | 0.0021162 | 6.048514 | 5.284400 | 5.839404 | 5.209793 | 6.984841 | 6.504942 |
| ENSG00000102755.12 FLT1 | 765.33467 | 0.6504951 | 0.1366575 | 4.760038 | 1.9e-06 | 0.0033008 | 9.274835 | 10.300166 | 9.337561 | 10.250167 | 8.628099 | 9.466493 |
| ENSG00000153721.19 CNKSR3 | 270.61973 | 0.9462839 | 0.1992663 | 4.748841 | 2.0e-06 | 0.0033008 | 7.282919 | 9.416852 | 8.112175 | 8.592940 | 7.154720 | 7.829637 |
| ENSG00000211660.3 IGLV2-23 | 127.58748 | -4.7808741 | 1.0122016 | -4.723243 | 2.3e-06 | 0.0034569 | 4.842938 | 4.439336 | 4.849674 | 4.439336 | 9.585722 | 6.197751 |
| ENSG00000142178.9 SIK1 | 83.51102 | 2.5309590 | 0.5568599 | 4.545056 | 5.5e-06 | 0.0075952 | 5.961863 | 6.699988 | 6.249406 | 8.681734 | 4.897052 | 4.439336 |
| ENSG00000270550.1 IGHV3-30 | 64.50448 | -3.2358005 | 0.7200406 | -4.493914 | 7.0e-06 | 0.0088325 | 5.133949 | 4.439336 | 4.729976 | 4.947109 | 8.607821 | 6.059306 |
| ENSG00000239922.2 AC092958.1 | 52.08835 | -1.0543928 | 0.2353789 | -4.479555 | 7.5e-06 | 0.0088325 | 6.746456 | 5.663994 | 6.231031 | 5.801409 | 7.346255 | 6.398354 |
| ENSG00000224799.1 AL139397.1 | 22.12071 | 3.9014647 | 0.8728998 | 4.469545 | 7.8e-06 | 0.0088325 | 4.439336 | 7.360643 | 5.200730 | 4.947109 | 4.439336 | 4.439336 |
| ENSG00000196436.8 NPIPB15 | 94.67156 | -2.2311487 | 0.5003033 | -4.459592 | 8.2e-06 | 0.0088325 | 7.322710 | 5.570767 | 4.849674 | 4.947109 | 8.510659 | 7.204041 |
| ENSG00000164035.10 EMCN | 921.78699 | 0.7818082 | 0.1766971 | 4.424567 | 9.7e-06 | 0.0092012 | 9.262829 | 10.632730 | 9.374555 | 10.648705 | 8.599629 | 9.858255 |
| ENSG00000120279.6 MYCT1 | 300.89766 | 0.8354745 | 0.1890875 | 4.418453 | 9.9e-06 | 0.0092012 | 8.014715 | 9.307625 | 8.354564 | 8.697591 | 6.997701 | 8.193281 |
dges <- dge
write.table(dges,file="deseq2_synovium.tsv",quote=FALSE,sep='\t')Here let’s look at some plots.
MA plot shows the avrage level and fold change of all detected genes. Volcano plot shows the fold change and the significance, as measured by -log(p-value). Significant genes are shown as red points.
There are heatmaps of the top ranked genes by p-value. Above the gene expression values there is a bar in yellow/orange colours. Yellow shows relatively low QuickDASH score and orange shows high score.
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(log2(de$baseMean),de$log2FoldChange,
xlab="log2 basemean", ylab="log2 foldchange",
pch=19, cex=0.5, col="dark gray",
main=contrast_name, cex.main=1)
points(log2(sig$baseMean),sig$log2FoldChange,
pch=19, cex=0.5, col="red")
mtext(SUBHEADER,cex = 1)
}
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="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(10,25), cexRow=1.5,
main=paste("Top ranked",n,"genes in",name) , ColSideColors = colcols )
}
maplot(dgeb,"bone")
make_volcano(dgeb,"bone")
make_heatmap(dgeb,"bone",ssb,xxb,n=20)
maplot(dgec,"capsule")
make_volcano(dgec,"capsule")
make_heatmap(dgec,"capsule",ssc,xxc,n=20)
maplot(dgef,"fat")
make_volcano(dgef,"fat")
make_heatmap(dgef,"fat",ssf,xxf,n=20)
maplot(dgem,"muscle")
make_volcano(dgem,"muscle")
make_heatmap(dgem,"muscle",ssm,xxm,n=20)
maplot(dges,"synovium")
make_volcano(dges,"synovium")
make_heatmap(dges,"synovium",sss,xxs,n=20)
Pathway enrichment
Here I’m using the mitch package and gene pathways from Reactome to understand the affected pathways separately for each tissue type.
Reactome pathways downloaded Sept 2023.
#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
genesets <- gmt_import("ReactomePathways_2023-09-01.gmt")
# gene table
gt <- as.data.frame(rownames(xx))
gt$gene <- sapply(strsplit(gt[,1]," "),"[[",2)
m <- mitch_import(dgeb, DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 20241## Note: no. genes in output = 20208## Note: estimated proportion of input genes in output = 0.998mb <- mitch_calc(m, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(mb$enrichment_result,20) %>% kbl(caption = "Top gene pathway correlates with QuickDASH in bone") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 206 | Classical antibody-mediated complement activation | 73 | 0 | -0.8232098 | 0 |
| 392 | FCGR activation | 80 | 0 | -0.7777213 | 0 |
| 1168 | Scavenging of heme from plasma | 74 | 0 | -0.8020023 | 0 |
| 580 | Initial triggering of complement | 86 | 0 | -0.7422825 | 0 |
| 234 | Creation of C4 and C2 activators | 78 | 0 | -0.7760098 | 0 |
| 131 | CD22 mediated BCR regulation | 61 | 0 | -0.7916283 | 0 |
| 1119 | Role of LAT2/NTAL/LAB on calcium mobilization | 81 | 0 | -0.6868502 | 0 |
| 1120 | Role of phospholipids in phagocytosis | 93 | 0 | -0.6394591 | 0 |
| 1046 | Regulation of Complement cascade | 103 | 0 | -0.5717725 | 0 |
| 218 | Complement cascade | 110 | 0 | -0.5526231 | 0 |
| 563 | Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell | 183 | 0 | -0.4221616 | 0 |
| 390 | FCERI mediated MAPK activation | 97 | 0 | -0.5722277 | 0 |
| 118 | Binding and Uptake of Ligands by Scavenger Receptors | 102 | 0 | -0.5575587 | 0 |
| 389 | FCERI mediated Ca+2 mobilization | 96 | 0 | -0.5598036 | 0 |
| 78 | Antigen activates B Cell Receptor (BCR) leading to generation of second messengers | 86 | 0 | -0.5858738 | 0 |
| 394 | FCGR3A-mediated phagocytosis | 127 | 0 | -0.4822991 | 0 |
| 652 | Leishmania phagocytosis | 127 | 0 | -0.4822991 | 0 |
| 869 | Parasite infection | 127 | 0 | -0.4822991 | 0 |
| 1077 | Regulation of actin dynamics for phagocytic cup formation | 129 | 0 | -0.4709487 | 0 |
| 393 | FCGR3A-mediated IL10 synthesis | 104 | 0 | -0.5240365 | 0 |
write.table(mb$enrichment_result,file="mitch_bone.tsv",quote=FALSE,sep='\t')
m <- mitch_import(dgec, DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 18726## Note: no. genes in output = 18694## Note: estimated proportion of input genes in output = 0.998mc <- mitch_calc(m, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(mc$enrichment_result,20) %>% kbl(caption = "Top gene pathway correlates with QuickDASH in capsule") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 579 | Initial triggering of complement | 50 | 0 | -0.8339519 | 0 |
| 207 | Classical antibody-mediated complement activation | 38 | 0 | -0.9170438 | 0 |
| 235 | Creation of C4 and C2 activators | 42 | 0 | -0.8687693 | 0 |
| 171 | Cell Cycle | 609 | 0 | 0.2265186 | 0 |
| 173 | Cell Cycle, Mitotic | 493 | 0 | 0.2383106 | 0 |
| 219 | Complement cascade | 73 | 0 | -0.6063172 | 0 |
| 1041 | Regulation of Complement cascade | 67 | 0 | -0.6273008 | 0 |
| 1161 | Scavenging of heme from plasma | 39 | 0 | -0.8067281 | 0 |
| 119 | Binding and Uptake of Ligands by Scavenger Receptors | 66 | 0 | -0.6090489 | 0 |
| 172 | Cell Cycle Checkpoints | 250 | 0 | 0.3121570 | 0 |
| 393 | FCGR activation | 45 | 0 | -0.7009479 | 0 |
| 132 | CD22 mediated BCR regulation | 36 | 0 | -0.7384024 | 0 |
| 660 | M Phase | 352 | 0 | 0.2363310 | 0 |
| 732 | Mitotic Prometaphase | 193 | 0 | 0.3136593 | 0 |
| 1113 | Role of LAT2/NTAL/LAB on calcium mobilization | 45 | 0 | -0.6377881 | 0 |
| 1100 | Resolution of Sister Chromatid Cohesion | 116 | 0 | 0.3889760 | 0 |
| 1177 | Separation of Sister Chromatids | 179 | 0 | 0.2984508 | 0 |
| 734 | Mitotic Spindle Checkpoint | 108 | 0 | 0.3758883 | 0 |
| 728 | Mitotic Anaphase | 222 | 0 | 0.2612593 | 0 |
| 731 | Mitotic Metaphase and Anaphase | 223 | 0 | 0.2604910 | 0 |
write.table(mc$enrichment_result,file="mitch_capsule.tsv",quote=FALSE,sep='\t')
m <- mitch_import(dgef, DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 19072## Note: no. genes in output = 19042## Note: estimated proportion of input genes in output = 0.998mf <- mitch_calc(m, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(mf$enrichment_result,20) %>% kbl(caption = "Top gene pathway correlates with QuickDASH in fat") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 699 | Metabolism of RNA | 717 | 0 | -0.2141551 | 0 |
| 587 | Initial triggering of complement | 67 | 0 | -0.6617576 | 0 |
| 207 | Classical antibody-mediated complement activation | 53 | 0 | -0.7390803 | 0 |
| 1502 | rRNA processing in the nucleus and cytosol | 190 | 0 | -0.3747719 | 0 |
| 1052 | Regulation of Complement cascade | 79 | 0 | -0.5769910 | 0 |
| 235 | Creation of C4 and C2 activators | 59 | 0 | -0.6660920 | 0 |
| 159 | Cap-dependent Translation Initiation | 118 | 0 | -0.4698177 | 0 |
| 383 | Eukaryotic Translation Initiation | 118 | 0 | -0.4698177 | 0 |
| 1458 | Viral Infection Pathways | 714 | 0 | -0.1926770 | 0 |
| 119 | Binding and Uptake of Ligands by Scavenger Receptors | 82 | 0 | -0.5585224 | 0 |
| 687 | Major pathway of rRNA processing in the nucleolus and cytosol | 180 | 0 | -0.3771975 | 0 |
| 1407 | Translation | 294 | 0 | -0.2948276 | 0 |
| 433 | Formation of a pool of free 40S subunits | 100 | 0 | -0.5013124 | 0 |
| 579 | Infectious disease | 894 | 0 | -0.1707816 | 0 |
| 382 | Eukaryotic Translation Elongation | 93 | 0 | -0.5171226 | 0 |
| 219 | Complement cascade | 87 | 0 | -0.5314238 | 0 |
| 646 | L13a-mediated translational silencing of Ceruloplasmin expression | 110 | 0 | -0.4716056 | 0 |
| 711 | Metabolism of proteins | 1688 | 0 | -0.1249480 | 0 |
| 879 | Peptide chain elongation | 88 | 0 | -0.5209634 | 0 |
| 481 | GTP hydrolysis and joining of the 60S ribosomal subunit | 111 | 0 | -0.4540515 | 0 |
write.table(mf$enrichment_result,file="mitch_fat.tsv",quote=FALSE,sep='\t')
m <- mitch_import(dgem, DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 15920## Note: no. genes in output = 15898## Note: estimated proportion of input genes in output = 0.999mm <- mitch_calc(m, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(mm$enrichment_result,20) %>% kbl(caption = "Top gene pathway correlates with QuickDASH in muscle") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 202 | Classical antibody-mediated complement activation | 49 | 0 | -0.8114089 | 0 |
| 230 | Creation of C4 and C2 activators | 53 | 0 | -0.7680001 | 0 |
| 1133 | Scavenging of heme from plasma | 49 | 0 | -0.7959518 | 0 |
| 384 | FCGR activation | 54 | 0 | -0.7529477 | 0 |
| 129 | CD22 mediated BCR regulation | 46 | 0 | -0.7752197 | 0 |
| 566 | Initial triggering of complement | 60 | 0 | -0.6569474 | 0 |
| 1015 | Regulation of Complement cascade | 73 | 0 | -0.5906996 | 0 |
| 1085 | Role of LAT2/NTAL/LAB on calcium mobilization | 56 | 0 | -0.6687535 | 0 |
| 214 | Complement cascade | 78 | 0 | -0.5595481 | 0 |
| 116 | Binding and Uptake of Ligands by Scavenger Receptors | 75 | 0 | -0.5564314 | 0 |
| 672 | Metabolism of RNA | 714 | 0 | 0.1780957 | 0 |
| 1086 | Role of phospholipids in phagocytosis | 65 | 0 | -0.5774045 | 0 |
| 77 | Antigen activates B Cell Receptor (BCR) leading to generation of second messengers | 69 | 0 | -0.5554243 | 0 |
| 381 | FCERI mediated Ca+2 mobilization | 70 | 0 | -0.5335915 | 0 |
| 385 | FCGR3A-mediated IL10 synthesis | 77 | 0 | -0.4871176 | 0 |
| 382 | FCERI mediated MAPK activation | 72 | 0 | -0.4854003 | 0 |
| 386 | FCGR3A-mediated phagocytosis | 102 | 0 | -0.3995340 | 0 |
| 635 | Leishmania phagocytosis | 102 | 0 | -0.3995340 | 0 |
| 844 | Parasite infection | 102 | 0 | -0.3995340 | 0 |
| 1046 | Regulation of actin dynamics for phagocytic cup formation | 103 | 0 | -0.3648033 | 0 |
write.table(mm$enrichment_result,file="mitch_muscle.tsv",quote=FALSE,sep='\t')
m <- mitch_import(dges, DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 19361## Note: no. genes in output = 19327## Note: estimated proportion of input genes in output = 0.998ms <- mitch_calc(m, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(ms$enrichment_result,20) %>% kbl(caption = "Top gene pathway correlates with QuickDASH in synovium") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 692 | Metabolism | 1826 | 0 | -0.1941481 | 0 |
| 585 | Innate Immune System | 943 | 0 | -0.2519377 | 0 |
| 706 | Metabolism of proteins | 1693 | 0 | -0.1710352 | 0 |
| 311 | Disease | 1559 | 0 | -0.1737713 | 0 |
| 565 | Immune System | 1770 | 0 | -0.1618208 | 0 |
| 208 | Classical antibody-mediated complement activation | 45 | 0 | -0.9327110 | 0 |
| 574 | Infectious disease | 886 | 0 | -0.2142286 | 0 |
| 582 | Initial triggering of complement | 58 | 0 | -0.7991593 | 0 |
| 395 | FCGR activation | 52 | 0 | -0.8438072 | 0 |
| 1453 | Viral Infection Pathways | 706 | 0 | -0.2277259 | 0 |
| 236 | Creation of C4 and C2 activators | 50 | 0 | -0.8391181 | 0 |
| 133 | CD22 mediated BCR regulation | 43 | 0 | -0.9024471 | 0 |
| 1045 | Regulation of Complement cascade | 75 | 0 | -0.6715327 | 0 |
| 796 | Neutrophil degranulation | 439 | 0 | -0.2797162 | 0 |
| 1404 | Translation | 294 | 0 | -0.3342605 | 0 |
| 220 | Complement cascade | 82 | 0 | -0.6258363 | 0 |
| 1118 | Role of LAT2/NTAL/LAB on calcium mobilization | 52 | 0 | -0.7725332 | 0 |
| 1167 | Scavenging of heme from plasma | 46 | 0 | -0.8196197 | 0 |
| 397 | FCGR3A-mediated phagocytosis | 99 | 0 | -0.5505103 | 0 |
| 653 | Leishmania phagocytosis | 99 | 0 | -0.5505103 | 0 |
write.table(ms$enrichment_result,file="mitch_synovium.tsv",quote=FALSE,sep='\t')Mitch 1D reports
Here I’m using the mitch package to identify the most enriched pathways across all tissues.
There is a heatmap showing the QuickDASH association correlation across tissues.
Then there is the results of the multi-contrast pathway analysis with mitch in the form of a table and heatmap of pathway enrichment.
x <- list("bone"=dgeb, "capsule"=dgec, "fat"=dgef , "muscle"=dgem, "synovium"=dges)
y <- mitch_import(x, DEtype="deseq2",geneTable=gt)## Note: Mean no. genes in input = 18664## Note: no. genes in output = 14927## Note: estimated proportion of input genes in output = 0.8heatmap.2(cor(y,method="s"),trace="none",scale="none",margins = c(10,10),main="Spearman correlation across tissues")
res <- mitch_calc(y, genesets, priority="significance")## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(res$enrichment_result,30) %>%
kbl(caption = "Top gene pathway correlates with QuickDASH across all tissues - prioritised by p-value") %>%
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.bone | s.capsule | s.fat | s.muscle | s.synovium | p.bone | p.capsule | p.fat | p.muscle | p.synovium | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 225 | Creation of C4 and C2 activators | 40 | 0 | -0.8338282 | -0.8742729 | -0.6907470 | -0.7715893 | -0.8583630 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.8080074 | 0.0752405 | 0 |
| 197 | Classical antibody-mediated complement activation | 37 | 0 | -0.8378632 | -0.9216234 | -0.7389432 | -0.7862124 | -0.9196776 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.8871615 | 0.0808043 | 0 |
| 559 | Initial triggering of complement | 47 | 0 | -0.7695236 | -0.8640872 | -0.7029027 | -0.6286519 | -0.8120253 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.6992568 | 0.0922494 | 0 |
| 379 | FCGR activation | 42 | 0 | -0.7698450 | -0.7425436 | -0.5179071 | -0.7138346 | -0.8512309 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.6267803 | 0.1235815 | 0 |
| 1123 | Scavenging of heme from plasma | 37 | 0 | -0.8241192 | -0.8112319 | -0.6084911 | -0.7658614 | -0.8069882 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.7161890 | 0.0892713 | 0 |
| 209 | Complement cascade | 64 | 0 | -0.5556415 | -0.6705767 | -0.5328732 | -0.5138746 | -0.6687601 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.3242627 | 0.0756986 | 0 |
| 1005 | Regulation of Complement cascade | 60 | 0 | -0.5417838 | -0.6839578 | -0.5471088 | -0.5257461 | -0.6978812 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.3505936 | 0.0841537 | 0 |
| 124 | CD22 mediated BCR regulation | 35 | 0 | -0.8250297 | -0.7362534 | -0.6167607 | -0.7391658 | -0.8807030 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.7103043 | 0.1003901 | 0 |
| 368 | Eukaryotic Translation Elongation | 92 | 0 | 0.4991633 | -0.2787430 | -0.4858765 | 0.3305168 | -0.2194766 | 0.0000000 | 0.0000039 | 0.0000000 | 0.0000000 | 0.0002775 | 0.8487337 | 0.4229598 | 0 |
| 411 | Formation of a pool of free 40S subunits | 99 | 0 | 0.4491721 | -0.2673171 | -0.4695021 | 0.3197268 | -0.1838400 | 0.0000000 | 0.0000044 | 0.0000000 | 0.0000000 | 0.0015887 | 0.7935166 | 0.3953044 | 0 |
| 841 | Peptide chain elongation | 87 | 0 | 0.4880240 | -0.2850838 | -0.4902903 | 0.3151842 | -0.2080940 | 0.0000000 | 0.0000044 | 0.0000000 | 0.0000004 | 0.0008011 | 0.8381342 | 0.4171241 | 0 |
| 562 | Innate Immune System | 880 | 0 | -0.1072497 | 0.0185628 | -0.0381631 | -0.0395685 | -0.2096494 | 0.0000001 | 0.3548771 | 0.0571536 | 0.0485880 | 0.0000000 | 0.2425326 | 0.0873737 | 0 |
| 617 | L13a-mediated translational silencing of Ceruloplasmin expression | 109 | 0 | 0.4182466 | -0.2189118 | -0.4377443 | 0.2994882 | -0.1539437 | 0.0000000 | 0.0000799 | 0.0000000 | 0.0000001 | 0.0055418 | 0.7265429 | 0.3626775 | 0 |
| 764 | Neutrophil degranulation | 422 | 0 | -0.1400896 | 0.1102937 | -0.0086007 | -0.0573543 | -0.2400114 | 0.0000009 | 0.0001092 | 0.7629265 | 0.0442587 | 0.0000000 | 0.3045632 | 0.1324873 | 0 |
| 111 | Binding and Uptake of Ligands by Scavenger Receptors | 63 | 0 | -0.4507065 | -0.6017586 | -0.5271178 | -0.4918414 | -0.5404621 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.1735034 | 0.0563916 | 0 |
| 456 | GTP hydrolysis and joining of the 60S ribosomal subunit | 110 | 0 | 0.4184469 | -0.2275556 | -0.4190383 | 0.2968004 | -0.1601576 | 0.0000000 | 0.0000380 | 0.0000000 | 0.0000001 | 0.0037456 | 0.7184799 | 0.3586568 | 0 |
| 1403 | Viral mRNA Translation | 87 | 0 | 0.4775459 | -0.2869350 | -0.4580165 | 0.3046349 | -0.2031214 | 0.0000000 | 0.0000038 | 0.0000000 | 0.0000009 | 0.0010669 | 0.8088397 | 0.4027150 | 0 |
| 149 | Cap-dependent Translation Initiation | 117 | 0 | 0.3779255 | -0.2055356 | -0.4359482 | 0.2979726 | -0.1628006 | 0.0000000 | 0.0001249 | 0.0000000 | 0.0000000 | 0.0023791 | 0.7002965 | 0.3489694 | 0 |
| 369 | Eukaryotic Translation Initiation | 117 | 0 | 0.3779255 | -0.2055356 | -0.4359482 | 0.2979726 | -0.1628006 | 0.0000000 | 0.0001249 | 0.0000000 | 0.0000000 | 0.0023791 | 0.7002965 | 0.3489694 | 0 |
| 370 | Eukaryotic Translation Termination | 91 | 0 | 0.4559254 | -0.2874416 | -0.4589312 | 0.3003105 | -0.2215060 | 0.0000000 | 0.0000022 | 0.0000000 | 0.0000007 | 0.0002630 | 0.8002249 | 0.3973038 | 0 |
| 1075 | Role of LAT2/NTAL/LAB on calcium mobilization | 44 | 0 | -0.5909763 | -0.6329094 | -0.4282250 | -0.6076854 | -0.7198237 | 0.0000000 | 0.0000000 | 0.0000009 | 0.0000000 | 0.0000000 | 1.3493079 | 0.1060717 | 0 |
| 1127 | Selenocysteine synthesis | 91 | 0 | 0.4587623 | -0.2733816 | -0.4583638 | 0.2926902 | -0.1903182 | 0.0000000 | 0.0000067 | 0.0000000 | 0.0000014 | 0.0017172 | 0.7856121 | 0.3909491 | 0 |
| 1355 | Translation | 293 | 0 | 0.2015718 | -0.0993684 | -0.2480786 | 0.2184846 | -0.2794040 | 0.0000000 | 0.0035321 | 0.0000000 | 0.0000000 | 0.0000000 | 0.4876991 | 0.2394251 | 0 |
| 161 | Cell Cycle | 557 | 0 | -0.1808123 | 0.2022281 | 0.0001157 | 0.0156525 | -0.0007756 | 0.0000000 | 0.0000000 | 0.9962980 | 0.5301736 | 0.9751856 | 0.2717258 | 0.1356188 | 0 |
| 772 | Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) | 93 | 0 | 0.4450384 | -0.2751656 | -0.4453950 | 0.2952879 | -0.2036893 | 0.0000000 | 0.0000046 | 0.0000000 | 0.0000009 | 0.0006938 | 0.7751363 | 0.3853799 | 0 |
| 676 | Metabolism of proteins | 1575 | 0 | 0.0914071 | 0.0127111 | -0.0920452 | 0.0508113 | -0.1318524 | 0.0000000 | 0.4086317 | 0.0000000 | 0.0009544 | 0.0000000 | 0.1922394 | 0.0948745 | 0 |
| 1106 | SRP-dependent cotranslational protein targeting to membrane | 110 | 0 | 0.4149313 | -0.2335033 | -0.3891721 | 0.2234583 | -0.2586501 | 0.0000000 | 0.0000237 | 0.0000000 | 0.0000523 | 0.0000028 | 0.7035483 | 0.3475545 | 0 |
| 1076 | Role of phospholipids in phagocytosis | 53 | 0 | -0.5413499 | -0.5436839 | -0.4366160 | -0.5052479 | -0.6212233 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.1918386 | 0.0670267 | 0 |
| 662 | Metabolism | 1683 | 0 | 0.0262351 | -0.0193113 | -0.0643987 | 0.0190041 | -0.1691716 | 0.0791085 | 0.1962020 | 0.0000162 | 0.2034139 | 0.0000000 | 0.1849015 | 0.0799461 | 0 |
| 1067 | Response of EIF2AK4 (GCN2) to amino acid deficiency | 99 | 0 | 0.4403776 | -0.1961550 | -0.3991357 | 0.2762423 | -0.1512658 | 0.0000000 | 0.0007528 | 0.0000000 | 0.0000021 | 0.0093691 | 0.7006495 | 0.3502608 | 0 |
pw <- res$enrichment_result[1:30,4:8]
rownames(pw) <- res$enrichment_result[1:30,1]
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(as.matrix(pw),trace="none",col=colfunc(25),scale="none", margins = c(10,25), cexRow=0.6,cexCol=0.7,
main="Top ranked pathways (p-val)")
res2 <- mitch_calc(y, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(subset(res2$enrichment_result,p.adjustMANOVA<0.05),30) %>%
kbl(caption = "Top gene pathway correlates with QuickDASH across all tissues - prioritised by effect size") %>%
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.bone | s.capsule | s.fat | s.muscle | s.synovium | p.bone | p.capsule | p.fat | p.muscle | p.synovium | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 197 | Classical antibody-mediated complement activation | 37 | 0.0000000 | -0.8378632 | -0.9216234 | -0.7389432 | -0.7862124 | -0.9196776 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.8871615 | 0.0808043 | 0.0000000 |
| 225 | Creation of C4 and C2 activators | 40 | 0.0000000 | -0.8338282 | -0.8742729 | -0.6907470 | -0.7715893 | -0.8583630 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.8080074 | 0.0752405 | 0.0000000 |
| 1123 | Scavenging of heme from plasma | 37 | 0.0000000 | -0.8241192 | -0.8112319 | -0.6084911 | -0.7658614 | -0.8069882 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.7161890 | 0.0892713 | 0.0000000 |
| 124 | CD22 mediated BCR regulation | 35 | 0.0000000 | -0.8250297 | -0.7362534 | -0.6167607 | -0.7391658 | -0.8807030 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.7103043 | 0.1003901 | 0.0000000 |
| 559 | Initial triggering of complement | 47 | 0.0000000 | -0.7695236 | -0.8640872 | -0.7029027 | -0.6286519 | -0.8120253 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.6992568 | 0.0922494 | 0.0000000 |
| 379 | FCGR activation | 42 | 0.0000000 | -0.7698450 | -0.7425436 | -0.5179071 | -0.7138346 | -0.8512309 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.6267803 | 0.1235815 | 0.0000000 |
| 1005 | Regulation of Complement cascade | 60 | 0.0000000 | -0.5417838 | -0.6839578 | -0.5471088 | -0.5257461 | -0.6978812 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.3505936 | 0.0841537 | 0.0000000 |
| 1075 | Role of LAT2/NTAL/LAB on calcium mobilization | 44 | 0.0000000 | -0.5909763 | -0.6329094 | -0.4282250 | -0.6076854 | -0.7198237 | 0.0000000 | 0.0000000 | 0.0000009 | 0.0000000 | 0.0000000 | 1.3493079 | 0.1060717 | 0.0000000 |
| 209 | Complement cascade | 64 | 0.0000000 | -0.5556415 | -0.6705767 | -0.5328732 | -0.5138746 | -0.6687601 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.3242627 | 0.0756986 | 0.0000000 |
| 1076 | Role of phospholipids in phagocytosis | 53 | 0.0000000 | -0.5413499 | -0.5436839 | -0.4366160 | -0.5052479 | -0.6212233 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.1918386 | 0.0670267 | 0.0000000 |
| 111 | Binding and Uptake of Ligands by Scavenger Receptors | 63 | 0.0000000 | -0.4507065 | -0.6017586 | -0.5271178 | -0.4918414 | -0.5404621 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 1.1735034 | 0.0563916 | 0.0000000 |
| 73 | Antigen activates B Cell Receptor (BCR) leading to generation of second messengers | 58 | 0.0000000 | -0.5297715 | -0.4545606 | -0.3073045 | -0.4907526 | -0.5633521 | 0.0000000 | 0.0000000 | 0.0000520 | 0.0000000 | 0.0000000 | 1.0676717 | 0.0992755 | 0.0000000 |
| 138 | CS/DS degradation | 10 | 0.0000014 | 0.6193605 | -0.1824227 | -0.5407790 | -0.3636388 | -0.4608567 | 0.0006948 | 0.3179067 | 0.0030648 | 0.0464776 | 0.0116217 | 1.0266204 | 0.4694645 | 0.0000090 |
| 1449 | tRNA processing in the mitochondrion | 23 | 0.0000000 | 0.7568263 | -0.2186574 | 0.3984107 | 0.4068940 | -0.0216983 | 0.0000000 | 0.0695465 | 0.0009422 | 0.0007310 | 0.8570853 | 0.9722971 | 0.3859870 | 0.0000000 |
| 380 | FCGR3A-mediated IL10 synthesis | 65 | 0.0000000 | -0.4003064 | -0.4674989 | -0.2951999 | -0.4105400 | -0.4953656 | 0.0000000 | 0.0000000 | 0.0000389 | 0.0000000 | 0.0000000 | 0.9380158 | 0.0771285 | 0.0000000 |
| 376 | FCERI mediated Ca+2 mobilization | 58 | 0.0000000 | -0.4159104 | -0.3873762 | -0.2732902 | -0.4579813 | -0.4571279 | 0.0000000 | 0.0000003 | 0.0003202 | 0.0000000 | 0.0000000 | 0.9035719 | 0.0759597 | 0.0000000 |
| 612 | Keratan sulfate degradation | 12 | 0.0000036 | 0.5497262 | 0.0293105 | -0.2762432 | -0.4061236 | -0.4913063 | 0.0009760 | 0.8604684 | 0.0975741 | 0.0148589 | 0.0032105 | 0.8863896 | 0.4227789 | 0.0000218 |
| 377 | FCERI mediated MAPK activation | 60 | 0.0000000 | -0.4395462 | -0.4314657 | -0.2388332 | -0.4024596 | -0.4299209 | 0.0000000 | 0.0000000 | 0.0013833 | 0.0000001 | 0.0000000 | 0.8849918 | 0.0847975 | 0.0000000 |
| 518 | Heme biosynthesis | 12 | 0.0000147 | -0.8333557 | 0.1049726 | -0.1106716 | -0.0865683 | -0.1790926 | 0.0000006 | 0.5289995 | 0.5068748 | 0.6036491 | 0.2827982 | 0.8702398 | 0.3582032 | 0.0000797 |
| 727 | NGF-stimulated transcription | 33 | 0.0000000 | 0.4084948 | 0.3189773 | 0.3918967 | 0.3397626 | 0.4381345 | 0.0000489 | 0.0015215 | 0.0000979 | 0.0007322 | 0.0000133 | 0.8541651 | 0.0491793 | 0.0000000 |
| 796 | Olfactory Signaling Pathway | 12 | 0.0000050 | -0.4127724 | -0.3156554 | 0.2614706 | -0.2497039 | 0.5666443 | 0.0132981 | 0.0583411 | 0.1168518 | 0.1342497 | 0.0006766 | 0.8496023 | 0.4234746 | 0.0000296 |
| 368 | Eukaryotic Translation Elongation | 92 | 0.0000000 | 0.4991633 | -0.2787430 | -0.4858765 | 0.3305168 | -0.2194766 | 0.0000000 | 0.0000039 | 0.0000000 | 0.0000000 | 0.0002775 | 0.8487337 | 0.4229598 | 0.0000000 |
| 381 | FCGR3A-mediated phagocytosis | 90 | 0.0000000 | -0.3701767 | -0.3616469 | -0.3180832 | -0.3311990 | -0.4899014 | 0.0000000 | 0.0000000 | 0.0000002 | 0.0000001 | 0.0000000 | 0.8477572 | 0.0681153 | 0.0000000 |
| 628 | Leishmania phagocytosis | 90 | 0.0000000 | -0.3701767 | -0.3616469 | -0.3180832 | -0.3311990 | -0.4899014 | 0.0000000 | 0.0000000 | 0.0000002 | 0.0000001 | 0.0000000 | 0.8477572 | 0.0681153 | 0.0000000 |
| 836 | Parasite infection | 90 | 0.0000000 | -0.3701767 | -0.3616469 | -0.3180832 | -0.3311990 | -0.4899014 | 0.0000000 | 0.0000000 | 0.0000002 | 0.0000001 | 0.0000000 | 0.8477572 | 0.0681153 | 0.0000000 |
| 1021 | Regulation of PTEN mRNA translation | 11 | 0.0184479 | 0.0686145 | 0.4350057 | 0.3791048 | -0.2659865 | 0.5558522 | 0.6935934 | 0.0124879 | 0.0294821 | 0.1266753 | 0.0014118 | 0.8469822 | 0.3325602 | 0.0473769 |
| 1443 | rRNA processing in the mitochondrion | 24 | 0.0000000 | 0.5328066 | -0.2006531 | 0.3328972 | 0.5298374 | 0.0075488 | 0.0000062 | 0.0889113 | 0.0047637 | 0.0000070 | 0.9489727 | 0.8460199 | 0.3265650 | 0.0000000 |
| 841 | Peptide chain elongation | 87 | 0.0000000 | 0.4880240 | -0.2850838 | -0.4902903 | 0.3151842 | -0.2080940 | 0.0000000 | 0.0000044 | 0.0000000 | 0.0000004 | 0.0008011 | 0.8381342 | 0.4171241 | 0.0000000 |
| 1036 | Regulation of actin dynamics for phagocytic cup formation | 91 | 0.0000000 | -0.3597427 | -0.3542764 | -0.3471093 | -0.2923569 | -0.4795197 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000015 | 0.0000000 | 0.8311577 | 0.0686325 | 0.0000000 |
| 263 | Defective B3GALT6 causes EDSP2 and SEMDJL1 | 15 | 0.0000054 | 0.5086463 | -0.3806509 | -0.4287554 | -0.0258941 | -0.3049803 | 0.0006478 | 0.0107025 | 0.0040417 | 0.8621787 | 0.0408732 | 0.8253065 | 0.3877301 | 0.0000320 |
pw <- subset(res2$enrichment_result,p.adjustMANOVA<0.05)[1:30,4:8]
rownames(pw) <- res2$enrichment_result[1:30,1]
heatmap.2(as.matrix(pw),trace="none",col=colfunc(25),scale="none", margins = c(10,25), cexRow=0.6,cexCol=0.7,
main="Top ranked pathways (ES)")
res3 <- mitch_calc(y, genesets, priority="SD")## Note: Prioritisation by SD after selecting sets with
## p.adjustMANOVA<=0.05.head(subset(res3$enrichment_result,p.adjustMANOVA<0.05),30) %>%
kbl(caption = "Top gene pathway correlates with QuickDASH across all tissues - prioritised by standard deviation") %>%
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.bone | s.capsule | s.fat | s.muscle | s.synovium | p.bone | p.capsule | p.fat | p.muscle | p.synovium | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 138 | CS/DS degradation | 10 | 0.0000014 | 0.6193605 | -0.1824227 | -0.5407790 | -0.3636388 | -0.4608567 | 0.0006948 | 0.3179067 | 0.0030648 | 0.0464776 | 0.0116217 | 1.0266204 | 0.4694645 | 0.0000090 |
| 796 | Olfactory Signaling Pathway | 12 | 0.0000050 | -0.4127724 | -0.3156554 | 0.2614706 | -0.2497039 | 0.5666443 | 0.0132981 | 0.0583411 | 0.1168518 | 0.1342497 | 0.0006766 | 0.8496023 | 0.4234746 | 0.0000296 |
| 368 | Eukaryotic Translation Elongation | 92 | 0.0000000 | 0.4991633 | -0.2787430 | -0.4858765 | 0.3305168 | -0.2194766 | 0.0000000 | 0.0000039 | 0.0000000 | 0.0000000 | 0.0002775 | 0.8487337 | 0.4229598 | 0.0000000 |
| 612 | Keratan sulfate degradation | 12 | 0.0000036 | 0.5497262 | 0.0293105 | -0.2762432 | -0.4061236 | -0.4913063 | 0.0009760 | 0.8604684 | 0.0975741 | 0.0148589 | 0.0032105 | 0.8863896 | 0.4227789 | 0.0000218 |
| 841 | Peptide chain elongation | 87 | 0.0000000 | 0.4880240 | -0.2850838 | -0.4902903 | 0.3151842 | -0.2080940 | 0.0000000 | 0.0000044 | 0.0000000 | 0.0000004 | 0.0008011 | 0.8381342 | 0.4171241 | 0.0000000 |
| 866 | Polo-like kinase mediated events | 14 | 0.0000021 | -0.5464552 | 0.5015279 | 0.0752843 | -0.3162054 | 0.0673334 | 0.0003998 | 0.0011579 | 0.6258119 | 0.0405349 | 0.6627421 | 0.8126074 | 0.4033550 | 0.0000129 |
| 1403 | Viral mRNA Translation | 87 | 0.0000000 | 0.4775459 | -0.2869350 | -0.4580165 | 0.3046349 | -0.2031214 | 0.0000000 | 0.0000038 | 0.0000000 | 0.0000009 | 0.0010669 | 0.8088397 | 0.4027150 | 0.0000000 |
| 370 | Eukaryotic Translation Termination | 91 | 0.0000000 | 0.4559254 | -0.2874416 | -0.4589312 | 0.3003105 | -0.2215060 | 0.0000000 | 0.0000022 | 0.0000000 | 0.0000007 | 0.0002630 | 0.8002249 | 0.3973038 | 0.0000000 |
| 411 | Formation of a pool of free 40S subunits | 99 | 0.0000000 | 0.4491721 | -0.2673171 | -0.4695021 | 0.3197268 | -0.1838400 | 0.0000000 | 0.0000044 | 0.0000000 | 0.0000000 | 0.0015887 | 0.7935166 | 0.3953044 | 0.0000000 |
| 1127 | Selenocysteine synthesis | 91 | 0.0000000 | 0.4587623 | -0.2733816 | -0.4583638 | 0.2926902 | -0.1903182 | 0.0000000 | 0.0000067 | 0.0000000 | 0.0000014 | 0.0017172 | 0.7856121 | 0.3909491 | 0.0000000 |
| 263 | Defective B3GALT6 causes EDSP2 and SEMDJL1 | 15 | 0.0000054 | 0.5086463 | -0.3806509 | -0.4287554 | -0.0258941 | -0.3049803 | 0.0006478 | 0.0107025 | 0.0040417 | 0.8621787 | 0.0408732 | 0.8253065 | 0.3877301 | 0.0000320 |
| 1081 | SARS-CoV-1 modulates host translation machinery | 36 | 0.0000000 | 0.4585395 | -0.2264418 | -0.4791410 | 0.2627799 | -0.2161858 | 0.0000019 | 0.0187516 | 0.0000007 | 0.0063777 | 0.0248350 | 0.7790379 | 0.3869315 | 0.0000000 |
| 1449 | tRNA processing in the mitochondrion | 23 | 0.0000000 | 0.7568263 | -0.2186574 | 0.3984107 | 0.4068940 | -0.0216983 | 0.0000000 | 0.0695465 | 0.0009422 | 0.0007310 | 0.8570853 | 0.9722971 | 0.3859870 | 0.0000000 |
| 772 | Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) | 93 | 0.0000000 | 0.4450384 | -0.2751656 | -0.4453950 | 0.2952879 | -0.2036893 | 0.0000000 | 0.0000046 | 0.0000000 | 0.0000009 | 0.0006938 | 0.7751363 | 0.3853799 | 0.0000000 |
| 301 | Diseases associated with glycosaminoglycan metabolism | 34 | 0.0000000 | 0.5610532 | -0.1117659 | -0.3610382 | -0.1263720 | -0.3298984 | 0.0000000 | 0.2595635 | 0.0002698 | 0.2023911 | 0.0008736 | 0.7631668 | 0.3726042 | 0.0000000 |
| 266 | Defective B4GALT7 causes EDS, progeroid type | 15 | 0.0000240 | 0.4819832 | -0.3505097 | -0.4143687 | -0.0270833 | -0.2932403 | 0.0012297 | 0.0187669 | 0.0054627 | 0.8559170 | 0.0492882 | 0.7833197 | 0.3677006 | 0.0001268 |
| 617 | L13a-mediated translational silencing of Ceruloplasmin expression | 109 | 0.0000000 | 0.4182466 | -0.2189118 | -0.4377443 | 0.2994882 | -0.1539437 | 0.0000000 | 0.0000799 | 0.0000000 | 0.0000001 | 0.0055418 | 0.7265429 | 0.3626775 | 0.0000000 |
| 456 | GTP hydrolysis and joining of the 60S ribosomal subunit | 110 | 0.0000000 | 0.4184469 | -0.2275556 | -0.4190383 | 0.2968004 | -0.1601576 | 0.0000000 | 0.0000380 | 0.0000000 | 0.0000001 | 0.0037456 | 0.7184799 | 0.3586568 | 0.0000000 |
| 518 | Heme biosynthesis | 12 | 0.0000147 | -0.8333557 | 0.1049726 | -0.1106716 | -0.0865683 | -0.1790926 | 0.0000006 | 0.5289995 | 0.5068748 | 0.6036491 | 0.2827982 | 0.8702398 | 0.3582032 | 0.0000797 |
| 265 | Defective B3GAT3 causes JDSSDHD | 15 | 0.0000407 | 0.4558834 | -0.3298194 | -0.4330562 | -0.0460658 | -0.3005007 | 0.0022368 | 0.0270096 | 0.0036874 | 0.7574446 | 0.0439282 | 0.7723800 | 0.3574714 | 0.0002080 |
| 1247 | Syndecan interactions | 26 | 0.0000000 | 0.6203610 | 0.3639714 | -0.1634790 | -0.1768234 | -0.0333225 | 0.0000000 | 0.0013182 | 0.1491674 | 0.1187044 | 0.7687459 | 0.7592268 | 0.3542007 | 0.0000000 |
| 268 | Defective EXT1 causes exostoses 1, TRPS2 and CHDS | 12 | 0.0022220 | 0.4765896 | -0.1671807 | -0.4021567 | 0.1323612 | -0.2592468 | 0.0042560 | 0.3160472 | 0.0158652 | 0.4273180 | 0.1199965 | 0.7081987 | 0.3506771 | 0.0078066 |
| 269 | Defective EXT2 causes exostoses 2 | 12 | 0.0022220 | 0.4765896 | -0.1671807 | -0.4021567 | 0.1323612 | -0.2592468 | 0.0042560 | 0.3160472 | 0.0158652 | 0.4273180 | 0.1199965 | 0.7081987 | 0.3506771 | 0.0078066 |
| 1067 | Response of EIF2AK4 (GCN2) to amino acid deficiency | 99 | 0.0000000 | 0.4403776 | -0.1961550 | -0.3991357 | 0.2762423 | -0.1512658 | 0.0000000 | 0.0007528 | 0.0000000 | 0.0000021 | 0.0093691 | 0.7006495 | 0.3502608 | 0.0000000 |
| 149 | Cap-dependent Translation Initiation | 117 | 0.0000000 | 0.3779255 | -0.2055356 | -0.4359482 | 0.2979726 | -0.1628006 | 0.0000000 | 0.0001249 | 0.0000000 | 0.0000000 | 0.0023791 | 0.7002965 | 0.3489694 | 0.0000000 |
| 369 | Eukaryotic Translation Initiation | 117 | 0.0000000 | 0.3779255 | -0.2055356 | -0.4359482 | 0.2979726 | -0.1628006 | 0.0000000 | 0.0001249 | 0.0000000 | 0.0000000 | 0.0023791 | 0.7002965 | 0.3489694 | 0.0000000 |
| 1106 | SRP-dependent cotranslational protein targeting to membrane | 110 | 0.0000000 | 0.4149313 | -0.2335033 | -0.3891721 | 0.2234583 | -0.2586501 | 0.0000000 | 0.0000237 | 0.0000000 | 0.0000523 | 0.0000028 | 0.7035483 | 0.3475545 | 0.0000000 |
| 1397 | VLDLR internalisation and degradation | 15 | 0.0051513 | 0.0439646 | -0.3315898 | -0.3123033 | 0.3951359 | -0.4425697 | 0.7681841 | 0.0261988 | 0.0362671 | 0.0080633 | 0.0030019 | 0.7492786 | 0.3455441 | 0.0157358 |
| 439 | G1/S-Specific Transcription | 23 | 0.0000001 | -0.4915692 | 0.4329623 | 0.0598147 | -0.2089255 | 0.0416871 | 0.0000449 | 0.0003254 | 0.6195870 | 0.0829056 | 0.7293549 | 0.6914204 | 0.3437110 | 0.0000005 |
| 237 | Cytosolic iron-sulfur cluster assembly | 12 | 0.0079977 | -0.3283160 | -0.0887473 | -0.5143144 | 0.3778299 | -0.2928595 | 0.0489454 | 0.5945681 | 0.0020363 | 0.0234469 | 0.0790211 | 0.7801981 | 0.3411114 | 0.0228438 |
pw <- subset(res3$enrichment_result,p.adjustMANOVA<0.05)[1:30,4:8]
rownames(pw) <- res3$enrichment_result[1:30,1]
heatmap.2(as.matrix(pw),trace="none",col=colfunc(25),scale="none", margins = c(10,25), cexRow=0.6,cexCol=0.7,
main="Top ranked pathways (SD)")
Mitch 5D report
Now I will look at the expression of the top ranked pathway genes (Creation of C4 and C2 activators) across the whole experiment. It will give us a good idea about what’s going on.
set <- head(res$enrichment_result$set,1)
genes <- genesets[[set]]
genes <- gt[match(genes,gt$gene),1]
xxn <- xx/rowSums(xx)*1000000
mx <- as.matrix(xxn[which(rownames(xxn) %in% genes),])
rownames(mx) <- sapply(strsplit(rownames(mx)," "),"[[",2)
colfunc <- colorRampPalette(c("blue", "white", "red"))
parti <- substr(colnames(mx), 1, 2)
values <- as.numeric(sss[match(parti,sss$ParticipantID),"quickdash"])
f <- colorRamp(c("yellow", "orange"))
rr <- range(values)
svals <- (values-rr[1])/diff(rr)
colcols <- rgb(f(svals)/255)
heatmap.2(mx,trace="none",col=colfunc(25),scale="row", margins = c(5,10),
cexRow=0.7, cexCol=0.7, main=set , ColSideColors = colcols )
mtext("yellow=low quickdash, orange=high quickdash")
It looks like expression of immunoglubulins in some of the bone samples could explain the observed pathway enrichment.
Conclusions
The number of differentially expressed genes was relatively small. This is to be expected given that the sample number (n) is very small. Still, there were some statistically significant (FDR<0.05) differentially expressed genes in each tissue type (except in muscle). Some of the top ranked genes appear to be involved in inflammation and response to stress including interleukins, chemokines, complement components, nuclear response factors and immediate early response genes. Some of these were up and down regulated, which suggests a shifting balance between inflammatory and antiinflammatory signaling pathways. Using a gene expression heatmap of one of the top pathways, we see it is dominated by expression of innumoglobulins in bone. This could be an important finding or it might be due to contamination of bone samples with blood.
Session information
For reproducibility.
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.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] beeswarm_0.4.0 gtools_3.9.5
## [5] echarts4r_0.4.5 kableExtra_1.4.0
## [7] topconfects_1.20.0 limma_3.60.3
## [9] eulerr_7.0.2 mitch_1.16.0
## [11] MASS_7.3-61 fgsea_1.30.0
## [13] gplots_3.1.3.1 DESeq2_1.44.0
## [15] SummarizedExperiment_1.34.0 Biobase_2.64.0
## [17] MatrixGenerics_1.16.0 matrixStats_1.3.0
## [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] lubridate_1.9.3 forcats_1.0.0
## [27] stringr_1.5.1 dplyr_1.1.4
## [29] purrr_1.0.2 readr_2.1.5
## [31] tidyr_1.3.1 tibble_3.2.1
## [33] ggplot2_3.5.1 tidyverse_2.0.0
## [35] zoo_1.8-12
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 gridExtra_2.3 rlang_1.1.4
## [4] magrittr_2.0.3 compiler_4.4.1 systemfonts_1.1.0
## [7] vctrs_0.6.5 pkgconfig_2.0.3 crayon_1.5.3
## [10] fastmap_1.2.0 XVector_0.44.0 labeling_0.4.3
## [13] caTools_1.18.2 utf8_1.2.4 promises_1.3.0
## [16] rmarkdown_2.27 tzdb_0.4.0 UCSC.utils_1.0.0
## [19] xfun_0.45 zlibbioc_1.50.0 cachem_1.1.0
## [22] jsonlite_1.8.8 progress_1.2.3 highr_0.11
## [25] later_1.3.2 DelayedArray_0.30.1 BiocParallel_1.38.0
## [28] prettyunits_1.2.0 parallel_4.4.1 R6_2.5.1
## [31] bslib_0.7.0 stringi_1.8.4 RColorBrewer_1.1-3
## [34] jquerylib_0.1.4 assertthat_0.2.1 Rcpp_1.0.12
## [37] knitr_1.48 httpuv_1.6.15 Matrix_1.7-0
## [40] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.16.0
## [43] abind_1.4-5 yaml_2.3.9 codetools_0.2-20
## [46] lattice_0.22-6 plyr_1.8.9 shiny_1.8.1.1
## [49] withr_3.0.0 evaluate_0.24.0 ggstats_0.6.0
## [52] xml2_1.3.6 pillar_1.9.0 KernSmooth_2.23-24
## [55] generics_0.1.3 hms_1.1.3 munsell_0.5.1
## [58] scales_1.3.0 xtable_1.8-4 glue_1.7.0
## [61] tools_4.4.1 data.table_1.15.4 locfit_1.5-9.10
## [64] fastmatch_1.1-4 cowplot_1.1.3 grid_4.4.1
## [67] colorspace_2.1-0 GenomeInfoDbData_1.2.12 cli_3.6.3
## [70] fansi_1.0.6 S4Arrays_1.4.1 viridisLite_0.4.2
## [73] svglite_2.1.3 gtable_0.3.5 sass_0.4.9
## [76] digest_0.6.36 SparseArray_1.4.8 farver_2.1.2
## [79] htmlwidgets_1.6.4 htmltools_0.5.8.1 lifecycle_1.0.4
## [82] httr_1.4.7 statmod_1.5.0 mime_0.12