Source: TBA

Introduction

Here we analyse the effect of two antimicrobial peptides on cancer cells.

Reads were mapped to the human transcriptome version 46 from GENCODE genes (gencode.v46.transcripts.fa).

Genes with mean counts > 10 are classified as detected.

Differential expression is conducted with DESeq2.

Pathway enrichment analysis is conducted with mitch.

Gene sets were obtained from the gene ontology database (q3 2023). Biological process sets were used.

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

Import read counts

tmp <- read.table("3col.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
xx$geneid = NULL
xx <- round(xx)
head(xx)
##                           control-1 control-2 control-3 P1-1 P1-2 P1-3 PA-1
## ENSG00000000003.16 TSPAN6      1737      1850      2001 1677 1614 1959 1979
## ENSG00000000005.6 TNMD            0         0         0    0    0    0    0
## ENSG00000000419.14 DPM1        1708      1717      1835 1750 1841 2027 1875
## ENSG00000000457.14 SCYL3        360       369       411  297  302  423  313
## ENSG00000000460.17 FIRRM        896       959      1157 1030  953 1219  942
## ENSG00000000938.13 FGR            3         1         0    1    0    0    0
##                           PA-2 PA-3
## ENSG00000000003.16 TSPAN6 2194 2205
## ENSG00000000005.6 TNMD       0    0
## ENSG00000000419.14 DPM1   1914 1825
## ENSG00000000457.14 SCYL3   269  283
## ENSG00000000460.17 FIRRM  1091 1180
## ENSG00000000938.13 FGR       0    0

Sample sheet

ss <- data.frame(colnames(xx))
ss$P1 <- as.numeric(grepl("P1",ss[,1]))
ss$PA <- as.numeric(grepl("PA",ss[,1]))
rownames(ss) <- ss[,1]
ss[,1] = NULL
ss
##           P1 PA
## control-1  0  0
## control-2  0  0
## control-3  0  0
## P1-1       1  0
## P1-2       1  0
## P1-3       1  0
## PA-1       0  1
## PA-2       0  1
## PA-3       0  1

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

colSums(xx)
## control-1 control-2 control-3      P1-1      P1-2      P1-3      PA-1      PA-2 
##  30590171  31289111  34035097  31852247  30788467  37034503  30507675  33557235 
##      PA-3 
##  34880005

MDS plot

All samples have sufficient number of reads.

All sample cluster with others from the same group.

P1-3 could be an outlier.

cols <- c(rep("lightgreen",3),rep("lightblue",3),rep("pink",3))

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=colnames(xx) )

Correlation heatmap

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

cor(xx) %>% kbl(caption = "Pearson correlation coefficients") %>% kable_paper("hover", full_width = F)
Pearson correlation coefficients
control-1 control-2 control-3 P1-1 P1-2 P1-3 PA-1 PA-2 PA-3
control-1 1.0000000 0.9959543 0.9961223 0.9576279 0.9681948 0.9895587 0.9440021 0.9139680 0.9605070
control-2 0.9959543 1.0000000 0.9950670 0.9683034 0.9756609 0.9822263 0.9603951 0.9361069 0.9723902
control-3 0.9961223 0.9950670 1.0000000 0.9691303 0.9782163 0.9904619 0.9557366 0.9291840 0.9707045
P1-1 0.9576279 0.9683034 0.9691303 1.0000000 0.9984858 0.9722928 0.9774447 0.9675656 0.9804198
P1-2 0.9681948 0.9756609 0.9782163 0.9984858 1.0000000 0.9811723 0.9783232 0.9655411 0.9832908
P1-3 0.9895587 0.9822263 0.9904619 0.9722928 0.9811723 1.0000000 0.9446017 0.9163373 0.9606774
PA-1 0.9440021 0.9603951 0.9557366 0.9774447 0.9783232 0.9446017 1.0000000 0.9954243 0.9967331
PA-2 0.9139680 0.9361069 0.9291840 0.9675656 0.9655411 0.9163373 0.9954243 1.0000000 0.9873859
PA-3 0.9605070 0.9723902 0.9707045 0.9804198 0.9832908 0.9606774 0.9967331 0.9873859 1.0000000 
cor(xx,method="spearman") %>% kbl(caption = "Spearman correlation coefficients") %>% kable_paper("hover", full_width = F)
Spearman correlation coefficients
control-1 control-2 control-3 P1-1 P1-2 P1-3 PA-1 PA-2 PA-3
control-1 1.0000000 0.9036573 0.9054184 0.8992221 0.9010131 0.9018963 0.9016323 0.8980050 0.8950872
control-2 0.9036573 1.0000000 0.9050045 0.8995620 0.8975293 0.9021868 0.9019837 0.8976937 0.8952984
control-3 0.9054184 0.9050045 1.0000000 0.9010410 0.9004592 0.9036801 0.9029744 0.9014315 0.8978321
P1-1 0.8992221 0.8995620 0.9010410 1.0000000 0.9005013 0.9033892 0.8974724 0.8964800 0.8937102
P1-2 0.9010131 0.8975293 0.9004592 0.9005013 1.0000000 0.9028771 0.8982444 0.8993328 0.8946727
P1-3 0.9018963 0.9021868 0.9036801 0.9033892 0.9028771 1.0000000 0.9024410 0.9003899 0.8963719
PA-1 0.9016323 0.9019837 0.9029744 0.8974724 0.8982444 0.9024410 1.0000000 0.9018240 0.8986201
PA-2 0.8980050 0.8976937 0.9014315 0.8964800 0.8993328 0.9003899 0.9018240 1.0000000 0.8965230
PA-3 0.8950872 0.8952984 0.8978321 0.8937102 0.8946727 0.8963719 0.8986201 0.8965230 1.0000000

Separate data by peptide

Also a good point to filter out any genes with low expression (average < 10 counts).

ss1 <- ss[grep("PA",rownames(ss),invert=TRUE),]
xx1 <- xx[,which(colnames(xx) %in% rownames(ss1))]
xx1 <- xx1[which(rowMeans(xx1)>10),]
dim(xx1)
## [1] 17146     6
ss2 <- ss[grep("P1",rownames(ss),invert=TRUE),]
xx2 <- xx[,which(colnames(xx) %in% rownames(ss2))]
xx2 <- xx2[which(rowMeans(xx2)>10),]
dim(xx2)
## [1] 17154     6
ss3 <- ss[grep("control",rownames(ss),invert=TRUE),]
xx3 <- xx[,which(colnames(xx) %in% rownames(ss3))]
xx3 <- xx3[which(rowMeans(xx3)>10),]
dim(xx3)
## [1] 16998     6

Analysis of differential gene expression

First we will look at control vs P1.

dds <- DESeqDataSetFromMatrix(countData = xx1 , colData = ss1, design = ~ P1 )
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
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 P1") %>% kable_paper("hover", full_width = F)
Top gene expression changes caused by P1
baseMean log2FoldChange lfcSE stat pvalue padj control-1 control-2 control-3 P1-1 P1-2 P1-3
ENSG00000124225.16 PMEPA1 8143.1899 -1.9293068 0.0844113 -22.85603 0 0 13.67197 14.00319 13.61005 12.066646 12.141998 12.130757
ENSG00000131746.13 TNS4 7155.2606 -1.0737202 0.0513338 -20.91645 0 0 13.39358 13.35150 13.44468 12.470638 12.440417 12.494875
ENSG00000140545.16 MFGE8 22565.4814 -1.2149936 0.0654488 -18.56404 0 0 14.94240 14.99461 15.04544 13.706080 13.825111 13.980465
ENSG00000140416.26 TPM1 36330.4584 -1.1591965 0.0648555 -17.87351 0 0 15.62419 15.77043 15.53542 14.437436 14.507786 14.623200
ENSG00000170525.21 PFKFB3 7666.1622 -1.0364511 0.0582596 -17.79023 0 0 13.42972 13.41801 13.57515 12.549707 12.554546 12.610695
ENSG00000146674.16 IGFBP3 3119.0789 -1.2449452 0.0704772 -17.66450 0 0 12.43271 12.44209 12.35649 11.504732 11.375675 11.552654
ENSG00000120708.17 TGFBI 39828.1867 -1.1124256 0.0638632 -17.41887 0 0 15.63721 15.86420 15.77601 14.581117 14.692049 14.767992
ENSG00000070182.21 SPTB 1900.4787 -1.4332987 0.0827732 -17.31599 0 0 11.93549 11.98739 11.74929 10.954666 10.975583 10.934915
ENSG00000113361.13 CDH6 1982.1075 -1.5690478 0.0913336 -17.17931 0 0 11.79342 12.00929 12.08089 10.863476 10.967376 10.965613
ENSG00000187634.13 SAMD11 612.1634 -1.9166579 0.1119968 -17.11351 0 0 10.89706 10.96890 10.85580 10.064108 10.031649 10.145959
ENSG00000183098.12 GPC6 1798.1506 -1.1394117 0.0677016 -16.82991 0 0 11.76954 11.78597 11.76348 10.975263 11.049420 11.058519
ENSG00000108821.14 COL1A1 712.4712 -1.9660417 0.1185062 -16.59021 0 0 10.89706 11.17299 11.02573 10.159588 10.154026 10.114317
ENSG00000134871.20 COL4A2 18327.9107 -0.8872857 0.0542645 -16.35113 0 0 14.56013 14.68468 14.56232 13.706408 13.825951 13.770841
ENSG00000112139.17 MDGA1 446.6408 -2.0201398 0.1244462 -16.23303 0 0 10.68723 10.74953 10.60370 9.916427 9.901668 9.952014
ENSG00000187498.17 COL4A1 10242.8007 -0.9983254 0.0620144 -16.09828 0 0 13.80909 13.95081 13.77814 12.863672 13.009250 12.971031
ENSG00000167244.22 IGF2 3255.7313 -1.3849107 0.0881265 -15.71504 0 0 12.53016 12.61268 12.32451 11.569911 11.376872 11.409439
ENSG00000038427.16 VCAN 8515.9821 -1.2443633 0.0820626 -15.16359 0 0 13.41886 13.77584 13.79594 12.584173 12.576560 12.576910
ENSG00000134531.10 EMP1 2691.7653 0.9608265 0.0640429 15.00285 0 0 11.50165 11.47130 11.42912 12.181601 12.218789 12.102848
ENSG00000138166.6 DUSP5 4048.8590 1.0033023 0.0674610 14.87233 0 0 11.88501 11.80702 11.89637 12.726223 12.703304 12.542086
ENSG00000137868.19 STRA6 937.8446 -1.3500772 0.0909544 -14.84345 0 0 11.17334 11.23652 11.14545 10.451194 10.431568 10.542422 
dge1 <- dge
write.table(dge1,file="deseq2_P1.tsv",quote=FALSE,sep='\t')

Now look at control vs PA.

dds <- DESeqDataSetFromMatrix(countData = xx2 , colData = ss2, design = ~ PA )
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
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 PA") %>% kable_paper("hover", full_width = F)
Top gene expression changes caused by PA
baseMean log2FoldChange lfcSE stat pvalue padj control-1 control-2 control-3 PA-1 PA-2 PA-3
ENSG00000167244.22 IGF2 3032.0533 -1.7184197 0.0846209 -20.30728 0 0 12.55700 12.63191 12.36147 11.36673 11.24199 11.33230
ENSG00000143127.13 ITGA10 1898.3449 -1.5700029 0.0782032 -20.07595 0 0 11.91085 12.02070 12.00113 11.05824 10.93114 11.01297
ENSG00000087074.9 PPP1R15A 5139.9232 1.2751594 0.0766952 16.62633 0 0 12.13083 12.05252 11.92334 12.96904 13.02449 13.17063
ENSG00000004399.13 PLXND1 2507.0654 -1.0146982 0.0634905 -15.98190 0 0 12.20537 12.15568 12.10306 11.45783 11.45689 11.45342
ENSG00000167552.16 TUBA1A 5672.8285 0.8789336 0.0566847 15.50566 0 0 12.35001 12.33919 12.35649 12.99144 13.11605 13.09966
ENSG00000162772.18 ATF3 638.8682 2.3145354 0.1498493 15.44576 0 0 10.20324 10.21657 10.05421 11.01772 10.96934 11.24261
ENSG00000070182.21 SPTB 1950.9585 -1.2370549 0.0821687 -15.05507 0 0 11.98132 12.02647 11.80609 11.14198 11.14926 11.16606
ENSG00000187720.15 THSD4 2597.2011 -1.0364188 0.0704640 -14.70849 0 0 12.25432 12.23288 12.10589 11.46236 11.45149 11.52119
ENSG00000131746.13 TNS4 7385.0358 -0.8871514 0.0629630 -14.09005 0 0 13.40190 13.35465 13.45698 12.55558 12.73888 12.64285
ENSG00000185477.6 GPRIN3 1428.0135 -1.1981448 0.0864362 -13.86160 0 0 11.66216 11.68705 11.51524 10.94005 10.88220 10.96385
ENSG00000103257.9 SLC7A5 36070.2863 0.8789024 0.0637873 13.77863 0 0 14.74286 14.68552 14.68372 15.41621 15.69001 15.53476
ENSG00000125148.7 MT2A 8229.9034 0.7384889 0.0541415 13.63998 0 0 12.82990 12.85523 12.88085 13.44090 13.56489 13.48007
ENSG00000128510.12 CPA4 419.7991 1.6999704 0.1289706 13.18107 0 0 10.19951 10.10097 10.08600 10.70041 10.73495 10.73755
ENSG00000001617.12 SEMA3F 786.3019 -1.1816806 0.0905709 -13.04703 0 0 11.14901 11.08756 11.07066 10.55577 10.55017 10.54407
ENSG00000089127.15 OAS1 2145.3616 -1.0724892 0.0827958 -12.95342 0 0 12.00097 11.87058 12.12585 11.29662 11.29937 11.29290
ENSG00000172572.7 PDE3A 1161.7406 -1.0105152 0.0782886 -12.90757 0 0 11.36947 11.42482 11.39492 10.85772 10.81285 10.86800
ENSG00000130433.8 CACNG6 2039.7500 -0.9378125 0.0732255 -12.80719 0 0 11.99714 11.90369 11.85111 11.31832 11.26538 11.33342
ENSG00000047457.14 CP 1503.1175 -1.0983780 0.0860599 -12.76295 0 0 11.70899 11.53680 11.70241 11.02979 11.02324 10.95239
ENSG00000079257.8 LXN 2579.2827 -0.9644283 0.0776956 -12.41290 0 0 12.28175 12.18264 12.04616 11.50956 11.47615 11.52216
ENSG00000079215.15 SLC1A3 2756.5387 -0.8400147 0.0680309 -12.34754 0 0 12.22340 12.14282 12.26988 11.66908 11.56941 11.61172 
dge2 <- dge
write.table(dge2,file="deseq2_PA.tsv",quote=FALSE,sep='\t')

Now look at P1 (ctrl) vs PA (trt).

dds <- DESeqDataSetFromMatrix(countData = xx3 , colData = ss3, design = ~ PA )
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
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 differences between P1 (ctrl) and PA (trt)") %>% kable_paper("hover", full_width = F)
Top gene expression differences between P1 (ctrl) and PA (trt)
baseMean log2FoldChange lfcSE stat pvalue padj P1-1 P1-2 P1-3 PA-1 PA-2 PA-3
ENSG00000124225.16 PMEPA1 6774.771 1.5691693 0.0438888 35.75328 0 0 12.23735 12.30871 12.29498 13.53543 13.52634 13.51016
ENSG00000130635.17 COL5A1 7069.383 1.9233954 0.0552575 34.80788 0 0 12.06545 12.11753 12.21366 13.58751 13.68312 13.67272
ENSG00000163735.7 CXCL5 54013.886 -1.1474498 0.0337063 -34.04264 0 0 16.18697 16.23300 16.22789 15.12502 15.09422 15.10097
ENSG00000113140.11 SPARC 13422.585 1.7847703 0.0562106 31.75152 0 0 12.81560 12.89174 12.99054 14.36164 14.51851 14.48467
ENSG00000108602.18 ALDH3A1 17844.694 -1.1449323 0.0411265 -27.83927 0 0 14.69132 14.72351 14.62453 13.60229 13.65910 13.65673
ENSG00000120708.17 TGFBI 42412.671 1.2260775 0.0464459 26.39798 0 0 14.62494 14.73720 14.80704 15.86822 15.92350 15.90576
ENSG00000213626.13 LBH 1295.755 2.5036710 0.0973782 25.71079 0 0 10.63946 10.64796 10.78013 11.81517 11.93653 11.82233
ENSG00000264230.9 ANXA8L1 2383.041 1.7634777 0.0693289 25.43641 0 0 11.20189 11.30226 11.35229 12.39147 12.33535 12.34758
ENSG00000125730.18 C3 1462.155 1.7126909 0.0680432 25.17065 0 0 10.96193 10.98166 11.02071 11.85888 11.88094 11.88293
ENSG00000113361.13 CDH6 2128.930 1.6889951 0.0686197 24.61386 0 0 11.17643 11.26738 11.26349 12.24841 12.16222 12.27820
ENSG00000146477.6 SLC22A3 1354.452 1.8520193 0.0753666 24.57349 0 0 10.85212 10.91176 10.93975 11.84299 11.84600 11.78513
ENSG00000140416.26 TPM1 35837.808 1.1130550 0.0468099 23.77819 0 0 14.48455 14.55689 14.66519 15.66089 15.62917 15.59918
ENSG00000168056.18 LTBP3 6754.093 1.0085822 0.0441605 22.83900 0 0 12.57072 12.55853 12.57422 13.43266 13.34696 13.35966
ENSG00000118257.17 NRP2 8105.134 0.9780253 0.0454191 21.53335 0 0 12.75510 12.84043 12.76121 13.58587 13.56802 13.63905
ENSG00000152377.15 SPOCK1 1867.992 1.6178651 0.0751857 21.51825 0 0 11.20189 11.22553 11.09091 12.12265 12.03210 12.11044
ENSG00000115414.21 FN1 13903.719 1.1516221 0.0539833 21.33294 0 0 13.24216 13.32970 13.40083 14.34777 14.26214 14.42833
ENSG00000166923.12 GREM1 8943.822 1.0838568 0.0517023 20.96342 0 0 12.74425 12.88022 12.89585 13.77384 13.77507 13.71006
ENSG00000095752.7 IL11 4157.910 1.3150971 0.0641814 20.49030 0 0 11.83690 11.90921 12.03719 12.87896 12.89514 12.86892
ENSG00000277443.3 MARCKS 5841.722 0.9801698 0.0481513 20.35604 0 0 12.36832 12.46830 12.44108 13.19434 13.16545 13.22684
ENSG00000089127.15 OAS1 2278.760 -1.1654483 0.0573409 -20.32491 0 0 12.24488 12.18477 12.16865 11.48473 11.48004 11.47657 
dge3 <- dge
write.table(dge3,file="deseq3_P1vsPA.tsv",quote=FALSE,sep='\t')

Make some plots.

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(7,15), cexRow=0.9, cexCol=1.4,
    main=paste("Top ranked",n,"genes in",name) )
}

maplot(dge1,"ctrl vs P1")

make_volcano(dge1,"ctrl vs P1")

make_heatmap(dge1,"ctrl vs P1",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

maplot(dge2,"ctrl vs PA")

make_volcano(dge2,"ctrl vs PA")

make_heatmap(dge2,"ctrl vs PA",ss2,xx2,n=30)
## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf
## Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

maplot(dge3,"P1 vs PA")

make_volcano(dge3,"P1 vs PA")

make_heatmap(dge3,"P1 vs PA",ss3,xx3,n=30)
## Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
## Inf
## Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
## -Inf

Pathway enrichment

Here I’m using the mitch package and gene pathways from Reactome to understand the affected pathways separately for each tissue type.

Gene ontology terms obtained from GO website and converted to list format, downloaded in Feb 2024 (GO_bp_2023q4.Rds).

First look at P1.

go <- readRDS("GO_bp_2023q4.Rds")
go <- go[which(lapply(go,length)>4)]
gobp <- go[grep(" BP ",names(go))]
length(gobp)
## [1] 5130
gt <- as.data.frame(rownames(xx))
gt$gene <- sapply(strsplit(gt[,1]," "),"[[",2)
head(gt)
##                rownames(xx)   gene
## 1 ENSG00000000003.16 TSPAN6 TSPAN6
## 2    ENSG00000000005.6 TNMD   TNMD
## 3   ENSG00000000419.14 DPM1   DPM1
## 4  ENSG00000000457.14 SCYL3  SCYL3
## 5  ENSG00000000460.17 FIRRM  FIRRM
## 6    ENSG00000000938.13 FGR    FGR
m1 <- mitch_import(dge1, 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 = 17146
## Note: no. genes in output = 17082
## Note: estimated proportion of input genes in output = 0.996
mres1 <- mitch_calc(m1, gobp, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(mres1$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences caused by P1") %>%
  kable_paper("hover", full_width = F)
Top gene pathway differences caused by P1
set setSize pANOVA s.dist p.adjustANOVA
1828 GO:0006122 BP mitochondrial electron transport, ubiquinol to cytochrome c 12 4.00e-07 0.8480570 0.0000259
1758 GO:0006123 BP mitochondrial electron transport, cytochrome c to oxygen 16 0.00e+00 0.8123022 0.0000019
1596 GO:1904851 BP positive regulation of establishment of protein localization to telomere 10 1.69e-05 0.7855787 0.0007655
1195 GO:0042776 BP proton motive force-driven mitochondrial ATP synthesis 54 0.00e+00 0.7844096 0.0000000
1911 GO:0000338 BP protein deneddylation 11 9.90e-06 0.7694974 0.0004981
1486 GO:0032543 BP mitochondrial translation 94 0.00e+00 0.7671610 0.0000000
1882 GO:0006120 BP mitochondrial electron transport, NADH to ubiquinone 38 0.00e+00 0.7639824 0.0000000
2210 GO:0030150 BP protein import into mitochondrial matrix 18 0.00e+00 0.7614406 0.0000022
1903 GO:0002181 BP cytoplasmic translation 89 0.00e+00 0.7523475 0.0000000
1598 GO:1904874 BP positive regulation of telomerase RNA localization to Cajal body 15 6.00e-07 0.7449425 0.0000405
1651 GO:0044331 BP cell-cell adhesion mediated by cadherin 13 5.10e-06 -0.7304876 0.0002895
2165 GO:0000463 BP maturation of LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 15 1.10e-06 0.7274038 0.0000718
1999 GO:0044571 BP [2Fe-2S] cluster assembly 11 2.98e-05 0.7268414 0.0011823
1597 GO:1904871 BP positive regulation of protein localization to Cajal body 11 3.30e-05 0.7228154 0.0012419
2077 GO:0000470 BP maturation of LSU-rRNA 12 1.59e-05 0.7193907 0.0007359
1017 GO:0003094 BP glomerular filtration 11 8.82e-05 -0.6826623 0.0026463
1759 GO:0045333 BP cellular respiration 34 0.00e+00 0.6788652 0.0000000
1194 GO:0015986 BP proton motive force-driven ATP synthesis 20 2.00e-07 0.6730981 0.0000159
593 GO:0009060 BP aerobic respiration 60 0.00e+00 0.6650589 0.0000000
1346 GO:0032981 BP mitochondrial respiratory chain complex I assembly 55 0.00e+00 0.6649856 0.0000000 
write.table(mres1$enrichment_result,file="mitch_P1.tsv",quote=FALSE,sep='\t')

par(mar=c(5,27,3,3))

top <- mres1$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 243
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.65,main="Ctrl vs P1",xlab="ES")
grid()

Now look at PA.

m2 <- mitch_import(dge2, 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 = 17154
## Note: no. genes in output = 17089
## Note: estimated proportion of input genes in output = 0.996
mres2 <- mitch_calc(m2, gobp, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(mres2$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences caused by PA") %>%
  kable_paper("hover", full_width = F)
Top gene pathway differences caused by PA
set setSize pANOVA s.dist p.adjustANOVA
1591 GO:1904851 BP positive regulation of establishment of protein localization to telomere 10 0.0000042 0.8404122 0.0001677
1717 GO:0048026 BP positive regulation of mRNA splicing, via spliceosome 19 0.0000000 0.8250794 0.0000001
1593 GO:1904874 BP positive regulation of telomerase RNA localization to Cajal body 15 0.0000000 0.8236227 0.0000026
2071 GO:0000470 BP maturation of LSU-rRNA 12 0.0000025 0.7848763 0.0001173
1864 GO:2000767 BP positive regulation of cytoplasmic translation 13 0.0000026 0.7524551 0.0001184
2159 GO:0000463 BP maturation of LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 15 0.0000008 0.7350592 0.0000505
1191 GO:0015986 BP proton motive force-driven ATP synthesis 20 0.0000000 0.7326205 0.0000013
1732 GO:0000387 BP spliceosomal snRNP assembly 28 0.0000000 0.7290353 0.0000000
2109 GO:0030490 BP maturation of SSU-rRNA 20 0.0000000 0.7261937 0.0000016
2160 GO:0000447 BP endonucleolytic cleavage in ITS1 to separate SSU-rRNA from 5.8S rRNA and LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 10 0.0000763 0.7223140 0.0019188
1463 GO:0000462 BP maturation of SSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 24 0.0000000 0.7206563 0.0000001
395 GO:0050850 BP positive regulation of calcium-mediated signaling 11 0.0000384 -0.7167754 0.0011328
1955 GO:0042274 BP ribosomal small subunit biogenesis 71 0.0000000 0.6881463 0.0000000
2085 GO:0000727 BP double-strand break repair via break-induced replication 12 0.0000395 0.6851418 0.0011489
2173 GO:1903241 BP U2-type prespliceosome assembly 25 0.0000000 0.6803094 0.0000004
2052 GO:0045040 BP protein insertion into mitochondrial outer membrane 14 0.0000142 0.6699770 0.0004684
850 GO:0009303 BP rRNA transcription 13 0.0000295 0.6690091 0.0009074
1584 GO:0006270 BP DNA replication initiation 25 0.0000000 0.6673136 0.0000007
1592 GO:1904871 BP positive regulation of protein localization to Cajal body 11 0.0001608 0.6570388 0.0035243
1482 GO:0032543 BP mitochondrial translation 94 0.0000000 0.6466545 0.0000000 
write.table(mres2$enrichment_result,file="mitch_PA.tsv",quote=FALSE,sep='\t')

par(mar=c(5,27,3,3))

top <- mres2$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 253
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.65,main="Ctrl vs PA",xlab="ES")
grid()

Now compare P1 (ctrl) vs PA (trt).

m3 <- mitch_import(dge3, 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 = 16998
## Note: no. genes in output = 16940
## Note: estimated proportion of input genes in output = 0.997
mres3 <- mitch_calc(m3, gobp, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(mres3$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences between P1 (ctrl) vs PA (trt)") %>%
  kable_paper("hover", full_width = F)
Top gene pathway differences between P1 (ctrl) vs PA (trt)
set setSize pANOVA s.dist p.adjustANOVA
262 GO:0045176 BP apical protein localization 10 0.0000082 0.8145304 0.0014961
259 GO:0034333 BP adherens junction assembly 11 0.0000646 0.6955842 0.0059286
1708 GO:0048026 BP positive regulation of mRNA splicing, via spliceosome 19 0.0000006 0.6627672 0.0001780
2147 GO:0000447 BP endonucleolytic cleavage in ITS1 to separate SSU-rRNA from 5.8S rRNA and LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 10 0.0009884 0.6014885 0.0453236
614 GO:0010001 BP glial cell differentiation 14 0.0002902 0.5593423 0.0199607
2073 GO:0000727 BP double-strand break repair via break-induced replication 12 0.0012200 0.5391757 0.0516377
1646 GO:0050807 BP regulation of synapse organization 10 0.0032226 0.5379327 0.0876798
1438 GO:0061158 BP 3’-UTR-mediated mRNA destabilization 16 0.0001966 0.5376093 0.0149184
2029 GO:0006405 BP RNA export from nucleus 17 0.0001284 0.5364401 0.0100905
577 GO:0050873 BP brown fat cell differentiation 24 0.0000072 -0.5291785 0.0014382
1843 GO:0000380 BP alternative mRNA splicing, via spliceosome 22 0.0000261 0.5178992 0.0040958
841 GO:0048711 BP positive regulation of astrocyte differentiation 10 0.0050125 0.5124749 0.1060813
1143 GO:0098761 BP cellular response to interleukin-7 11 0.0040176 0.5009210 0.0990624
548 GO:0044030 BP regulation of DNA methylation 11 0.0040507 0.5004699 0.0990624
1379 GO:0043068 BP positive regulation of programmed cell death 11 0.0046825 -0.4924471 0.1044054
2108 GO:0019985 BP translesion synthesis 13 0.0021469 0.4916088 0.0750059
1456 GO:0000462 BP maturation of SSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 24 0.0000340 0.4887680 0.0043969
2150 GO:0006376 BP mRNA splice site recognition 16 0.0007617 0.4860775 0.0399142
1639 GO:0044331 BP cell-cell adhesion mediated by cadherin 13 0.0026337 0.4817474 0.0816461
2197 GO:0036159 BP inner dynein arm assembly 11 0.0057585 0.4808049 0.1159190 
write.table(mres3$enrichment_result,file="mitch_P1vsPA.tsv",quote=FALSE,sep='\t')

par(mar=c(5,27,3,3))

top <- mres3$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 50
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.65,main="P1 vs PA",xlab="ES")
grid()

2D Pathway enrichment

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

l <- list("CTL_v_P1"=dge1,"CTL_v_PA"=dge2)
mm <- mitch_import(l, DEtype="deseq2",geneTable=gt)
## Note: Mean no. genes in input = 17150
## Note: no. genes in output = 16825
## Note: estimated proportion of input genes in output = 0.981
head(mm)
##             CTL_v_P1   CTL_v_PA
## 5_8S_rRNA  5.2119983  4.4662250
## A1BG       0.5419210 -0.2639765
## A1BG-AS1   0.2337854  0.9360647
## A1CF      -1.5885330 -1.1467790
## A2M-AS1   -0.4343253  0.3703636
## A4GALT    -1.7775305 -2.2559323
mmres <- mitch_calc(mm, gobp, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
head(mmres$enrichment_result,20) %>%
  kbl(caption = "Top pathways for P1 and PA vs control") %>%
  kable_paper("hover", full_width = F)
Top pathways for P1 and PA vs control
set setSize pMANOVA s.CTL_v_P1 s.CTL_v_PA p.CTL_v_P1 p.CTL_v_PA s.dist SD p.adjustMANOVA
1585 GO:1904851 BP positive regulation of establishment of protein localization to telomere 10 7.00e-07 0.7836456 0.8382753 1.77e-05 0.0000044 1.1475216 0.0386291 0.0000308
1587 GO:1904874 BP positive regulation of telomerase RNA localization to Cajal body 15 0.00e+00 0.7437636 0.8221059 6.00e-07 0.0000000 1.1086219 0.0553963 0.0000002
2064 GO:0000470 BP maturation of LSU-rRNA 12 5.00e-07 0.7181050 0.7828763 1.65e-05 0.0000026 1.0623419 0.0458002 0.0000211
1816 GO:0006122 BP mitochondrial electron transport, ubiquinol to cytochrome c 12 5.00e-07 0.8465671 0.6176371 4.00e-07 0.0002114 1.0479273 0.1618780 0.0000211
1747 GO:0006123 BP mitochondrial electron transport, cytochrome c to oxygen 15 0.00e+00 0.8488519 0.5985723 0.00e+00 0.0000596 1.0386714 0.1769744 0.0000009
2152 GO:0000463 BP maturation of LSU-rRNA from tricistronic rRNA transcript (SSU-rRNA, 5.8S rRNA, LSU-rRNA) 15 0.00e+00 0.7259726 0.7338251 1.10e-06 0.0000009 1.0322478 0.0055525 0.0000020
1476 GO:0032543 BP mitochondrial translation 94 0.00e+00 0.7659269 0.6451904 0.00e+00 0.0000000 1.0014563 0.0853736 0.0000000
1185 GO:0015986 BP proton motive force-driven ATP synthesis 20 0.00e+00 0.6718239 0.7316037 2.00e-07 0.0000000 0.9932730 0.0422707 0.0000001
1186 GO:0042776 BP proton motive force-driven mitochondrial ATP synthesis 54 0.00e+00 0.7825545 0.6113551 0.00e+00 0.0000000 0.9930492 0.1210563 0.0000000
2196 GO:0030150 BP protein import into mitochondrial matrix 18 0.00e+00 0.7600405 0.6249512 0.00e+00 0.0000044 0.9839845 0.0955225 0.0000004
1891 GO:0002181 BP cytoplasmic translation 89 0.00e+00 0.7510769 0.6239594 0.00e+00 0.0000000 0.9764434 0.0898856 0.0000000
1586 GO:1904871 BP positive regulation of protein localization to Cajal body 11 1.28e-05 0.7213145 0.6567363 3.43e-05 0.0001620 0.9754984 0.0456637 0.0003633
1726 GO:0000387 BP spliceosomal snRNP assembly 28 0.00e+00 0.6255837 0.7277320 0.00e+00 0.0000000 0.9596608 0.0722298 0.0000000
1899 GO:0000338 BP protein deneddylation 11 1.87e-05 0.7670880 0.5505044 1.05e-05 0.0015690 0.9441818 0.1531477 0.0004625
2102 GO:0030490 BP maturation of SSU-rRNA 20 0.00e+00 0.5948468 0.7247307 4.10e-06 0.0000000 0.9375912 0.0918418 0.0000004
2045 GO:0045040 BP protein insertion into mitochondrial outer membrane 14 2.90e-06 0.6316697 0.6692727 4.26e-05 0.0000145 0.9202894 0.0265893 0.0001013
1986 GO:0044571 BP [2Fe-2S] cluster assembly 11 3.83e-05 0.7248938 0.5650486 3.13e-05 0.0011738 0.9191033 0.1130276 0.0008796
1870 GO:0006120 BP mitochondrial electron transport, NADH to ubiquinone 38 0.00e+00 0.7621186 0.4989230 0.00e+00 0.0000001 0.9109055 0.1861073 0.0000000
1949 GO:0042274 BP ribosomal small subunit biogenesis 71 0.00e+00 0.5577798 0.6871548 0.00e+00 0.0000000 0.8850424 0.0914820 0.0000000
2020 GO:0042273 BP ribosomal large subunit biogenesis 35 0.00e+00 0.5824794 0.6440943 0.00e+00 0.0000000 0.8684121 0.0435683 0.0000000 
write.table(mmres$enrichment_result,file="mitch_multi.tsv",quote=FALSE,sep='\t')

Generate detailed HTML reports

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

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

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

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

Venn diagrams

Gene level.

dge1_up <- rownames(subset(dge1,padj<0.05 & log2FoldChange>0))
dge1_dn <- rownames(subset(dge1,padj<0.05 & log2FoldChange<0))

dge2_up <- rownames(subset(dge2,padj<0.05 & log2FoldChange>0))
dge2_dn <- rownames(subset(dge2,padj<0.05 & log2FoldChange<0))

v1 <- list("P1 up"=dge1_up,"P1 dn"=dge1_dn,"PA up"=dge2_up,"PA dn"=dge2_dn)

plot(euler(v1),quantities = TRUE, main="Gene level")

Pathway level.

m1up <- subset(mres1$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
m1dn <- subset(mres1$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

m2up <- subset(mres2$enrichment_result,p.adjustANOVA < 0.05 & s.dist > 0)$set
m2dn <- subset(mres2$enrichment_result,p.adjustANOVA < 0.05 & s.dist < 0)$set

v1 <- list("P1 up"=m1up,"P1 dn"=m1dn,"PA up"=m2up,"PA dn"=m2dn)

plot(euler(v1),quantities = TRUE, main="Pathway level")

Session information

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_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] kableExtra_1.4.0            eulerr_7.0.2               
##  [3] mitch_1.16.0                MASS_7.3-61                
##  [5] gplots_3.1.3.1              DESeq2_1.44.0              
##  [7] SummarizedExperiment_1.34.0 Biobase_2.64.0             
##  [9] MatrixGenerics_1.16.0       matrixStats_1.3.0          
## [11] GenomicRanges_1.56.0        GenomeInfoDb_1.40.0        
## [13] IRanges_2.38.0              S4Vectors_0.42.0           
## [15] BiocGenerics_0.50.0         reshape2_1.4.4             
## [17] dplyr_1.1.4                 zoo_1.8-12                 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        viridisLite_0.4.2       bitops_1.0-7           
##  [4] fastmap_1.2.0           GGally_2.2.1            promises_1.3.0         
##  [7] digest_0.6.36           mime_0.12               lifecycle_1.0.4        
## [10] polylabelr_0.2.0        magrittr_2.0.3          compiler_4.4.1         
## [13] sass_0.4.9              rlang_1.1.4             tools_4.4.1            
## [16] yaml_2.3.9              utf8_1.2.4              knitr_1.48             
## [19] S4Arrays_1.4.0          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-5             BiocParallel_1.38.0     KernSmooth_2.23-24     
## [28] purrr_1.0.2             polyclip_1.10-6         grid_4.4.1             
## [31] fansi_1.0.6             caTools_1.18.2          xtable_1.8-4           
## [34] colorspace_2.1-0        ggplot2_3.5.1           scales_1.3.0           
## [37] gtools_3.9.5            cli_3.6.3               rmarkdown_2.27         
## [40] crayon_1.5.3            generics_0.1.3          rstudioapi_0.16.0      
## [43] httr_1.4.7              cachem_1.1.0            stringr_1.5.1          
## [46] zlibbioc_1.50.0         parallel_4.4.1          XVector_0.44.0         
## [49] vctrs_0.6.5             Matrix_1.7-0            jsonlite_1.8.8         
## [52] echarts4r_0.4.5         beeswarm_0.4.0          systemfonts_1.1.0      
## [55] locfit_1.5-9.10         jquerylib_0.1.4         tidyr_1.3.1            
## [58] glue_1.7.0              ggstats_0.6.0           codetools_0.2-20       
## [61] stringi_1.8.4           gtable_0.3.5            later_1.3.2            
## [64] UCSC.utils_1.0.0        munsell_0.5.1           tibble_3.2.1           
## [67] pillar_1.9.0            htmltools_0.5.8.1       GenomeInfoDbData_1.2.12
## [70] R6_2.5.1                shiny_1.8.1.1           evaluate_0.24.0        
## [73] lattice_0.22-6          highr_0.11              bslib_0.7.0            
## [76] httpuv_1.6.15           Rcpp_1.0.12             svglite_2.1.3          
## [79] gridExtra_2.3           SparseArray_1.4.3       xfun_0.45              
## [82] pkgconfig_2.0.3