Source codes: https://github.com/markziemann/tohidul_rnaseq

Background

Here we have n=3 control (H2O; “H”) and n=3 seaweed based fertiliser treatments:

Reads underwent quality trimming using Skewer (Jiang et al, 2014). I mapped the reads to the Arabidopsis transcriptome (TAIR10/Ensembl47) using Kallisto (Bray et al, 2016). Expression counts were loaded into R and then DE analysis was performed with DESeq2 (Love et al, 2014). Enrichment analysis was performed using Plant Reactome genesets with the Mitch package (Kaspi & Ziemann 2020).

suppressPackageStartupMessages({
  library("reshape2")
  library("gplots")
  library("DESeq2")
  library("mitch")
  library("kableExtra")
  library("eulerr")
  library("UpSetR")

})

Load data

Here we load the data in from the aligner. Take note that the columns are arranged in order.

tmp <- read.table("3col.tsv.gz")
x <- as.data.frame(acast(tmp, V2~V1, value.var="V3"))
x$gene <- sapply(strsplit(rownames(x),"\\."),"[[",1)
xx <- aggregate(. ~ gene, x, sum)
rownames(xx) <- xx$gene
xx$gene = NULL
xx <- xx[,order(colnames(xx))]

Sample sheet

The sample sheet is read in and put in order as well.

ss <- read.table("samplesheet.tsv",header=TRUE)
ss <- ss[order(ss$Run),]
ss$label <- gsub("-0hpi","",ss$SampleName)
# check that names are in order
colnames(xx) == ss$Run
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# change data header
colnames(xx) <- ss$label
kbl(ss) %>% kable_styling()
Run SampleName Treatment label
6 SRR11404271 S80-0hpi_R3 Seaweed_16280 S80_R3
5 SRR11404272 S80-0hpi_R2 Seaweed_16280 S80_R2
4 SRR11404273 S80-0hpi_R1 Seaweed_16280 S80_R1
12 SRR11404287 S94-0hpi_R3 Seaweed_15294 S94_R3
11 SRR11404288 S94-0hpi_R2 Seaweed_15294 S94_R2
10 SRR11404289 S94-0hpi_R1 Seaweed_15294 S94_R1
3 SRR11404301 H-0hpi_R3 Water H_R3
9 SRR11404304 S93-0hpi_R3 Seaweed_15293 S93_R3
8 SRR11404305 S93-0hpi_R2 Seaweed_15293 S93_R2
7 SRR11404306 S93-0hpi_R1 Seaweed_15293 S93_R1
2 SRR11404312 H-0hpi_R2 Water H_R2
1 SRR11404313 H-0hpi_R1 Water H_R1

MDS

MDS is just like PCA. The more similar (correlated) the data sets are the closer they will appear on the scatterplot.

cols <- as.numeric(factor(ss$Treatment))
plot(cmdscale(dist(t(xx))), xlab="Coordinate 1", ylab="Coordinate 2", 
  col=cols, cex=4 , main="MDS plot")
text(cmdscale(dist(t(xx))), labels=ss$label )

#colfunc <- colorRampPalette(c("white", "yellow", "orange" ,  "red", "darkred"))
colfunc <- colorRampPalette(c("white","blue"))

heatmap.2(cor(xx,method="pearson"),col=colfunc(25),trace="none",
  scale="none",margin=c(10,10),main="Pearson correlation")

heatmap.2(cor(xx,method="spearman"),col=colfunc(25),trace="none",
  scale="none",margin=c(10,10),main="Spearman correlation")

Functions

run_de <- function(ss,xx){
xx <- xx[which(rowMeans(xx)>10),]
y <- round(xx)
# MDS
mds <- cmdscale(dist(t(y)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds , xlab="Coordinate 1", ylab="Coordinate 2",
  type = "n" , xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(y) )
# DE
dds <- DESeqDataSetFromMatrix(countData=y, colData = ss, design = ~ trt)
dds <- DESeq(dds)
de <- DESeq2::results(dds)
de <- de[order(de$pvalue),]
up <- rownames(subset(de, log2FoldChange>0 & padj<0.05 ))
dn <- rownames(subset(de, log2FoldChange<0 & padj<0.05 ))
str(up)
str(dn)
# MA plot
sig <-subset(de, padj < 0.05 )
GENESUP <- length(up)
GENESDN <- length(dn)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down")
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="smear plot")
points(log2(sig$baseMean),sig$log2FoldChange,
       pch=19, cex=0.5, col="red")
mtext(SUBHEADER)
# heatmap
yn <- y/colSums(y)*1000000
yf <- yn[which(rownames(yn) %in% rownames(de)[1:50]),]
mycols <- gsub("0","yellow",ss$trt)
mycols <- gsub("1","orange",mycols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(  as.matrix(yf), col=colfunc(25),scale="row",
    ColSideColors =mycols ,trace="none",
    margin = c(10,10), cexRow=0.6, cexCol=0.8 , main="Top 50 genes by p-val")
mtext("yellow=ctrl, orange=trt")
return(de)
}


mitch_barplot <- function(res){
  sig <- res$enrichment_result
  sig <- head( sig[order(sig$pANOVA),] ,30)
  sig <- sig[order(sig$s.dist),]
  par(mar=c(3,25,1,1)); barplot(sig$s.dist,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6,
    names.arg=sig$set,main="Enrichment score") ;grid()
}

Split data into contrast groups

ss93 <- subset(ss,Treatment=="Seaweed_15293"|Treatment=="Water")
ss93$trt <- as.numeric(grepl("Seaweed",ss93$Treatment))
xx93 <- xx[,which(colnames(xx) %in% ss93$label)]

ss94 <- subset(ss,Treatment=="Seaweed_15294"|Treatment=="Water")
ss94$trt <- as.numeric(grepl("Seaweed",ss94$Treatment))
xx94 <- xx[,which(colnames(xx) %in% ss94$label)]

ss80 <- subset(ss,Treatment=="Seaweed_16280"|Treatment=="Water")
ss80$trt <- as.numeric(grepl("Seaweed",ss80$Treatment))
xx80 <- xx[,which(colnames(xx) %in% ss80$label)]

DE

Here, were using DESeq2 to perform differential expression analysis to understand gene expression changes caused by the different formulations. Enrichment analysis is performed at the end using mitch with all 4 timepoints.

de93 <- run_de(ss93,xx93)
## 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
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

##  chr [1:131] "ATCG00650" "ATCG00270" "ATCG00450" "ATCG00090" "AT2G29380" ...
##  chr [1:74] "AT1G07550" "AT4G15393" "AT1G61480" "AT4G13420" "AT5G42600" ...

de94 <- run_de(ss94,xx94)
## 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
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

##  chr [1:48] "AT1G05523" "AT1G07893" "AT5G47450" "AT2G29380" "ATCG01000" ...
##  chr [1:63] "AT3G20380" "AT4G13420" "AT2G46390" "AT1G67148" "AT5G28630" ...

de80 <- run_de(ss80,xx80)
## 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
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

##  chr [1:729] "AT1G65970" "AT4G06665" "AT3G11510" "AT5G64100" "AT1G24996" ...
##  chr [1:730] "AT2G47450" "AT1G55673" "AT2G22510" "AT2G32870" "AT2G18328" ...

de93_top <- head(as.data.frame(de93),40)
kbl(de93_top[order(-de93_top$log2FoldChange),]) %>% kable_styling()
baseMean log2FoldChange lfcSE stat pvalue padj
AT2G05250 20.06891 7.6582985 1.3561246 5.647194 0e+00 0.0000138
ATCG00030 156.69794 2.4690560 0.4318902 5.716861 0e+00 0.0000108
ATCG00460 81.51351 2.2885593 0.3389210 6.752486 0e+00 0.0000000
ATCG00150 613.03109 2.2836146 0.3304977 6.909622 0e+00 0.0000000
AT1G08440 148.04467 2.2193028 0.3056109 7.261857 0e+00 0.0000000
AT2G29440 66.63172 2.1181942 0.3736903 5.668314 0e+00 0.0000127
ATCG00300 475.91230 2.0292842 0.3685969 5.505430 0e+00 0.0000289
ATCG00270 5911.36124 1.9394232 0.2467682 7.859292 0e+00 0.0000000
ATCG00650 399.31947 1.8721366 0.1957331 9.564740 0e+00 0.0000000
AT2G29380 103.12033 1.8664464 0.2549292 7.321432 0e+00 0.0000000
ATCG00160 207.44469 1.8551026 0.3500219 5.299961 1e-07 0.0000747
ATCG00450 107.87068 1.8385628 0.2376897 7.735140 0e+00 0.0000000
AT5G05250 194.13408 1.7381615 0.2879731 6.035847 0e+00 0.0000019
ATCG00550 194.49345 1.6110753 0.2224401 7.242739 0e+00 0.0000000
ATCG00090 271.96536 1.3955804 0.1811235 7.705130 0e+00 0.0000000
ATCG00440 3614.44781 1.3786288 0.2343762 5.882121 0e+00 0.0000042
AT2G38600 313.27568 1.3159796 0.1959970 6.714286 0e+00 0.0000000
AT1G43800 226.58505 1.2514585 0.2083557 6.006355 0e+00 0.0000022
ATCG00080 244.13459 1.1923254 0.2025920 5.885352 0e+00 0.0000042
AT3G12500 683.07210 0.9640148 0.1390948 6.930631 0e+00 0.0000000
AT5G41080 753.96991 0.8650560 0.1587080 5.450615 1e-07 0.0000355
ATCG00820 1138.90003 0.8390081 0.1649416 5.086699 4e-07 0.0002103
ATCG00830 3961.83710 0.7438248 0.1206828 6.163472 0e+00 0.0000010
ATCG01310 3961.83710 0.7438248 0.1206828 6.163472 0e+00 0.0000010
AT1G25220 1104.70865 0.6507050 0.1296480 5.019012 5e-07 0.0002921
AT2G30860 19565.46921 0.6150604 0.1134059 5.423529 1e-07 0.0000401
AT2G01890 931.51406 -0.6675194 0.1292639 -5.164005 2e-07 0.0001433
AT1G60095 1453.37407 -0.7102906 0.1297534 -5.474157 0e+00 0.0000321
AT1G74450 595.40288 -0.8671836 0.1665104 -5.207986 2e-07 0.0001163
AT4G15393 1989.17973 -0.9901515 0.1485213 -6.666732 0e+00 0.0000000
AT4G19975 199.66962 -0.9984882 0.1802379 -5.539835 0e+00 0.0000246
AT2G17050 256.44657 -1.0689654 0.2032928 -5.258256 1e-07 0.0000912
AT2G14510 235.54066 -1.0899585 0.2187095 -4.983591 6e-07 0.0003418
AT1G07550 358.84301 -1.1546757 0.1496251 -7.717123 0e+00 0.0000000
AT1G74810 301.28682 -1.2219915 0.2292445 -5.330516 1e-07 0.0000651
AT1G61480 144.72070 -1.3241403 0.2122536 -6.238483 0e+00 0.0000007
AT5G08275 80.48553 -1.4168300 0.2585898 -5.479064 0e+00 0.0000321
AT5G42600 192.60018 -1.5297996 0.2693292 -5.680037 0e+00 0.0000127
AT4G13420 4964.41642 -1.9741099 0.3251097 -6.072134 0e+00 0.0000016
AT4G01985 583.32836 -3.1637997 0.5575694 -5.674270 0e+00 0.0000127 
de94_top <- head(as.data.frame(de94),40)
kbl(de94_top[order(-de94_top$log2FoldChange),]) %>% kable_styling()
baseMean log2FoldChange lfcSE stat pvalue padj
AT4G06665 176.949923 11.0825144 1.9998867 5.541571 0.0e+00 0.0000499
AT2G08750 80.716371 9.9494840 2.1877668 4.547781 5.4e-06 0.0032604
AT1G05523 50.263784 9.2684255 1.2794841 7.243877 0.0e+00 0.0000000
ATCG01000 34.161919 8.7077199 1.5371206 5.664955 0.0e+00 0.0000318
AT2G09975 24.725558 8.2394817 1.5397308 5.351248 1.0e-07 0.0001051
AT1G07893 291.385881 3.0347990 0.4746791 6.393370 0.0e+00 0.0000009
AT2G29380 109.958811 2.3300527 0.3813925 6.109330 0.0e+00 0.0000031
AT1G08440 85.917146 1.5807682 0.3547177 4.456412 8.3e-06 0.0045109
AT5G47450 2398.597002 1.3591791 0.2170536 6.261951 0.0e+00 0.0000014
AT1G57590 504.091543 1.1338581 0.2509785 4.517749 6.3e-06 0.0035609
AT4G22530 356.754173 0.9359380 0.1871747 5.000345 6.0e-07 0.0005387
AT3G57010 687.200399 0.8450270 0.1742118 4.850574 1.2e-06 0.0010251
AT4G11570 2785.154584 0.6875898 0.1498746 4.587768 4.5e-06 0.0028586
AT1G35580 5760.956899 0.6528030 0.1340270 4.870682 1.1e-06 0.0009631
AT1G06430 3533.787890 0.6351913 0.1402416 4.529264 5.9e-06 0.0034634
AT4G10960 1006.088382 0.6234614 0.1334511 4.671835 3.0e-06 0.0020848
AT5G58200 1709.302202 -0.6244188 0.1392258 -4.484937 7.3e-06 0.0040489
AT2G40765 1293.208229 -0.7169874 0.1407630 -5.093577 4.0e-07 0.0003622
AT5G65040 903.352822 -0.8117109 0.1473786 -5.507657 0.0e+00 0.0000562
AT2G46390 970.073218 -0.9195554 0.1451276 -6.336184 0.0e+00 0.0000010
AT1G54950 280.407408 -0.9595267 0.2047506 -4.686320 2.8e-06 0.0020074
AT3G62040 730.548089 -0.9724651 0.1865740 -5.212222 2.0e-07 0.0002126
AT5G65980 480.343952 -1.0101416 0.1812039 -5.574614 0.0e+00 0.0000448
AT5G53880 729.219249 -1.0167014 0.2032885 -5.001273 6.0e-07 0.0005387
AT4G15393 1615.603995 -1.0464929 0.2172873 -4.816171 1.5e-06 0.0011734
AT3G52420 190.667544 -1.1073950 0.2233719 -4.957629 7.0e-07 0.0006437
AT1G49310 818.702027 -1.1100436 0.2342177 -4.739367 2.1e-06 0.0016005
AT3G52310 273.218762 -1.1626381 0.2544375 -4.569444 4.9e-06 0.0030249
AT3G48180 168.568433 -1.1761801 0.2552235 -4.608433 4.1e-06 0.0027449
AT3G07195 364.438346 -1.2081644 0.2240949 -5.391307 1.0e-07 0.0000891
AT3G20380 671.602736 -1.2220039 0.1760420 -6.941546 0.0e+00 0.0000000
AT5G52260 182.006398 -1.5566513 0.2886810 -5.392290 1.0e-07 0.0000891
AT1G67148 280.866425 -1.5632862 0.2578229 -6.063411 0.0e+00 0.0000036
AT1G71030 815.503034 -1.7013740 0.3144715 -5.410265 1.0e-07 0.0000891
AT5G28630 199.457731 -1.7434598 0.2958335 -5.893382 0.0e+00 0.0000091
AT4G13420 4050.402883 -2.0339078 0.3126679 -6.505009 0.0e+00 0.0000006
AT3G20760 85.636689 -2.9389409 0.5748671 -5.112383 3.0e-07 0.0003444
AT2G07495 9.177598 -6.4736086 1.4111899 -4.587340 4.5e-06 0.0028586
AT3G08230 11.757532 -6.8235725 1.4243698 -4.790591 1.7e-06 0.0012858
AT5G02195 17.480899 -7.3974662 1.3109886 -5.642662 0.0e+00 0.0000330 
de80_top <- head(as.data.frame(de80),40)
kbl(de80_top[order(-de80_top$log2FoldChange),]) %>% kable_styling()
baseMean log2FoldChange lfcSE stat pvalue padj
AT4G06665 132.02756 10.6703428 1.2214435 8.735846 0 0
AT2G08750 37.84520 8.8747367 1.2919905 6.869042 0 0
AT1G24996 323.82079 4.3818892 0.5618073 7.799630 0 0
AT3G20810 1371.83363 3.1462473 0.4775527 6.588273 0 0
AT5G06760 1048.96459 2.9414243 0.4414667 6.662845 0 0
AT4G33905 330.48540 2.1883867 0.3352305 6.528005 0 0
AT1G65970 319.83244 1.8879852 0.2028938 9.305287 0 0
AT3G48510 71.68292 1.7527497 0.2646986 6.621681 0 0
AT1G19250 177.80497 1.3027098 0.2005740 6.494908 0 0
AT3G14440 696.62229 1.1073515 0.1674480 6.613105 0 0
AT5G03345 831.78940 0.9323214 0.1398791 6.665195 0 0
AT5G64100 3886.71544 0.9061570 0.1091747 8.300067 0 0
AT4G12600 1344.31545 0.8738919 0.1317509 6.632912 0 0
AT3G11510 6114.57364 0.8224084 0.0964760 8.524483 0 0
AT5G42090 3764.73635 0.7636518 0.1042944 7.322079 0 0
AT1G14980 1578.18877 0.7608785 0.1126658 6.753414 0 0
AT3G22230 3405.70403 0.7473244 0.0963363 7.757453 0 0
AT3G09735 1220.96912 0.7320855 0.0981377 7.459776 0 0
AT2G44310 1560.18430 0.6872796 0.0950144 7.233422 0 0
AT3G18740 3328.73251 0.6726190 0.1004285 6.697488 0 0
AT4G09320 6585.58229 0.6659963 0.0894218 7.447803 0 0
AT4G20150 3164.43679 0.6605573 0.0950525 6.949394 0 0
AT5G02870 7649.88412 0.6494672 0.0860326 7.549084 0 0
AT1G55330 2376.62233 0.6279042 0.0917265 6.845395 0 0
AT3G53430 5234.77439 0.6126416 0.0824943 7.426474 0 0
AT5G57660 1896.01283 -0.7450316 0.1135806 -6.559499 0 0
AT2G32870 1214.28389 -0.8573352 0.1021669 -8.391513 0 0
AT3G47470 34652.01704 -0.8809525 0.1335732 -6.595277 0 0
AT3G15353 23919.22730 -0.9793861 0.1511211 -6.480802 0 0
AT2G18328 689.16212 -1.0454080 0.1455697 -7.181493 0 0
AT4G39720 211.17876 -1.0949131 0.1686647 -6.491655 0 0
AT2G47450 3594.05175 -1.1612794 0.1044117 -11.122115 0 0
AT5G48490 1040.34652 -1.2382333 0.1779417 -6.958647 0 0
AT5G54270 14956.45719 -1.2741526 0.1886138 -6.755352 0 0
AT1G64370 5582.52780 -1.2768304 0.1829606 -6.978719 0 0
AT4G38080 3041.08385 -1.4147394 0.1978268 -7.151404 0 0
AT2G22510 1187.84821 -1.5712553 0.1672951 -9.392116 0 0
AT3G02400 210.64282 -2.0573297 0.3175210 -6.479350 0 0
AT1G15830 157.00814 -2.1186271 0.3203694 -6.613077 0 0
AT1G55673 735.08530 -12.7838910 1.2309372 -10.385494 0 0 
write.table(de93,file="de93.tsv",quote=FALSE,sep="\t")
write.table(de94,file="de94.tsv",quote=FALSE,sep="\t")
write.table(de80,file="de80.tsv",quote=FALSE,sep="\t")

de93_up <- rownames(subset(de93,padj<0.05&log2FoldChange>0))
de93_dn <- rownames(subset(de93,padj<0.05&log2FoldChange<0))
de94_up <- rownames(subset(de94,padj<0.05&log2FoldChange>0))
de94_dn <- rownames(subset(de94,padj<0.05&log2FoldChange<0))
de80_up <- rownames(subset(de80,padj<0.05&log2FoldChange>0))
de80_dn <- rownames(subset(de80,padj<0.05&log2FoldChange<0))

up <- sapply(list(de93_up,de94_up,de80_up),length)
dn <- sapply(list(de93_dn,de94_dn,de80_dn),length)
res_summary <- data.frame(up,dn,row.names=c("S93","S94","S80"))
res_summary$de <- c("93","94","80")

kbl(res_summary) %>% kable_styling()
up dn de
S93 131 74 93
S94 48 63 94
S80 729 730 80 
barplot(res_summary$up,ylim=c(-750,750), axes=TRUE, main="number of DEGs up/down")

barplot(-res_summary$dn, add = TRUE, names.arg = rownames(res_summary), axes = FALSE)

Venn diagram plots of overlap

# upregulated genes
v1 <- list("93 up"=de93_up, "94 up"=de94_up, "80 up" =de80_up )
plot(euler(v1),quantities = TRUE)

intersect(de93_up,de94_up)
## [1] "AT2G29380" "AT1G08440" "AT2G29440" "AT1G59940" "AT3G57010" "AT1G17180"
## [7] "AT5G59090" "AT1G24996"
intersect(de93_up,de80_up)
##  [1] "ATCG00650" "ATCG00270" "ATCG00090" "AT2G29380" "ATCG00550" "AT3G12500"
##  [7] "ATCG00460" "ATCG00440" "ATCG00030" "AT2G29440" "ATCG00300" "AT5G41080"
## [13] "ATCG00820" "AT3G48510" "ATCG00640" "AT1G59940" "ATCG00620" "AT3G57010"
## [19] "AT1G17180" "AT4G37070" "ATCG01040" "AT3G56970" "ATCG00210" "AT1G19540"
## [25] "ATCG00360" "AT5G63350" "ATCG00670" "AT5G59090" "AT1G24996" "ATCG00040"
## [31] "AT3G56980" "AT1G61255" "ATCG00050" "AT3G15670" "AT3G48100" "ATCG00380"
## [37] "AT3G18280" "ATCG00370" "AT1G14980" "AT1G07985"
intersect(de93_up,de80_up)
##  [1] "ATCG00650" "ATCG00270" "ATCG00090" "AT2G29380" "ATCG00550" "AT3G12500"
##  [7] "ATCG00460" "ATCG00440" "ATCG00030" "AT2G29440" "ATCG00300" "AT5G41080"
## [13] "ATCG00820" "AT3G48510" "ATCG00640" "AT1G59940" "ATCG00620" "AT3G57010"
## [19] "AT1G17180" "AT4G37070" "ATCG01040" "AT3G56970" "ATCG00210" "AT1G19540"
## [25] "ATCG00360" "AT5G63350" "ATCG00670" "AT5G59090" "AT1G24996" "ATCG00040"
## [31] "AT3G56980" "AT1G61255" "ATCG00050" "AT3G15670" "AT3G48100" "ATCG00380"
## [37] "AT3G18280" "ATCG00370" "AT1G14980" "AT1G07985"
intersect(de94_up,de80_up)
##  [1] "AT1G07893" "AT2G29380" "AT4G06665" "AT4G22530" "AT1G35580" "AT3G57010"
##  [7] "AT4G10960" "AT2G08750" "AT1G57590" "AT5G64100" "AT4G02280" "AT2G36830"
## [13] "AT4G17340" "AT1G43910" "AT3G56800" "AT1G24996" "AT2G02870" "AT2G39040"
## [19] "AT1G17180" "AT4G23630" "AT1G59940" "AT3G59500" "AT2G29440" "AT5G59090"
## [25] "AT5G60790" "AT5G24080"
names(which(table(c(de93_up,de94_up,de80_up))==3))
## [1] "AT1G17180" "AT1G24996" "AT1G59940" "AT2G29380" "AT2G29440" "AT3G57010"
## [7] "AT5G59090"
# downregulated genes
v2 <- list("93 dn"=de93_dn, "94 dn"=de94_dn, "80 dn" =de80_dn )
plot(euler(v2),quantities = TRUE)

intersect(de93_dn,de94_dn)
## [1] "AT1G07550" "AT4G15393" "AT4G13420" "AT4G19975" "AT5G65980" "AT3G52310"
## [7] "AT3G10320"
intersect(de93_dn,de80_dn)
##  [1] "AT4G15393" "AT4G13420" "AT4G01985" "AT5G08275" "AT5G56270" "AT3G52310"
##  [7] "AT4G11300" "AT3G10320" "AT2G02610" "AT1G62760" "AT2G26650" "AT3G02620"
## [13] "AT3G13760"
intersect(de93_dn,de80_dn)
##  [1] "AT4G15393" "AT4G13420" "AT4G01985" "AT5G08275" "AT5G56270" "AT3G52310"
##  [7] "AT4G11300" "AT3G10320" "AT2G02610" "AT1G62760" "AT2G26650" "AT3G02620"
## [13] "AT3G13760"
intersect(de94_dn,de80_dn)
##  [1] "AT3G20380" "AT4G13420" "AT1G67148" "AT5G28630" "AT1G71030" "AT5G52260"
##  [7] "AT3G07195" "AT3G62040" "AT2G40765" "AT5G53880" "AT3G52420" "AT4G15393"
## [13] "AT1G54950" "AT3G52310" "AT5G58200" "AT5G19151" "AT3G43800" "AT5G57340"
## [19] "AT2G28110" "AT5G20110" "AT3G61210" "AT2G45660" "AT4G16240" "AT1G78780"
## [25] "AT1G70470" "AT5G48560" "AT3G10320" "AT5G57565" "AT5G44578" "AT1G56200"
## [31] "AT4G01390" "AT3G53235" "AT3G58550" "AT2G16750" "AT1G54410" "AT5G19090"
## [37] "AT5G44480" "AT3G14020" "AT1G10690"
names(which(table(c(de93_dn,de94_dn,de80_dn))==3))
## [1] "AT3G10320" "AT3G52310" "AT4G13420" "AT4G15393"
# up and downregulated genes
v3 <- list("93 up"=de93_up, "94 up"=de94_up, "80 up" =de80_up ,
  "93 dn"=de93_dn, "94 dn"=de94_dn, "80 dn" =de80_dn )
plot(euler(v3),quantities = TRUE)

intersect(de93_up,de80_dn)
## [1] "AT2G05510"

DEG summary

Here is a heatmap of all the differentially expressed genes. It looks a bit random.

degs <- unique(c(de93_up,de93_dn,de94_up,de94_dn,de80_up,de80_dn))
xxx <- xx/colSums(xx)*1e6
deg_mx <- as.matrix(xxx[which(rownames(xxx) %in% degs),])
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(deg_mx,col=colfunc(25),trace="none",scale="row",margin=c(10,10))

Pathway analysis

Mapman pathways last modified in 2012 and used in the previous RNA-seq analysis.

genesets <- gmt_import("../baseline/ref/Ath_AGI_LOCUS_TAIR10_Aug2012.txt.gmt")

gt <- read.table("../baseline/ref/Arabidopsis_thaliana.TAIR10.46.geneaccession2symbol.tsv",
    fill=TRUE) 

l <- list("S93"=de93,"S94"=de94,"S80"=de80)

m <- mitch_import(x=l,DEtype="deseq2",geneTable=gt)
## Note: Mean no. genes in input = 21817.6666666667
## Note: no. genes in output = 21283
## Note: estimated proportion of input genes in output = 0.975
# significance
res_sig <- mitch_calc(x=m,genesets=genesets,priority="significance")
## Note: When prioritising by significance (ie: small 
##             p-values), large effect sizes might be missed.
top <- head(subset(res_sig$enrichment_result,p.adjustMANOVA<0.05),50)
kbl(head(top,30)) %>% kable_styling()
set setSize pMANOVA s.S93 s.S94 s.S80 p.S93 p.S94 p.S80 s.dist SD p.adjustMANOVA
90 NOT_ASSIGNED_NO_ONTOLOGY_PENTATRICOPEPTIDE_(PPR)_REPEAT-CONTAINING_PROTEIN 417 0.0e+00 -0.3577936 0.2891597 0.1931390 0.0000000 0.0000000 0.0000000 0.4989312 0.3491168 0.00e+00
93 NOT_ASSIGNED_UNKNOWN 4521 0.0e+00 0.0411330 -0.1138630 -0.0073580 0.0000212 0.0000000 0.4469948 0.1212883 0.0792868 0.00e+00
12 CELL_WALL_CELL_WALL_PROTEINS_AGPS_AGP 39 0.0e+00 0.5658792 -0.1886530 -0.1315871 0.0000000 0.0414797 0.1550328 0.6108391 0.4201259 0.00e+00
247 TRANSPORT_ABC_TRANSPORTERS_AND_MULTIDRUG_RESISTANCE_SYSTEMS 110 0.0e+00 -0.1271843 0.2524158 -0.0458457 0.0211976 0.0000048 0.4061819 0.2863414 0.1998633 0.00e+00
111 PROTEIN_DEGRADATION_UBIQUITIN_E3_SCF_FBOX 250 0.0e+00 -0.0642514 0.1762971 0.2054725 0.0802577 0.0000016 0.0000000 0.2782586 0.1480235 0.00e+00
25 DNA_SYNTHESIS/CHROMATIN_STRUCTURE 204 0.0e+00 -0.1564962 -0.0362562 -0.2357039 0.0001165 0.3720839 0.0000000 0.2852401 0.1004248 0.00e+00
241 STRESS_BIOTIC_PR-PROTEINS 170 0.0e+00 -0.1955422 -0.0866714 -0.2698443 0.0000109 0.0512409 0.0000000 0.3443322 0.0921285 0.00e+00
91 NOT_ASSIGNED_NO_ONTOLOGY_PROLINE_RICH_FAMILY 38 0.0e+00 -0.0392835 -0.5140553 -0.5498953 0.6751904 0.0000000 0.0000000 0.7537778 0.2850196 0.00e+00
113 PROTEIN_DEGRADATION_UBIQUITIN_PROTEASOM 59 0.0e+00 0.2370549 0.2544928 0.4930459 0.0016348 0.0007213 0.0000000 0.6033705 0.1430286 0.00e+00
174 RNA_REGULATION_OF_TRANSCRIPTION_G2-LIKE_TRANSCRIPTION_FACTOR_FAMILY,_GARP 37 0.0e+00 -0.0823863 -0.5920351 -0.4031105 0.3858352 0.0000000 0.0000220 0.7209655 0.2576492 0.00e+00
85 NOT_ASSIGNED_NO_ONTOLOGY_DC1_DOMAIN_CONTAINING_PROTEIN 101 0.0e+00 -0.3728974 -0.1910010 -0.1320400 0.0000000 0.0009097 0.0218507 0.4392817 0.1255488 1.00e-07
8 CELL_ORGANISATION 344 0.0e+00 -0.0593010 -0.0476817 -0.1865938 0.0588195 0.1286907 0.0000000 0.2015128 0.0770660 1.00e-07
144 PS_LIGHTREACTION_PHOTOSYSTEM_II_LHC-II 14 0.0e+00 0.1861126 0.3946925 -0.5209661 0.2279330 0.0105548 0.0007371 0.6795776 0.4799126 1.00e-07
185 RNA_REGULATION_OF_TRANSCRIPTION_MYB_DOMAIN_TRANSCRIPTION_FACTOR_FAMILY 121 0.0e+00 0.0812981 -0.2043840 -0.2470646 0.1224722 0.0001030 0.0000027 0.3307916 0.1785394 3.00e-07
108 PROTEIN_DEGRADATION_UBIQUITIN_E2 38 0.0e+00 0.1794552 -0.2956634 0.1368396 0.0555798 0.0016095 0.1443744 0.3719490 0.2628728 4.00e-07
88 NOT_ASSIGNED_NO_ONTOLOGY_GLYCINE_RICH_PROTEINS 55 0.0e+00 -0.2697929 -0.3713466 -0.4505542 0.0005375 0.0000019 0.0000000 0.6431839 0.0906105 4.00e-07
165 RNA_REGULATION_OF_TRANSCRIPTION_BHLH,BASIC_HELIX-LOOP-HELIX_FAMILY 120 0.0e+00 0.1152845 -0.1998834 -0.1914765 0.0291695 0.0001554 0.0002911 0.2998452 0.1795846 6.00e-07
52 MICRO_RNA,_NATURAL_ANTISENSE_ETC 201 2.0e-07 0.0749339 -0.1662268 -0.0648100 0.0670504 0.0000484 0.1132170 0.1935117 0.1210869 3.30e-06
104 PROTEIN_DEGRADATION_METALLOPROTEASE 37 3.0e-07 -0.3346538 0.2183788 -0.0101463 0.0004266 0.0215180 0.9149490 0.3997316 0.2779010 3.60e-06
10 CELL_WALL_CELLULOSE_SYNTHESIS 14 5.0e-07 0.5263999 0.5927406 -0.1163061 0.0006482 0.0001227 0.4511721 0.8012274 0.3916247 6.40e-06
78 MITOCHONDRIAL_ELECTRON_TRANSPORT_/_ATP_SYNTHESIS_F1-ATPASE 20 6.0e-07 0.6517189 0.1681371 0.3626064 0.0000004 0.1930114 0.0049926 0.7645202 0.2433296 7.60e-06
168 RNA_REGULATION_OF_TRANSCRIPTION_C2C2(ZN)_DOF_ZINC_FINGER_FAMILY 31 9.0e-07 0.2344037 -0.2650134 -0.4000291 0.0238919 0.0106530 0.0001156 0.5340417 0.3342035 1.10e-05
186 RNA_REGULATION_OF_TRANSCRIPTION_MYB-RELATED_TRANSCRIPTION_FACTOR_FAMILY 42 1.0e-06 0.0546069 -0.1793903 -0.4518395 0.5403245 0.0442579 0.0000004 0.4892052 0.2534664 1.11e-05
173 RNA_REGULATION_OF_TRANSCRIPTION_CHROMATIN_REMODELING_FACTORS 35 1.0e-06 -0.3666201 -0.0808627 -0.3783294 0.0001740 0.4077425 0.0001071 0.5329936 0.1684640 1.11e-05
118 PROTEIN_POSTRANSLATIONAL_MODIFICATION 505 1.1e-06 -0.0773232 0.0529974 -0.0324519 0.0029407 0.0415321 0.2120316 0.0992005 0.0662048 1.13e-05
11 CELL_WALL_CELLULOSE_SYNTHESIS_CELLULOSE_SYNTHASE 24 1.3e-06 -0.0120851 0.4942652 0.0070088 0.9183745 0.0000276 0.9526060 0.4944626 0.2869884 1.32e-05
79 MITOCHONDRIAL_ELECTRON_TRANSPORT_/_ATP_SYNTHESIS_NADH-DH_LOCALISATION_NOT_CLEAR 33 1.3e-06 0.4985982 0.2220720 0.3395594 0.0000007 0.0272532 0.0007349 0.6428194 0.1387824 1.32e-05
134 PROTEIN_TARGETING_CHLOROPLAST 36 1.7e-06 -0.4328925 0.0106237 0.1362885 0.0000069 0.9121674 0.1570381 0.4539640 0.2990165 1.57e-05
110 PROTEIN_DEGRADATION_UBIQUITIN_E3_RING 374 1.8e-06 -0.1628086 -0.0663701 -0.0559343 0.0000001 0.0275576 0.0633083 0.1845001 0.0589229 1.60e-05
13 CELL_WALL_CELL_WALL_PROTEINS_HRGP 14 2.3e-06 -0.1955227 -0.6041254 -0.8034698 0.2052754 0.0000906 0.0000002 1.0240900 0.3099178 2.03e-05 
rownames(top) <- top[,1]
top <- top[,4:6]
colfunc <- colorRampPalette(c("blue", "white", "red"))
colnames(top) <- gsub("s\\.","",colnames(top))

heatmap.2(  as.matrix(top), col=colfunc(25), 
    scale="none",Colv=FALSE,trace="none", dendrogram="row",
    margins = c(10,30), cexCol=0.5 , cexRow=0.5, main="Top genes sets by p-value")

# effect size
res_eff <- mitch_calc(x=m,genesets=genesets,priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
top <- head(subset(res_eff$enrichment_result,p.adjustMANOVA<0.05),50)
kbl(head(top,30)) %>% kable_styling()
set setSize pMANOVA s.S93 s.S94 s.S80 p.S93 p.S94 p.S80 s.dist SD p.adjustMANOVA
172 RNA_REGULATION_OF_TRANSCRIPTION_CCAAT_BOX_BINDING_FACTOR_FAMILY,_HAP2 10 0.0000738 -0.2732572 -0.6907347 -0.8002632 0.1345778 0.0001551 0.0000117 1.0918815 0.2780945 0.0003579
153 REDOX_THIOREDOXIN_PDIL 13 0.0000557 0.4920328 0.6378648 0.6511085 0.0021265 0.0000681 0.0000479 1.0358137 0.0882680 0.0003013
13 CELL_WALL_CELL_WALL_PROTEINS_HRGP 14 0.0000023 -0.1955227 -0.6041254 -0.8034698 0.2052754 0.0000906 0.0000002 1.0240900 0.3099178 0.0000203
251 TRANSPORT_MAJOR_INTRINSIC_PROTEINS_PIP 13 0.0002546 0.4479549 0.6066761 0.5947488 0.0051625 0.0001519 0.0002044 0.9604403 0.0883960 0.0010882
20 DEVELOPMENT_LATE_EMBRYOGENESIS_ABUNDANT 18 0.0000743 0.2187946 0.5334065 0.5854848 0.1080369 0.0000891 0.0000170 0.8216970 0.1983912 0.0003579
10 CELL_WALL_CELLULOSE_SYNTHESIS 14 0.0000005 0.5263999 0.5927406 -0.1163061 0.0006482 0.0001227 0.4511721 0.8012274 0.3916247 0.0000064
128 PROTEIN_SYNTHESIS_RIBOSOMAL_RNA 15 0.0000613 0.6871544 0.2616325 0.2850793 0.0000040 0.0793575 0.0559191 0.7886082 0.2391941 0.0003183
188 RNA_REGULATION_OF_TRANSCRIPTION_NUCLEOSOME/CHROMATIN_ASSEMBLY_FACTOR_GROUP 13 0.0008450 -0.0674912 -0.5582221 -0.5318216 0.6735141 0.0004915 0.0008987 0.7739517 0.2760183 0.0029083
78 MITOCHONDRIAL_ELECTRON_TRANSPORT_/_ATP_SYNTHESIS_F1-ATPASE 20 0.0000006 0.6517189 0.1681371 0.3626064 0.0000004 0.1930114 0.0049926 0.7645202 0.2433296 0.0000076
91 NOT_ASSIGNED_NO_ONTOLOGY_PROLINE_RICH_FAMILY 38 0.0000000 -0.0392835 -0.5140553 -0.5498953 0.6751904 0.0000000 0.0000000 0.7537778 0.2850196 0.0000000
26 DNA_SYNTHESIS/CHROMATIN_STRUCTURE_HISTONE_CORE_H2A 14 0.0006407 0.5562959 0.2184265 0.4055399 0.0003129 0.1570537 0.0086042 0.7222450 0.1692604 0.0023581
174 RNA_REGULATION_OF_TRANSCRIPTION_G2-LIKE_TRANSCRIPTION_FACTOR_FAMILY,_GARP 37 0.0000000 -0.0823863 -0.5920351 -0.4031105 0.3858352 0.0000000 0.0000220 0.7209655 0.2576492 0.0000000
39 HORMONE_METABOLISM_JASMONATE_INDUCED-REGULATED-RESPONSIVE-ACTIVATED 13 0.0013418 0.3910455 0.5808904 0.1711186 0.0146321 0.0002868 0.2854031 0.7208549 0.2050698 0.0043900
138 PROTEIN_TARGETING_SECRETORY_PATHWAY_ER 14 0.0006586 0.4545650 0.2192997 0.5105956 0.0032288 0.1554019 0.0009390 0.7179342 0.1545652 0.0023584
75 MITOCHONDRIAL_ELECTRON_TRANSPORT_/_ATP_SYNTHESIS_CYTOCHROME_C 10 0.0180976 0.4681239 0.2835895 0.4501763 0.0103633 0.1204538 0.0136949 0.7086760 0.1017564 0.0361612
87 NOT_ASSIGNED_NO_ONTOLOGY_FORMIN_HOMOLOGY_2_DOMAIN-CONTAINING_PROTEIN 16 0.0007986 -0.1948030 -0.3068310 -0.5883881 0.1773123 0.0335886 0.0000459 0.6915881 0.2027864 0.0028217
209 SIGNALLING_14-3-3_PROTEINS 11 0.0073188 0.4501350 0.2105200 0.4671784 0.0097319 0.2266695 0.0072941 0.6820527 0.1435150 0.0179582
14 CELL_WALL_CELL_WALL_PROTEINS_LRR 17 0.0010390 -0.1931445 -0.3283396 -0.5628537 0.1679790 0.0190838 0.0000586 0.6796440 0.1870648 0.0035300
144 PS_LIGHTREACTION_PHOTOSYSTEM_II_LHC-II 14 0.0000000 0.1861126 0.3946925 -0.5209661 0.2279330 0.0105548 0.0007371 0.6795776 0.4799126 0.0000001
160 RNA_REGULATION_OF_TRANSCRIPTION_ARF,_AUXIN_RESPONSE_FACTOR_FAMILY 16 0.0026024 -0.2512167 -0.3240995 -0.5359595 0.0818967 0.0247948 0.0002054 0.6748354 0.1479161 0.0078369
242 STRESS_BIOTIC_PR-PROTEINS_PLANT_DEFENSINS 13 0.0001118 0.4622835 0.0020397 0.4478247 0.0038994 0.9898405 0.0051755 0.6436281 0.2616479 0.0005198
88 NOT_ASSIGNED_NO_ONTOLOGY_GLYCINE_RICH_PROTEINS 55 0.0000000 -0.2697929 -0.3713466 -0.4505542 0.0005375 0.0000019 0.0000000 0.6431839 0.0906105 0.0000004
79 MITOCHONDRIAL_ELECTRON_TRANSPORT_/_ATP_SYNTHESIS_NADH-DH_LOCALISATION_NOT_CLEAR 33 0.0000013 0.4985982 0.2220720 0.3395594 0.0000007 0.0272532 0.0007349 0.6428194 0.1387824 0.0000132
163 RNA_REGULATION_OF_TRANSCRIPTION_AUX/IAA_FAMILY 22 0.0000493 0.5582094 0.3107696 0.0665666 0.0000058 0.0116213 0.5888659 0.6423446 0.2458232 0.0002722
227 SIGNALLING_RECEPTOR_KINASES_LEUCINE_RICH_REPEAT_X 12 0.0002626 -0.2642095 0.4842900 0.3117233 0.1130200 0.0036721 0.0615144 0.6336520 0.3919459 0.0010896
48 LIPID_METABOLISM_PHOSPHOLIPID_SYNTHESIS 16 0.0132641 0.2237622 0.4541778 0.3682525 0.1212280 0.0016573 0.0107595 0.6260646 0.1164417 0.0276771
124 PROTEIN_SYNTHESIS_ELONGATION 29 0.0000070 0.0022584 0.4575626 0.4268957 0.9832063 0.0000199 0.0000690 0.6257863 0.2544796 0.0000514
183 RNA_REGULATION_OF_TRANSCRIPTION_MADS_BOX_TRANSCRIPTION_FACTOR_FAMILY 28 0.0001533 -0.1825688 -0.4780858 -0.3586854 0.0944994 0.0000119 0.0010181 0.6249421 0.1486628 0.0006885
82 NOT_ASSIGNED_NO_ONTOLOGY_AGENET_DOMAIN-CONTAINING_PROTEIN 15 0.0088425 -0.4661024 -0.2449752 -0.3204251 0.0017733 0.1004408 0.0316574 0.6163899 0.1124069 0.0202004
196 RNA_REGULATION_OF_TRANSCRIPTION_TRIHELIX,_TRIPLE-HELIX_TRANSCRIPTION_FACTOR_FAMILY 29 0.0000272 -0.1438464 -0.2879523 -0.5232703 0.1800034 0.0072730 0.0000011 0.6143452 0.1915305 0.0001676 
rownames(top) <- top[,1]
top <- top[,4:6]
colfunc <- colorRampPalette(c("blue", "white", "red"))
colnames(top) <- gsub("s\\.","",colnames(top))

heatmap.2(  as.matrix(top), col=colfunc(25),
    scale="none",Colv=FALSE,trace="none", dendrogram="row",
    margins = c(10,30), cexCol=0.5 , cexRow=0.5, main="Top genes sets by effect size (FDR<0.05)")

Mitch reports

unlink("mapman_report_sig.html")
capture.output(
    mitch_report(res_sig,outfile=paste("mapman_report_sig.html"))
    , file = "/dev/null", append = FALSE,
    type = c("output", "message"), split = FALSE)
## Dataset saved as " /tmp/RtmpE6FsNj/mapman_report_sig.rds ".
## 
## 
## processing file: mitch.Rmd
## output file: /mnt/bfx4/tohidul/baseline3/mitch.knit.md
## 
## Output created: /tmp/RtmpE6FsNj/mitch_report.html
unlink("mapman_report_eff.html")
capture.output(
    mitch_report(res_eff,outfile=paste("mapman_report_eff.html"))
    , file = "/dev/null", append = FALSE,
    type = c("output", "message"), split = FALSE)
## Dataset saved as " /tmp/RtmpE6FsNj/mapman_report_eff.rds ".
## 
## 
## processing file: mitch.Rmd
## output file: /mnt/bfx4/tohidul/baseline3/mitch.knit.md
## 
## Output created: /tmp/RtmpE6FsNj/mitch_report.html

Session information

So you know what version of R and packages was used.

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.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       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pkgload_1.1.0               GGally_2.0.0               
##  [3] ggplot2_3.3.2               beeswarm_0.2.3             
##  [5] gtools_3.8.2                tibble_3.0.4               
##  [7] dplyr_1.0.2                 echarts4r_0.3.3            
##  [9] UpSetR_1.4.0                eulerr_6.1.0               
## [11] kableExtra_1.3.4            mitch_1.2.2                
## [13] DESeq2_1.30.0               SummarizedExperiment_1.20.0
## [15] Biobase_2.50.0              MatrixGenerics_1.2.0       
## [17] matrixStats_0.57.0          GenomicRanges_1.42.0       
## [19] GenomeInfoDb_1.26.0         IRanges_2.24.0             
## [21] S4Vectors_0.28.0            BiocGenerics_0.36.0        
## [23] gplots_3.1.0                reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6           bit64_4.0.5            webshot_0.5.2         
##  [4] RColorBrewer_1.1-2     httr_1.4.2             rprojroot_1.3-2       
##  [7] backports_1.2.0        tools_4.0.3            R6_2.5.0              
## [10] KernSmooth_2.23-18     DBI_1.1.0              colorspace_2.0-0      
## [13] withr_2.3.0            tidyselect_1.1.0       gridExtra_2.3         
## [16] bit_4.0.4              compiler_4.0.3         rvest_0.3.6           
## [19] xml2_1.3.2             desc_1.2.0             DelayedArray_0.16.0   
## [22] labeling_0.4.2         caTools_1.18.0         scales_1.1.1          
## [25] genefilter_1.72.0      systemfonts_1.0.1      stringr_1.4.0         
## [28] digest_0.6.27          rmarkdown_2.5          svglite_2.0.0         
## [31] XVector_0.30.0         pkgconfig_2.0.3        htmltools_0.5.0       
## [34] highr_0.8              fastmap_1.0.1          htmlwidgets_1.5.2     
## [37] rlang_0.4.8            rstudioapi_0.12        RSQLite_2.2.1         
## [40] shiny_1.5.0            farver_2.0.3           generics_0.1.0        
## [43] jsonlite_1.7.1         BiocParallel_1.24.1    RCurl_1.98-1.2        
## [46] magrittr_1.5           GenomeInfoDbData_1.2.4 Matrix_1.2-18         
## [49] Rcpp_1.0.5             munsell_0.5.0          lifecycle_0.2.0       
## [52] yaml_2.2.1             stringi_1.5.3          MASS_7.3-53           
## [55] zlibbioc_1.36.0        plyr_1.8.6             grid_4.0.3            
## [58] blob_1.2.1             promises_1.1.1         crayon_1.3.4          
## [61] lattice_0.20-41        splines_4.0.3          annotate_1.68.0       
## [64] polylabelr_0.2.0       locfit_1.5-9.4         knitr_1.30            
## [67] pillar_1.4.6           geneplotter_1.68.0     XML_3.99-0.5          
## [70] glue_1.4.2             evaluate_0.14          vctrs_0.3.4           
## [73] httpuv_1.5.4           testthat_3.0.0         polyclip_1.10-0       
## [76] gtable_0.3.0           purrr_0.3.4            assertthat_0.2.1      
## [79] reshape_0.8.8          xfun_0.19              mime_0.9              
## [82] xtable_1.8-4           later_1.1.0.1          survival_3.2-7        
## [85] viridisLite_0.3.0      AnnotationDbi_1.52.0   memoise_1.1.0         
## [88] ellipsis_0.3.1

References

Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification [published correction appears in Nat Biotechnol. 2016 Aug 9;34(8):888]. Nat Biotechnol. 2016;34(5):525-527. doi:10.1038/nbt.3519

Jiang H, Lei R, Ding SW, Zhu S. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinformatics. 2014;15:182. Published 2014 Jun 12. doi:10.1186/1471-2105-15-182

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. doi:10.1186/s13059-014-0550-8

Kaspi A, Ziemann M. mitch: multi-contrast pathway enrichment for multi-omics and single-cell profiling data. BMC Genomics. 2020;21(1):447. Published 2020 Jun 29. doi:10.1186/s12864-020-06856-9