RELIEF: RNA expression analysis
Mark Ziemann
2020-08-14
Introduction
Packages
suppressPackageStartupMessages({
library("tidyverse")
library("reshape2")
library("DESeq2")
library("gplots")
library("fgsea")
library("MASS")
library("mitch")
})Import read counts
tmp <- read.table("3col.tsv",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)curate the samplesheet
Here we're manipulating the sample sheet. Making sure that the samplesheet matches the names of the data sets and vice versa. Ensuring that each patient has a baseline and postop sample.
samplesheet<-read.csv("samplesheet.tsv",sep="\t",stringsAsFactors=FALSE)
mysamples <- colnames(xx)
baseline <- mysamples[grep("T0",mysamples)]
postop <- mysamples[grep("POD",mysamples)]
baseline <- sapply(strsplit(baseline,"-"),"[[",1)
postop <- sapply(strsplit(postop,"-"),"[[",1)
mysamples <- intersect(baseline,postop)
samplesheet2 <- samplesheet
rownames(samplesheet2) <- paste(samplesheet2$RG_number,"T0",sep="-")
samplesheet3 <- samplesheet2
rownames(samplesheet2) <- paste(samplesheet2$RG_number,"T0",sep="-")
rownames(samplesheet3) <- paste(samplesheet2$RG_number,"POD",sep="-")
samplesheet2$timepoint <- 0
samplesheet3$timepoint <- 1
samplesheet <- rbind(samplesheet2,samplesheet3)
# filter the data matrix
xx <- xx[,which(colnames(xx) %in% rownames(samplesheet))]
samplesheet <- samplesheet[which(rownames(samplesheet) %in% colnames(xx)),]
samplesheet <- samplesheet[order(rownames(samplesheet)),]Blood cell composition correction
Here the idea is to use the concept of Monaco et al, 2019 to deconvolute blood RNA-seq signature.
DOI:https://doi.org/10.1016/j.celrep.2019.01.041
xn <- xx
gt <- as.data.frame(sapply(strsplit(rownames(xn)," "),"[[",2) )
rownames(gt) <- rownames(xx)
colnames(gt) = "genesymbol"
gt$geneID <- rownames(xx)
blood <- read.table("https://raw.githubusercontent.com/giannimonaco/ABIS/master/data/sigmatrixRNAseq.txt")
blood2 <- merge(gt,blood,by.x="genesymbol",by.y=0)
blood2 <- blood2[which(!duplicated(blood2$genesymbol)),]
rownames(blood2) <- blood2$geneID
blood2 <- blood2[,c(3:ncol(blood2))]
genes <- intersect(rownames(xx), rownames(blood2))
dec <- apply(xx[genes, , drop=F], 2, function(x) coef(rlm( as.matrix(blood2[genes,]), x, maxit =100 ))) *100## Warning in rlm.default(as.matrix(blood2[genes, ]), x, maxit = 100): 'rlm' failed
## to converge in 100 steps
## Warning in rlm.default(as.matrix(blood2[genes, ]), x, maxit = 100): 'rlm' failed
## to converge in 100 stepsdec <- t(dec/colSums(dec)*100)
dec <- signif(dec, 3)
cs <- colSums(dec)
order(-cs)## [1] 1 4 8 13 2 6 3 14 10 12 5 16 17 9 15 11 7cv <- apply(dec,2,function(x) {sd(x)/mean(x)} )
cv[order(-cv)]## T.CD8.Naive T.CD4.Naive pDCs Monocytes.C B.Naive
## 12.309866 6.935020 6.768844 5.828486 5.773624
## MAIT Monocytes.NC.I NK Basophils.LD B.Memory
## 5.611708 5.598135 5.555667 5.480702 5.163804
## T.CD8.Memory Neutrophils.LD T.CD4.Memory T.gd.non.Vd2 mDCs
## 3.524342 2.972512 -5.555388 -6.493188 -7.442438
## T.gd.Vd2 Plasmablasts
## -7.869309 -18.095614plot(cs,cv,main="cell abundance versus SD")
text(cs,cv+0.5,label=names(cs))
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( dec, col=colfunc(25),scale="row",
trace="none",margins = c(10,15), cexRow=.7, cexCol=.8, main="cell type abundances")
samplesheet <- cbind(samplesheet,dec)
dec_t0_lo <- as.data.frame(dec[which(rownames(dec) %in% rownames(samplesheet)[which(samplesheet$CrpGroup==0 & samplesheet$timepoint==0 )] ),])
dec_po_lo <- as.data.frame(dec[which(rownames(dec) %in% rownames(samplesheet)[which(samplesheet$CrpGroup==0 & samplesheet$timepoint==1 )] ),])
dec_t0_hi <- as.data.frame(dec[which(rownames(dec) %in% rownames(samplesheet)[which(samplesheet$CrpGroup==1 & samplesheet$timepoint==0 )] ),])
dec_po_hi <- as.data.frame(dec[which(rownames(dec) %in% rownames(samplesheet)[which(samplesheet$CrpGroup==1 & samplesheet$timepoint==1 )] ),])
boxplot(dec_t0_lo$Monocytes.C , dec_t0_hi$Monocytes.C, dec_po_lo$Monocytes.C, dec_po_hi$Monocytes.C ,
ylim=c(0,100),names=c("t0 lo","t0 hi","po lo","po hi"),main="Monocytes.C")
t.test(dec_t0_hi$Monocytes.C , dec_po_hi$Monocytes.C)##
## Welch Two Sample t-test
##
## data: dec_t0_hi$Monocytes.C and dec_po_hi$Monocytes.C
## t = -2.7175, df = 25.773, p-value = 0.0116
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -55.076844 -7.627109
## sample estimates:
## mean of x mean of y
## 44.13850 75.49048boxplot(dec_t0_lo$NK , dec_t0_hi$NK, dec_po_lo$NK, dec_po_hi$NK ,
ylim=c(0,15),names=c("t0 lo","t0 hi","po lo","po hi"),main="NK")
t.test(dec_t0_lo$NK , dec_po_lo$NK)##
## Welch Two Sample t-test
##
## data: dec_t0_lo$NK and dec_po_lo$NK
## t = 2.4688, df = 47.994, p-value = 0.01717
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.4820947 4.7136093
## sample estimates:
## mean of x mean of y
## 6.707104 4.109252boxplot(dec_t0_lo$T.CD8.Memory , dec_t0_hi$T.CD8.Memory, dec_po_lo$T.CD8.Memory, dec_po_hi$T.CD8.Memory ,
ylim=c(0,40),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.CD8.Memory")
boxplot(dec_t0_lo$T.CD8.Memory , dec_t0_hi$T.CD8.Memory, dec_po_lo$T.CD8.Memory, dec_po_hi$T.CD8.Memory ,
ylim=c(0,40),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.CD8.Memory")
boxplot(dec_t0_lo$T.CD4.Naive , dec_t0_hi$T.CD4.Naive, dec_po_lo$T.CD4.Naive, dec_po_hi$T.CD4.Naive ,
ylim=c(0,25),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.CD4.Naive")
boxplot(dec_t0_lo$T.CD8.Naive , dec_t0_hi$T.CD8.Naive, dec_po_lo$T.CD8.Naive, dec_po_hi$T.CD8.Naive ,
ylim=c(0,15),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.CD8.Naive")
boxplot(dec_t0_lo$B.Naive , dec_t0_hi$B.Naive, dec_po_lo$B.Naive, dec_po_hi$B.Naive ,
ylim=c(0,10),names=c("t0 lo","t0 hi","po lo","po hi"),main="B.Naive")
boxplot(dec_t0_lo$T.CD4.Memory , dec_t0_hi$T.CD4.Memory, dec_po_lo$T.CD4.Memory, dec_po_hi$T.CD4.Memory ,
ylim=c(0,10),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.CD4.Memory")
boxplot(dec_t0_lo$MAIT , dec_t0_hi$MAIT, dec_po_lo$MAIT, dec_po_hi$MAIT ,
ylim=c(0,15),names=c("t0 lo","t0 hi","po lo","po hi"),main="MAIT")
boxplot(dec_t0_lo$T.gd.Vd2 , dec_t0_hi$T.gd.Vd2 , dec_po_lo$T.gd.Vd2 , dec_po_hi$T.gd.Vd2 ,
ylim=c(0,5),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.gd.Vd2")
boxplot(dec_t0_lo$Neutrophils.LD , dec_t0_hi$Neutrophils.LD , dec_po_lo$Neutrophils.LD , dec_po_hi$Neutrophils.LD ,
ylim=c(0,50),names=c("t0 lo","t0 hi","po lo","po hi"),main="Neutrophils.LD")
boxplot(dec_t0_lo$T.gd.non.Vd2 , dec_t0_hi$T.gd.non.Vd2 , dec_po_lo$T.gd.non.Vd2 , dec_po_hi$T.gd.non.Vd2 ,
ylim=c(0,1),names=c("t0 lo","t0 hi","po lo","po hi"),main="T.gd.non.Vd2")
boxplot(dec_t0_lo$Basophils.LD , dec_t0_hi$Basophils.LD , dec_po_lo$Basophils.LD , dec_po_hi$Basophils.LD ,
ylim=c(0,4),names=c("t0 lo","t0 hi","po lo","po hi"),main="Basophils.LD")
boxplot(dec_t0_lo$Monocytes.NC.I , dec_t0_hi$Monocytes.NC.I , dec_po_lo$Monocytes.NC.I , dec_po_hi$Monocytes.NC.I ,
ylim=c(0,20),names=c("t0 lo","t0 hi","po lo","po hi"),main="Monocytes.NC.I")
boxplot(dec_t0_lo$B.Memory , dec_t0_hi$B.Memory , dec_po_lo$B.Memory , dec_po_hi$B.Memory ,
ylim=c(0,20),names=c("t0 lo","t0 hi","po lo","po hi"),main="B.Memory")
boxplot(dec_t0_lo$mDCs , dec_t0_hi$mDCs , dec_po_lo$mDCs , dec_po_hi$mDCs ,
ylim=c(0,1),names=c("t0 lo","t0 hi","po lo","po hi"),main="mDCs")
boxplot(dec_t0_lo$pDCs , dec_t0_hi$pDCs , dec_po_lo$pDCs , dec_po_hi$pDCs ,
ylim=c(0,1),names=c("t0 lo","t0 hi","po lo","po hi"),main="pDCs")
boxplot(dec_t0_lo$Plasmablasts , dec_t0_hi$Plasmablasts , dec_po_lo$Plasmablasts , dec_po_hi$Plasmablasts ,
ylim=c(0,1),names=c("t0 lo","t0 hi","po lo","po hi"),main="Plasmablasts")
MDS analysis
plot(cmdscale(dist(t(xx))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p", col="gray", pch=19, cex.axis=1.3,cex.lab=1.3, bty='l')
text(cmdscale(dist(t(xx))), labels=colnames(xx),cex=1.3)
heatmap.2(cor(xx),trace="none",margins = c(10,10), main="pearson correlation heatmap")
heatmap.2(cor(xx,method="s"),trace="none",margins = c(10,10), main="spearman correlation heatmap")
DE analysis baseline versus post-op
# filter genes with fewer than 10 reads per sample
NAME = "t0_v_pod"
# no correction for cell types
y <- xx[which(rowSums(xx)/ncol(xx)>=(10)),]
dds <- DESeqDataSetFromMatrix(countData = y , colData = samplesheet,
design = ~ Sex + timepoint )## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## 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() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 157 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## 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),])
dge[1:10,1:6]## baseMean log2FoldChange lfcSE stat
## ENSG00000087116.16 ADAMTS2 634.35845 4.896390 0.3810110 12.85105
## ENSG00000108950.12 FAM20A 792.00129 2.730657 0.2322472 11.75754
## ENSG00000132170.21 PPARG 70.86749 2.706325 0.2315206 11.68935
## ENSG00000170439.7 METTL7B 126.30501 4.034815 0.3456263 11.67393
## ENSG00000121316.11 PLBD1 10689.06659 1.523023 0.1413316 10.77624
## ENSG00000112290.13 WASF1 113.84190 1.280362 0.1209149 10.58895
## ENSG00000163221.9 S100A12 8507.68558 2.298324 0.2194034 10.47533
## ENSG00000156414.19 TDRD9 595.16883 1.691611 0.1624131 10.41549
## ENSG00000137869.15 CYP19A1 28.87688 4.457505 0.4287471 10.39658
## ENSG00000136160.17 EDNRB 34.38674 2.860093 0.2758407 10.36864
## pvalue padj
## ENSG00000087116.16 ADAMTS2 8.485416e-38 1.861700e-33
## ENSG00000108950.12 FAM20A 6.458640e-32 7.085128e-28
## ENSG00000132170.21 PPARG 1.444849e-31 9.502188e-28
## ENSG00000170439.7 METTL7B 1.732395e-31 9.502188e-28
## ENSG00000121316.11 PLBD1 4.457457e-27 1.955932e-23
## ENSG00000112290.13 WASF1 3.353495e-26 1.226261e-22
## ENSG00000163221.9 S100A12 1.121403e-25 3.514798e-22
## ENSG00000156414.19 TDRD9 2.107196e-25 5.778985e-22
## ENSG00000137869.15 CYP19A1 2.569828e-25 6.264670e-22
## ENSG00000136160.17 EDNRB 3.443944e-25 7.556014e-22sig <- subset(dge,padj<0.05)
SIG = nrow(sig)
DN = nrow(subset(sig,log2FoldChange<0))
UP = nrow(subset(sig,log2FoldChange>0))
HEADER = paste(NAME, SIG , "DGEs,", UP ,"upregulated,", DN, "downregulated")
# smear
plot(log2(dge$baseMean),dge$log2FoldChange,cex=0.6,cex.axis=1.2,cex.lab=1.3,
xlab="log2 base mean",
ylim=c(-3,3),ylab="log2 fold change" ,pch=19,col="#838383")
points(log2(sig$baseMean),sig$log2FoldChange,cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
t0_v_pod <- dge
write.table(dge,file="t0_v_pod_rna.tsv",sep="\t",quote=FALSE)
# correct for cell types
y <- xx[which(rowSums(xx)/ncol(xx)>=(10)),]
dds <- DESeqDataSetFromMatrix(countData = y , colData = samplesheet,
design = ~ Sex + Monocytes.C + NK + T.CD8.Memory + T.CD4.Naive + T.CD8.Naive +
B.Naive + T.CD4.Memory + MAIT + T.gd.Vd2 + Neutrophils.LD + T.gd.non.Vd2 +
Basophils.LD + Monocytes.NC.I + B.Memory + mDCs + pDCs + Plasmablasts + timepoint )## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## 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## the design formula contains one or more numeric variables that have mean or
## standard deviation larger than 5 (an arbitrary threshold to trigger this message).
## it is generally a good idea to center and scale numeric variables in the design
## to improve GLM convergence.res <- 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),])
dge[1:10,1:6]## baseMean log2FoldChange lfcSE stat
## ENSG00000087116.16 ADAMTS2 634.35845 4.7023159 0.4039865 11.639784
## ENSG00000108950.12 FAM20A 792.00129 2.5404730 0.2511944 10.113573
## ENSG00000132170.21 PPARG 70.86749 2.4081016 0.2499049 9.636070
## ENSG00000170439.7 METTL7B 126.30501 3.5296298 0.3775196 9.349528
## ENSG00000121316.11 PLBD1 10689.06659 1.4286025 0.1547672 9.230652
## ENSG00000156414.19 TDRD9 595.16883 1.5089966 0.1660858 9.085646
## ENSG00000163221.9 S100A12 8507.68558 1.8947868 0.2118030 8.945985
## ENSG00000164125.15 GASK1B 3752.55250 0.9897545 0.1112864 8.893762
## ENSG00000136052.9 SLC41A2 23.75792 -0.8607596 0.1004404 -8.569856
## ENSG00000163220.11 S100A9 114258.61851 1.3111839 0.1530976 8.564367
## pvalue padj
## ENSG00000087116.16 ADAMTS2 2.586708e-31 5.675237e-27
## ENSG00000108950.12 FAM20A 4.809732e-24 5.276276e-20
## ENSG00000132170.21 PPARG 5.630207e-22 4.117558e-18
## ENSG00000170439.7 METTL7B 8.804006e-21 4.828997e-17
## ENSG00000121316.11 PLBD1 2.689886e-20 1.180322e-16
## ENSG00000156414.19 TDRD9 1.030851e-19 3.769479e-16
## ENSG00000163221.9 S100A12 3.686457e-19 1.155441e-15
## ENSG00000164125.15 GASK1B 5.907390e-19 1.620102e-15
## ENSG00000136052.9 SLC41A2 1.036147e-17 1.881890e-14
## ENSG00000163220.11 S100A9 1.086708e-17 1.881890e-14sig <- subset(dge,padj<0.05)
SIG = nrow(sig)
DN = nrow(subset(sig,log2FoldChange<0))
UP = nrow(subset(sig,log2FoldChange>0))
HEADER = paste(NAME, SIG , "DGEs,", UP ,"upregulated,", DN, "downregulated")
# smear
plot(log2(dge$baseMean),dge$log2FoldChange,cex=0.6,cex.axis=1.2,cex.lab=1.3,
xlab="log2 base mean",
ylim=c(-3,3),ylab="log2 fold change" ,pch=19,col="#838383")
points(log2(sig$baseMean),sig$log2FoldChange,cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
top <- head(sig,20)
# volcano
plot(dge$log2FoldChange, -log2(dge$pvalue) ,cex=0.6, cex.lab=1.3,cex.axis=1.2,
xlim=c(-3,3),xlab="log2 fold change", ylab="-log2 p-value" ,pch=19,col="#838383")
points(sig$log2FoldChange, -log2(sig$pvalue),cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
# top N gene heatmap
colCols <- as.character(samplesheet$timepoint)
colCols <- gsub("0","gray",colCols)
colCols <- gsub("1","black",colCols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( as.matrix(dge[1:50,c(7:ncol(dge))]), col=colfunc(25),scale="row",
ColSideColors=colCols ,
trace="none",margins = c(6,10), cexRow=.6, cexCol=.5, main="Top 50 genes")
#t0_v_pod <- dgeDE analysis post-op low versus high CRP
# filter genes with fewer than 10 reads per sample
NAME = "pod_crp"
ss <- subset(samplesheet,timepoint==1)
xxx <- xx[,colnames(xx) %in% rownames(ss)]
y <- xxx[which(rowSums(xxx)/ncol(xxx)>=(10)),]
dds <- DESeqDataSetFromMatrix(countData = y , colData = ss, design = ~ Sex + CrpGroup )## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## 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() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 123 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## 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),])
dge[1:10,1:6]## baseMean log2FoldChange lfcSE stat
## ENSG00000135424.17 ITGA7 387.1708 2.4417089 0.2463446 9.911760
## ENSG00000108950.12 FAM20A 1327.9835 2.3039461 0.2634579 8.745025
## ENSG00000110079.18 MS4A4A 673.1434 1.8345385 0.2298361 7.981943
## ENSG00000143365.19 RORC 146.5522 -1.4711386 0.1911887 -7.694696
## ENSG00000123643.13 SLC36A1 1610.6451 1.0619669 0.1405202 7.557396
## ENSG00000121316.11 PLBD1 15404.6149 1.3033344 0.1739599 7.492155
## ENSG00000183578.8 TNFAIP8L3 38.9341 3.0442523 0.4078996 7.463239
## ENSG00000014257.16 ACPP 863.6512 0.8117360 0.1088652 7.456341
## ENSG00000018280.17 SLC11A1 17883.6005 1.1672332 0.1586895 7.355455
## ENSG00000166446.15 CDYL2 620.3575 0.8830566 0.1204562 7.330934
## pvalue padj
## ENSG00000135424.17 ITGA7 3.700687e-23 8.016798e-19
## ENSG00000108950.12 FAM20A 2.229669e-18 2.415066e-14
## ENSG00000110079.18 MS4A4A 1.440474e-15 1.040166e-11
## ENSG00000143365.19 RORC 1.418311e-14 7.681219e-11
## ENSG00000123643.13 SLC36A1 4.112185e-14 1.781645e-10
## ENSG00000121316.11 PLBD1 6.775166e-14 2.408881e-10
## ENSG00000183578.8 TNFAIP8L3 8.442088e-14 2.408881e-10
## ENSG00000014257.16 ACPP 8.895834e-14 2.408881e-10
## ENSG00000018280.17 SLC11A1 1.902775e-13 4.579980e-10
## ENSG00000166446.15 CDYL2 2.285536e-13 4.951156e-10sig <- subset(dge,padj<0.05)
SIG = nrow(sig)
DN = nrow(subset(sig,log2FoldChange<0))
UP = nrow(subset(sig,log2FoldChange>0))
HEADER = paste(NAME, SIG , "DGEs,", UP ,"upregulated,", DN, "downregulated")
# smear
plot(log2(dge$baseMean),dge$log2FoldChange,cex=0.6,cex.axis=1.2,cex.lab=1.3, xlab="log2 base mean",
ylim=c(-3,3),ylab="log2 fold change" ,pch=19,col="#838383")
points(log2(sig$baseMean),sig$log2FoldChange,cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
top <- head(sig,20)
# volcano
plot(dge$log2FoldChange, -log2(dge$pvalue) ,cex=0.6, cex.lab=1.3,cex.axis=1.2,
xlim=c(-3,3),xlab="log2 fold change", ylab="-log2 p-value" ,pch=19,col="#838383")
points(sig$log2FoldChange, -log2(sig$pvalue),cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
# top N gene heatmap
colCols <- as.character(ss$CrpGroup)
colCols <- gsub("1","orange",colCols)
colCols <- gsub("0","yellow",colCols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( as.matrix(dge[1:50,c(7:ncol(dge))]), col=colfunc(25),scale="row",
ColSideColors=colCols ,
trace="none",margins = c(6,10), cexRow=.6, cexCol=.5, main="Top 50 genes")
pod_crp <- dge
write.table(dge,file="pod_crp_rna.tsv",sep="\t",quote=FALSE)DE analysis baseline low versus high CRP
# filter genes with fewer than 10 reads per sample
NAME = "t0_crp"
ss <- subset(samplesheet,timepoint==0)
xxx <- xx[,colnames(xx) %in% rownames(ss)]
y <- xxx[which(rowSums(xxx)/ncol(xxx)>=(10)),]
dds <- DESeqDataSetFromMatrix(countData = y , colData = ss, design = ~ Sex + CrpGroup )## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## 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() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 144 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## 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),])
dge[1:10,1:6]## baseMean log2FoldChange lfcSE stat
## ENSG00000129824.16 RPS4Y1 2044.11884 6.1788455 1.0886520 5.675685
## ENSG00000067048.17 DDX3Y 1432.20310 4.9349386 0.9943400 4.963030
## ENSG00000165949.12 IFI27 67.17018 2.3925962 0.4848216 4.935004
## ENSG00000131002.12 TXLNGY 1100.23784 3.3195287 0.6847198 4.848011
## ENSG00000012817.15 KDM5D 1897.13606 2.5522742 0.5268486 4.844417
## ENSG00000169245.6 CXCL10 83.86172 1.4970790 0.3113083 4.808992
## ENSG00000279358.1 AP002833.3 37.34539 0.6674952 0.1393493 4.790085
## ENSG00000183878.15 UTY 427.24603 3.3066722 0.7065698 4.679895
## ENSG00000225886.3 AL445490.1 40.31291 1.0620089 0.2312029 4.593406
## ENSG00000137959.16 IFI44L 2247.61806 1.9332896 0.4322416 4.472706
## pvalue padj
## ENSG00000129824.16 RPS4Y1 1.381350e-08 0.0003060796
## ENSG00000067048.17 DDX3Y 6.940201e-07 0.0052771052
## ENSG00000165949.12 IFI27 8.014928e-07 0.0052771052
## ENSG00000131002.12 TXLNGY 1.247057e-06 0.0052771052
## ENSG00000012817.15 KDM5D 1.269841e-06 0.0052771052
## ENSG00000169245.6 CXCL10 1.516934e-06 0.0052771052
## ENSG00000279358.1 AP002833.3 1.667106e-06 0.0052771052
## ENSG00000183878.15 UTY 2.870223e-06 0.0079498013
## ENSG00000225886.3 AL445490.1 4.360700e-06 0.0107360444
## ENSG00000137959.16 IFI44L 7.723591e-06 0.0171139329sig <- subset(dge,padj<0.05)
SIG = nrow(sig)
DN = nrow(subset(sig,log2FoldChange<0))
UP = nrow(subset(sig,log2FoldChange>0))
HEADER = paste(NAME, SIG , "DGEs,", UP ,"upregulated,", DN, "downregulated")
# smear
plot(log2(dge$baseMean),dge$log2FoldChange,cex=0.6,cex.axis=1.2,cex.lab=1.3, xlab="log2 base mean",
ylim=c(-3,3),ylab="log2 fold change" ,pch=19,col="#838383")
points(log2(sig$baseMean),sig$log2FoldChange,cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
top <- head(sig,20)
# volcano
plot(dge$log2FoldChange, -log2(dge$pvalue) ,cex=0.6, cex.lab=1.3,cex.axis=1.2,
xlim=c(-3,3),xlab="log2 fold change", ylab="-log2 p-value" ,pch=19,col="#838383")
points(sig$log2FoldChange, -log2(sig$pvalue),cex=0.6,pch=19,col="red")
mtext((HEADER),cex=1.2)
# top N gene heatmap
colCols <- as.character(ss$CrpGroup)
colCols <- gsub("1","orange",colCols)
colCols <- gsub("0","yellow",colCols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( as.matrix(dge[1:50,c(7:ncol(dge))]), col=colfunc(25),scale="row",
ColSideColors=colCols ,
trace="none",margins = c(6,10), cexRow=.6, cexCol=.5, main="Top 50 genes")
t0_crp <- dge
write.table(dge,file="t0_crp_rna.tsv",sep="\t",quote=FALSE)One dimensional enrichment analysis
Not correcting for blood compostion changes.
# reactome
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
unzip("ReactomePathways.gmt.zip")
genesets <- gmt_import("ReactomePathways.gmt")
# set up the gene table
gt <- as.data.frame(sapply(strsplit(rownames(xx)," "),"[[",2))
gt$id <- rownames(xx)
colnames(gt) <- c("gene name","id")
# t0_v_pod
# baseline versus post-op
y <- mitch_import(t0_v_pod, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="significance",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set setSize pANOVA
## 736 Neutrophil degranulation 458 3.423751e-103
## 536 Innate Immune System 968 8.433314e-85
## 519 Immune System 1894 1.875234e-46
## 631 Membrane Trafficking 559 7.393463e-36
## 1312 Vesicle-mediated transport 650 3.619480e-32
## 1069 Signal Transduction 1894 8.203706e-28
## 496 Hemostasis 550 4.637545e-24
## 634 Metabolism 1772 3.215227e-22
## 1131 Signaling by Receptor Tyrosine Kinases 414 1.076170e-20
## 829 Platelet activation, signaling and aggregation 221 1.255287e-20
## 1102 Signaling by Interleukins 386 1.856845e-19
## 649 Metabolism of proteins 1719 6.346521e-19
## 841 Post-translational protein modification 1194 1.259555e-18
## 1010 Response to elevated platelet cytosolic Ca2+ 110 1.589264e-16
## 831 Platelet degranulation 106 7.213819e-16
## 1245 Toll-like Receptor Cascades 143 3.059522e-15
## 86 Asparagine N-linked glycosylation 269 8.718083e-15
## 284 Disease 1303 1.028999e-14
## 643 Metabolism of lipids 625 1.259450e-14
## 1288 Transport of small molecules 561 2.089128e-14
## s.dist p.adjustANOVA
## 736 0.5850203 4.669996e-100
## 536 0.3687993 5.751520e-82
## 519 0.1981788 8.526066e-44
## 631 0.3086979 2.521171e-33
## 1312 0.2710023 9.873940e-30
## 1069 0.1515211 1.864976e-25
## 496 0.2519735 9.036588e-22
## 634 0.1385352 5.481961e-20
## 1131 0.2669787 1.630996e-18
## 829 0.3631448 1.712212e-18
## 1102 0.2672402 2.302487e-17
## 649 0.1287882 7.213879e-17
## 841 0.1512473 1.321564e-16
## 1010 0.4549425 1.548397e-14
## 831 0.4531539 6.559766e-14
## 1245 0.3818579 2.608242e-13
## 86 0.2745672 6.994979e-13
## 284 0.1274988 7.797529e-13
## 643 0.1805352 9.041523e-13
## 1288 0.1886709 1.424785e-12unlink("t0_v_pod_mitch.html")
capture.output(
mitch_report(res, "t0_v_pod_mitch.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# pod_crp
# postop low versus high CRP
y <- mitch_import(pod_crp, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="significance",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 734 Neutrophil degranulation
## 534 Innate Immune System
## 1353 rRNA processing in the nucleus and cytosol
## 1351 rRNA processing
## 624 Major pathway of rRNA processing in the nucleolus and cytosol
## 629 Membrane Trafficking
## 517 Immune System
## 349 Eukaryotic Translation Elongation
## 391 Formation of a pool of free 40S subunits
## 806 Peptide chain elongation
## 1058 Selenocysteine synthesis
## 634 Metabolism of RNA
## 351 Eukaryotic Translation Termination
## 1313 Viral mRNA Translation
## 742 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 587 L13a-mediated translational silencing of Ceruloplasmin expression
## 433 GTP hydrolysis and joining of the 60S ribosomal subunit
## 526 Influenza Viral RNA Transcription and Replication
## 1266 Translation
## 1006 Response of EIF2AK4 (GCN2) to amino acid deficiency
## setSize pANOVA s.dist p.adjustANOVA
## 734 458 2.738906e-76 0.5020879 3.727651e-73
## 534 967 4.229759e-52 0.2877620 2.878351e-49
## 1353 190 1.067181e-43 -0.5819485 4.841444e-41
## 1351 217 6.164864e-42 -0.5333709 1.730183e-39
## 624 180 6.356293e-42 -0.5850279 1.730183e-39
## 629 558 1.154911e-32 0.2942613 2.619723e-30
## 517 1893 3.953957e-32 0.1636513 7.687623e-30
## 349 93 7.196091e-32 -0.7037933 1.224235e-29
## 391 100 5.497514e-31 -0.6688533 8.313463e-29
## 806 88 1.446938e-30 -0.7076882 1.969283e-28
## 1058 92 3.335066e-30 -0.6878524 4.126386e-28
## 634 685 9.581179e-30 -0.2535715 1.086665e-27
## 351 92 6.712858e-29 -0.6720014 7.027846e-27
## 1313 88 1.864969e-28 -0.6814373 1.813017e-26
## 742 94 2.177046e-27 -0.6461972 1.975306e-25
## 587 110 4.862350e-27 -0.5935235 4.136037e-25
## 433 111 5.727887e-27 -0.5900315 4.585679e-25
## 526 135 1.026196e-26 -0.5326421 7.759182e-25
## 1266 295 1.470316e-26 -0.3605457 1.053210e-24
## 1006 100 9.854840e-26 -0.6061679 6.706218e-24unlink("pod_crp_mitch.html")
capture.output(
mitch_report(res, "pod_crp_mitch.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# t0_crp
# baseline low versus high CRP
y <- mitch_import(t0_crp, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="significance",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 352 Eukaryotic Translation Elongation
## 809 Peptide chain elongation
## 394 Formation of a pool of free 40S subunits
## 354 Eukaryotic Translation Termination
## 1316 Viral mRNA Translation
## 1061 Selenocysteine synthesis
## 590 L13a-mediated translational silencing of Ceruloplasmin expression
## 436 GTP hydrolysis and joining of the 60S ribosomal subunit
## 745 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 145 Cap-dependent Translation Initiation
## 353 Eukaryotic Translation Initiation
## 1009 Response of EIF2AK4 (GCN2) to amino acid deficiency
## 1269 Translation
## 744 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)
## 746 Nonsense-Mediated Decay (NMD)
## 1041 SRP-dependent cotranslational protein targeting to membrane
## 529 Influenza Viral RNA Transcription and Replication
## 1355 rRNA processing
## 1060 Selenoamino acid metabolism
## 982 Regulation of expression of SLITs and ROBOs
## setSize pANOVA s.dist p.adjustANOVA
## 352 93 1.219276e-46 -0.8584146 1.664312e-43
## 809 88 1.043234e-44 -0.8632455 5.851941e-42
## 394 100 1.286141e-44 -0.8091627 5.851941e-42
## 354 92 1.490608e-42 -0.8229867 5.086699e-40
## 1316 88 6.813190e-42 -0.8346033 1.860001e-39
## 1061 92 5.234678e-41 -0.8073168 1.190889e-38
## 590 110 7.750141e-41 -0.7370186 1.511277e-38
## 436 111 1.284443e-39 -0.7222223 2.191581e-37
## 745 94 1.130965e-38 -0.7747095 1.438121e-36
## 145 118 1.158925e-38 -0.6917297 1.438121e-36
## 353 118 1.158925e-38 -0.6917297 1.438121e-36
## 1009 100 1.335083e-34 -0.7087695 1.518657e-32
## 1269 295 2.649744e-34 -0.4126274 2.782231e-32
## 744 114 1.002228e-33 -0.6551703 9.120273e-32
## 746 114 1.002228e-33 -0.6551703 9.120273e-32
## 1041 111 1.070142e-33 -0.6636255 9.129653e-32
## 529 135 3.517380e-33 -0.5972154 2.824250e-31
## 1355 219 8.947395e-32 -0.4592240 6.785108e-30
## 1060 114 1.590530e-30 -0.6216895 1.142670e-28
## 982 162 8.656005e-30 -0.5154065 5.907723e-28unlink("t0_crp_mitch.html")
capture.output(
mitch_report(res, "t0_crp_mitch.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
## prioritise by effect size
# t0_v_pod
# baseline versus post-op
y <- mitch_import(t0_v_pod, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="effect",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set setSize
## 513 IRAK4 deficiency (TLR2/4) 10
## 686 MyD88 deficiency (TLR2/4) 10
## 1331 alpha-linolenic (omega3) and linoleic (omega6) acid metabolism 12
## 1332 alpha-linolenic acid (ALA) metabolism 12
## 1278 Translocation of ZAP-70 to Immunological synapse 24
## 968 Regulation of TLR by endogenous ligand 11
## 505 Hyaluronan uptake and degradation 12
## 345 Erythrocytes take up carbon dioxide and release oxygen 11
## 765 O2/CO2 exchange in erythrocytes 11
## 1306 Uptake and function of anthrax toxins 10
## 1324 WNT5A-dependent internalization of FZD2, FZD5 and ROR2 11
## 1343 p130Cas linkage to MAPK signaling for integrins 11
## 61 Advanced glycosylation endproduct receptor signaling 12
## 1304 Unwinding of DNA 12
## 902 RNA Polymerase I Promoter Opening 19
## 889 RHO GTPases Activate WASPs and WAVEs 35
## 925 ROS and RNS production in phagocytes 31
## 736 Neutrophil degranulation 458
## 779 PD-1 signaling 28
## 434 GRB2:SOS provides linkage to MAPK signaling for Integrins 12
## pANOVA s.dist p.adjustANOVA
## 513 8.804924e-07 0.8976572 1.115504e-05
## 686 8.804924e-07 0.8976572 1.115504e-05
## 1331 2.312440e-06 0.7873944 2.483597e-05
## 1332 2.312440e-06 0.7873944 2.483597e-05
## 1278 2.762794e-10 -0.7439336 8.374336e-09
## 968 2.021659e-05 0.7420783 1.477979e-04
## 505 9.152072e-06 0.7394611 7.520137e-05
## 345 3.557086e-05 0.7198077 2.401913e-04
## 765 3.557086e-05 0.7198077 2.401913e-04
## 1306 1.130846e-04 0.7049367 6.337001e-04
## 1324 9.664711e-05 0.6787956 5.585875e-04
## 1343 1.393845e-04 0.6632012 7.498444e-04
## 61 9.451769e-05 0.6508107 5.509493e-04
## 1304 1.039696e-04 -0.6469590 5.983734e-04
## 902 2.758607e-06 0.6210925 2.838780e-05
## 889 5.295008e-10 0.6062964 1.536679e-08
## 925 7.976405e-09 0.5984434 1.754809e-07
## 736 3.423751e-103 0.5850203 4.669996e-100
## 779 1.170708e-07 -0.5782779 1.971414e-06
## 434 5.274371e-04 0.5778488 2.398081e-03unlink("t0_v_pod_mitche.html")
capture.output(
mitch_report(res, "t0_v_pod_mitche.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# pod_crp
# postop low versus high CRP
y <- mitch_import(pod_crp, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="effect",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 511 IRAK4 deficiency (TLR2/4)
## 684 MyD88 deficiency (TLR2/4)
## 1328 alpha-linolenic (omega3) and linoleic (omega6) acid metabolism
## 1329 alpha-linolenic acid (ALA) metabolism
## 1276 Translocation of ZAP-70 to Immunological synapse
## 806 Peptide chain elongation
## 349 Eukaryotic Translation Elongation
## 1058 Selenocysteine synthesis
## 1313 Viral mRNA Translation
## 966 Regulation of TLR by endogenous ligand
## 351 Eukaryotic Translation Termination
## 391 Formation of a pool of free 40S subunits
## 503 Hyaluronan uptake and degradation
## 1321 WNT5A-dependent internalization of FZD2, FZD5 and ROR2
## 742 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 886 RHO GTPases Activate ROCKs
## 1304 Uptake and function of anthrax toxins
## 618 MET activates RAP1 and RAC1
## 1308 VLDLR internalisation and degradation
## 705 NOTCH4 Activation and Transmission of Signal to the Nucleus
## setSize pANOVA s.dist p.adjustANOVA
## 511 10 1.738461e-06 0.8730308 2.190782e-05
## 684 10 1.738461e-06 0.8730308 2.190782e-05
## 1328 12 4.980696e-06 0.7610050 5.422981e-05
## 1329 12 4.980696e-06 0.7610050 5.422981e-05
## 1276 24 3.577379e-10 -0.7392126 9.186438e-09
## 806 88 1.446938e-30 -0.7076882 1.969283e-28
## 349 93 7.196091e-32 -0.7037933 1.224235e-29
## 1058 92 3.335066e-30 -0.6878524 4.126386e-28
## 1313 88 1.864969e-28 -0.6814373 1.813017e-26
## 966 11 1.083937e-04 0.6739467 7.764415e-04
## 351 92 6.712858e-29 -0.6720014 7.027846e-27
## 391 100 5.497514e-31 -0.6688533 8.313463e-29
## 503 12 8.491297e-05 0.6551182 6.456232e-04
## 1321 11 1.720915e-04 0.6540745 1.131481e-03
## 742 94 2.177046e-27 -0.6461972 1.975306e-25
## 886 18 2.749013e-06 0.6381961 3.225350e-05
## 1304 10 5.013788e-04 0.6354563 2.740468e-03
## 618 10 5.664645e-04 0.6294621 3.071547e-03
## 1308 12 2.062624e-04 0.6186456 1.305689e-03
## 705 10 7.044365e-04 0.6186300 3.659306e-03unlink("pod_crp_mitche.html")
capture.output(
mitch_report(res, "pod_crp_mitche.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# t0_crp
# baseline low versus high CRP
y <- mitch_import(t0_crp, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(y, genesets, priority="effect",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 809 Peptide chain elongation
## 352 Eukaryotic Translation Elongation
## 1316 Viral mRNA Translation
## 354 Eukaryotic Translation Termination
## 394 Formation of a pool of free 40S subunits
## 1061 Selenocysteine synthesis
## 745 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 590 L13a-mediated translational silencing of Ceruloplasmin expression
## 436 GTP hydrolysis and joining of the 60S ribosomal subunit
## 1009 Response of EIF2AK4 (GCN2) to amino acid deficiency
## 145 Cap-dependent Translation Initiation
## 353 Eukaryotic Translation Initiation
## 1041 SRP-dependent cotranslational protein targeting to membrane
## 744 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)
## 746 Nonsense-Mediated Decay (NMD)
## 400 Formation of the ternary complex, and subsequently, the 43S complex
## 1060 Selenoamino acid metabolism
## 529 Influenza Viral RNA Transcription and Replication
## 1270 Translation initiation complex formation
## 50 Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S
## setSize pANOVA s.dist p.adjustANOVA
## 809 88 1.043234e-44 -0.8632455 5.851941e-42
## 352 93 1.219276e-46 -0.8584146 1.664312e-43
## 1316 88 6.813190e-42 -0.8346033 1.860001e-39
## 354 92 1.490608e-42 -0.8229867 5.086699e-40
## 394 100 1.286141e-44 -0.8091627 5.851941e-42
## 1061 92 5.234678e-41 -0.8073168 1.190889e-38
## 745 94 1.130965e-38 -0.7747095 1.438121e-36
## 590 110 7.750141e-41 -0.7370186 1.511277e-38
## 436 111 1.284443e-39 -0.7222223 2.191581e-37
## 1009 100 1.335083e-34 -0.7087695 1.518657e-32
## 145 118 1.158925e-38 -0.6917297 1.438121e-36
## 353 118 1.158925e-38 -0.6917297 1.438121e-36
## 1041 111 1.070142e-33 -0.6636255 9.129653e-32
## 744 114 1.002228e-33 -0.6551703 9.120273e-32
## 746 114 1.002228e-33 -0.6551703 9.120273e-32
## 400 51 1.280989e-15 -0.6467870 4.995858e-14
## 1060 114 1.590530e-30 -0.6216895 1.142670e-28
## 529 135 3.517380e-33 -0.5972154 2.824250e-31
## 1270 58 4.367546e-15 -0.5950456 1.611270e-13
## 50 59 2.907221e-15 -0.5938195 1.102321e-13unlink("t0_crp_mitche.html")
capture.output(
mitch_report(res, "t0_crp_mitche.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# 3D analysys significance
z <- list("t0_v_pod"=t0_v_pod, "pod_crp"=pod_crp, "t0_crp"=t0_crp)
zz <- mitch_import(z, DEtype="DESeq2", geneTable=gt)
capture.output(
res <- mitch_calc(zz, genesets, priority="significance",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 734 Neutrophil degranulation
## 534 Innate Immune System
## 391 Formation of a pool of free 40S subunits
## 349 Eukaryotic Translation Elongation
## 587 L13a-mediated translational silencing of Ceruloplasmin expression
## 1265 Translation
## 433 GTP hydrolysis and joining of the 60S ribosomal subunit
## 806 Peptide chain elongation
## 142 Cap-dependent Translation Initiation
## 350 Eukaryotic Translation Initiation
## 351 Eukaryotic Translation Termination
## 1352 rRNA processing in the nucleus and cytosol
## 1312 Viral mRNA Translation
## 1058 Selenocysteine synthesis
## 1350 rRNA processing
## 624 Major pathway of rRNA processing in the nucleolus and cytosol
## 742 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 1006 Response of EIF2AK4 (GCN2) to amino acid deficiency
## 517 Immune System
## 634 Metabolism of RNA
## setSize pMANOVA s.t0_v_pod s.pod_crp s.t0_crp p.t0_v_pod
## 734 458 2.958869e-111 0.58479156 0.5071212 -0.12455422 4.384529e-103
## 534 966 5.488806e-90 0.36891702 0.2930445 -0.09023284 1.302536e-84
## 391 100 6.506201e-71 -0.21172313 -0.6656200 -0.80776843 2.532751e-04
## 349 93 6.848264e-69 -0.27690192 -0.7009300 -0.85714704 3.912367e-06
## 587 110 4.551002e-67 -0.14656489 -0.5899641 -0.73549562 7.914645e-03
## 1265 295 5.628229e-67 -0.05871998 -0.3551876 -0.41008815 8.278961e-02
## 433 111 1.689386e-65 -0.15065842 -0.5864369 -0.72072157 6.103043e-03
## 806 88 2.051258e-65 -0.28523245 -0.7048806 -0.86200437 3.731219e-06
## 142 118 2.139637e-63 -0.13507988 -0.5525330 -0.69008790 1.125976e-02
## 350 118 2.139637e-63 -0.13507988 -0.5525330 -0.69008790 1.125976e-02
## 351 92 4.919890e-62 -0.26775922 -0.6690044 -0.82164804 9.015423e-06
## 1352 190 1.367402e-61 -0.30939749 -0.5783771 -0.47102270 1.865520e-13
## 1312 88 1.894932e-61 -0.26588337 -0.6785649 -0.83331803 1.618729e-05
## 1058 92 2.275161e-60 -0.30282224 -0.6849465 -0.80593775 5.134641e-07
## 1350 217 2.991768e-60 -0.29097363 -0.5293168 -0.45628521 1.461443e-13
## 624 180 8.401547e-59 -0.30873491 -0.5814270 -0.46842864 8.809051e-13
## 742 94 4.899324e-58 -0.24781882 -0.6431885 -0.77333684 3.281129e-05
## 1006 100 1.766879e-53 -0.22533903 -0.6027217 -0.70729571 9.847802e-05
## 517 1888 3.670481e-52 0.19752101 0.1685980 -0.06640310 6.263035e-46
## 634 685 7.677541e-52 -0.08970780 -0.2478499 -0.23379009 6.294602e-05
## p.pod_crp p.t0_crp s.dist SD p.adjustMANOVA
## 734 8.676194e-78 4.922943e-06 0.7840069 0.38906255 4.024062e-108
## 534 6.748707e-54 2.059618e-06 0.4797049 0.24612895 3.732388e-87
## 391 1.057544e-30 1.810813e-44 1.0678795 0.31131403 2.949478e-68
## 349 1.268958e-31 1.654543e-46 1.1413495 0.30024656 2.328410e-66
## 587 9.819067e-27 1.129424e-40 0.9541974 0.30676313 1.237873e-64
## 1265 8.085094e-26 6.857521e-34 0.5456909 0.18901794 1.275732e-64
## 433 1.167321e-26 1.854683e-39 0.9413001 0.29802329 3.282236e-63
## 806 2.456611e-30 1.386306e-44 1.1494632 0.29817733 3.487138e-63
## 142 3.071357e-25 1.745814e-38 0.8942933 0.28902778 2.909907e-61
## 350 3.071357e-25 1.745814e-38 0.8942933 0.28902778 2.909907e-61
## 351 1.176274e-28 2.026207e-42 1.0928712 0.28609165 6.082773e-60
## 1352 3.515546e-43 3.591287e-29 0.8075328 0.13539924 1.549722e-59
## 1312 3.141898e-28 9.065273e-42 1.1070516 0.29332481 1.982391e-59
## 1058 5.819109e-30 7.141846e-41 1.1001767 0.26260968 2.210157e-58
## 1350 2.535060e-41 4.106094e-31 0.7569928 0.12211265 2.712536e-58
## 624 1.997079e-41 1.915500e-27 0.8079604 0.13701075 7.141315e-57
## 742 3.786007e-27 1.532860e-38 1.0359322 0.27368624 3.919459e-56
## 1006 1.854357e-25 1.836975e-34 0.9561999 0.25352042 1.334975e-51
## 517 6.839889e-34 1.817579e-06 0.2680470 0.14475152 2.627291e-50
## 634 1.767643e-28 1.632814e-25 0.3523278 0.08752742 5.220728e-50unlink("3d_mitch.html")
capture.output(
mitch_report(res, "3d_mitch.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
# 3D analysys effect size
capture.output(
res <- mitch_calc(zz, genesets, priority="effect",resrows=100)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)
head(res$enrichment_result,20)## set
## 511 IRAK4 deficiency (TLR2/4)
## 684 MyD88 deficiency (TLR2/4)
## 806 Peptide chain elongation
## 349 Eukaryotic Translation Elongation
## 1312 Viral mRNA Translation
## 1058 Selenocysteine synthesis
## 1327 alpha-linolenic (omega3) and linoleic (omega6) acid metabolism
## 1328 alpha-linolenic acid (ALA) metabolism
## 351 Eukaryotic Translation Termination
## 1275 Translocation of ZAP-70 to Immunological synapse
## 1320 WNT5A-dependent internalization of FZD2, FZD5 and ROR2
## 391 Formation of a pool of free 40S subunits
## 742 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 966 Regulation of TLR by endogenous ligand
## 1303 Uptake and function of anthrax toxins
## 503 Hyaluronan uptake and degradation
## 1006 Response of EIF2AK4 (GCN2) to amino acid deficiency
## 587 L13a-mediated translational silencing of Ceruloplasmin expression
## 433 GTP hydrolysis and joining of the 60S ribosomal subunit
## 1307 VLDLR internalisation and degradation
## setSize pMANOVA s.t0_v_pod s.pod_crp s.t0_crp p.t0_v_pod
## 511 10 7.837284e-06 0.8977867 0.8758224 0.234972626 8.773124e-07
## 684 10 7.837284e-06 0.8977867 0.8758224 0.234972626 8.773124e-07
## 806 88 2.051258e-65 -0.2852324 -0.7048806 -0.862004372 3.731219e-06
## 349 93 6.848264e-69 -0.2769019 -0.7009300 -0.857147042 3.912367e-06
## 1312 88 1.894932e-61 -0.2658834 -0.6785649 -0.833318035 1.618729e-05
## 1058 92 2.275161e-60 -0.3028222 -0.6849465 -0.805937747 5.134641e-07
## 1327 12 1.992636e-05 0.7872073 0.7652908 -0.005514218 2.325350e-06
## 1328 12 1.992636e-05 0.7872073 0.7652908 -0.005514218 2.325350e-06
## 351 92 4.919890e-62 -0.2677592 -0.6690044 -0.821648044 9.015423e-06
## 1275 24 6.089873e-10 -0.7436984 -0.7368076 -0.285265565 2.799114e-10
## 1320 11 1.132446e-06 0.6777130 0.6590122 -0.528004424 9.916521e-05
## 391 100 6.506201e-71 -0.2117231 -0.6656200 -0.807768432 2.532751e-04
## 742 94 4.899324e-58 -0.2478188 -0.6431885 -0.773336838 3.281129e-05
## 966 11 3.182361e-04 0.7417195 0.6773982 0.061283873 2.040495e-05
## 1303 10 1.726190e-04 0.7039633 0.6406813 -0.300977961 1.155837e-04
## 503 10 4.528317e-04 0.7266108 0.6675776 -0.138187263 6.909383e-05
## 1006 100 1.766879e-53 -0.2253390 -0.6027217 -0.707295710 9.847802e-05
## 587 110 4.551002e-67 -0.1465649 -0.5899641 -0.735495617 7.914645e-03
## 433 111 1.689386e-65 -0.1506584 -0.5864369 -0.720721568 6.103043e-03
## 1307 11 2.053616e-06 0.4667035 0.5913302 -0.558650615 7.353721e-03
## p.pod_crp p.t0_crp s.dist SD p.adjustMANOVA
## 511 1.610878e-06 1.982174e-01 1.2760478 0.3764955 4.368322e-05
## 684 1.610878e-06 1.982174e-01 1.2760478 0.3764955 4.368322e-05
## 806 2.456611e-30 1.386306e-44 1.1494632 0.2981773 3.487138e-63
## 349 1.268958e-31 1.654543e-46 1.1413495 0.3002466 2.328410e-66
## 1312 3.141898e-28 9.065273e-42 1.1070516 0.2933248 1.982391e-59
## 1058 5.819109e-30 7.141846e-41 1.1001767 0.2626097 2.210157e-58
## 1327 4.404524e-06 9.736155e-01 1.0979052 0.4514843 1.003698e-04
## 1328 4.404524e-06 9.736155e-01 1.0979052 0.4514843 1.003698e-04
## 351 1.176274e-28 2.026207e-42 1.0928712 0.2860917 6.082773e-60
## 1275 4.078332e-10 1.555367e-02 1.0850572 0.2627097 6.844816e-09
## 1320 1.536000e-04 2.425186e-03 1.0827653 0.6907861 8.280252e-06
## 391 1.057544e-30 1.810813e-44 1.0678795 0.3113140 2.949478e-68
## 742 3.786007e-27 1.532860e-38 1.0359322 0.2736862 3.919459e-56
## 966 9.990741e-05 7.248845e-01 1.0063657 0.3756609 1.151067e-03
## 1303 4.504122e-04 9.933326e-02 0.9983109 0.5628253 6.824471e-04
## 503 2.562060e-04 4.492496e-01 0.9963528 0.4831524 1.591347e-03
## 1006 1.854357e-25 1.836975e-34 0.9561999 0.2535204 1.334975e-51
## 587 9.819067e-27 1.129424e-40 0.9541974 0.3067631 1.237873e-64
## 433 1.167321e-26 1.854683e-39 0.9413001 0.2980233 3.282236e-63
## 1307 6.829866e-04 1.333938e-03 0.9378561 0.6310492 1.352994e-05unlink("3d_mitche.html")
capture.output(
mitch_report(res, "3d_mitche.html")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)