Ee are interested in comparing DEGs between the following comparison groups:
Latently- vs productively-infected cells (groups 3 vs 4).
Latently-infected vs bystander cells (2 vs 3).
Productively-infected vs mock (4 vs 1)
Mock vs bystander (1 vs 2).
As discussed, I think we will do these comparisons separately for both MDM and AlvMDM samples initially. Then, it would be of interest to compare DEGs for MDM vs AlvMDM for each of the four infection groups.
suppressPackageStartupMessages({
library("kableExtra")
library("ggplot2")
library("plyr")
library("Seurat")
library("hdf5r")
library("SingleCellExperiment")
library("parallel")
library("stringi")
library("beeswarm")
library("muscat")
library("DESeq2")
library("mitch")
library("harmony")
library("celldex")
library("SingleR")
library("limma")
library("gplots")
})
Sample | Patient | Group |
---|---|---|
mock0 | MDMP236 | mock |
mock1 | MDMV236 | mock |
mock2 | MDMO236 | mock |
mock3 | MDMP239 | mock |
mock4 | AlvP239 | mock |
mock5 | AlvM239 | mock |
mock6 | AlvN239 | mock |
active0 | MDMP236 | active |
active1 | MDMV236 | active |
active2 | MDMO236 | active |
active3 | MDMP239 | active |
active4 | AlvP239 | active |
active5 | AlvM239 | active |
active6 | AlvN239 | active |
latent0 | MDMP236 | latent |
latent1 | MDMV236 | latent |
latent2 | MDMO236 | latent |
latent3 | MDMP239 | latent |
latent4 | AlvP239 | latent |
latent5 | AlvM239 | latent |
latent6 | AlvN239 | latent |
bystander0 | MDMP236 | bystander |
bystander1 | MDMV236 | bystander |
bystander2 | MDMO236 | bystander |
bystander3 | MDMP239 | bystander |
bystander4 | AlvP239 | bystander |
bystander5 | AlvM239 | bystander |
bystander6 | AlvN239 | bystander |
react6 | AlvN239 | reactivated |
exclude react6
ss <- read.table("samplesheet.tsv",header=TRUE,row.names=1)
ss %>% kbl(caption="sample sheet") %>% kable_paper("hover", full_width = F)
Patient | Group | |
---|---|---|
mock0 | MDMP236 | mock |
mock1 | MDMV236 | mock |
mock2 | MDMO236 | mock |
mock3 | MDMP239 | mock |
mock4 | AlvP239 | mock |
mock5 | AlvM239 | mock |
mock6 | AlvN239 | mock |
active0 | MDMP236 | active |
active1 | MDMV236 | active |
active2 | MDMO236 | active |
active3 | MDMP239 | active |
active4 | AlvP239 | active |
active5 | AlvM239 | active |
active6 | AlvN239 | active |
latent0 | MDMP236 | latent |
latent1 | MDMV236 | latent |
latent2 | MDMO236 | latent |
latent3 | MDMP239 | latent |
latent4 | AlvP239 | latent |
latent5 | AlvM239 | latent |
latent6 | AlvN239 | latent |
bystander0 | MDMP236 | bystander |
bystander1 | MDMV236 | bystander |
bystander2 | MDMO236 | bystander |
bystander3 | MDMP239 | bystander |
bystander4 | AlvP239 | bystander |
bystander5 | AlvM239 | bystander |
bystander6 | AlvN239 | bystander |
react6 | AlvN239 | reactivated |
mylist <- readRDS("macrophage_counts_nofilter.rds")
comb <- do.call(cbind,mylist)
sce <- SingleCellExperiment(list(counts=comb))
sce
## class: SingleCellExperiment
## dim: 36603 28971
## metadata(0):
## assays(1): counts
## rownames(36603): gene-HIV1-1 gene-HIV1-2 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(28971): mdm_mock1|AAACGAATCACATACG mdm_mock1|AAACGCTCATCAGCGC
## ... react6|TTTGTTGTCTGAACGT react6|TTTGTTGTCTTGGCTC
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
cellmetadata <- data.frame(colnames(comb) ,sapply(strsplit(colnames(comb),"\\|"),"[[",1))
colnames(cellmetadata) <- c("cell","sample")
comb <- CreateSeuratObject(comb, project = "mac", assay = "RNA", meta.data = cellmetadata)
comb <- NormalizeData(comb)
## Normalizing layer: counts
comb <- FindVariableFeatures(comb, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
comb <- ScaleData(comb)
## Centering and scaling data matrix
comb <- RunPCA(comb, features = VariableFeatures(object = comb))
## PC_ 1
## Positive: MT-ATP8, MTRNR2L1, FABP4, AC105402.3, FABP3, gene-HIV1-2, HAMP, C1orf56, PVALB, AL136097.2
## GAPDH, BIRC7, UTS2, S100A9, C1QC, MT1G, MT1H, AC025154.2, CARD16, COX7A1
## BEX3, MT1A, AL450311.1, PPP1R1A, NUPR1, LILRA1, MT2A, CCL23, FXYD2, CCL15
## Negative: ARL15, NEAT1, DOCK3, EXOC4, MALAT1, RASAL2, FTX, LRMDA, DPYD, PLXDC2
## JMJD1C, TMEM117, FHIT, ZEB2, MITF, DENND4C, ATG7, FRMD4B, VPS13B, ZNF438
## MAML2, COP1, TPRG1, FMNL2, ATXN1, JARID2, ASAP1, ELMO1, FNDC3B, TCF12
## PC_ 2
## Positive: CYSTM1, FAH, CD164, FDX1, GDF15, PSAT1, ATP6V1D, SAT1, TXN, BCAT1
## CCPG1, SLAMF9, CSTA, CTSL, HES2, HEBP2, S100A10, NUPR1, PHGDH, STMN1
## TCEA1, GARS, MARCO, SCD, RILPL2, SNHG5, RHOQ, CD48, PLIN2, SNHG32
## Negative: SPRED1, RCBTB2, KCNMA1, GCLC, FGL2, HLA-DRA, MRC1, HLA-DRB1, RTN1, SLCO2B1
## HLA-DPA1, AC020656.1, CD74, HLA-DPB1, PDGFC, LINC02345, DPYD, MERTK, C1QA, CCDC102B
## STAC, ATP8B4, NFATC3, RNF128, VAMP5, MX1, ATG7, TGFBI, XYLT1, CRIP1
## PC_ 3
## Positive: NCAPH, CRABP2, TM4SF19, GAL, DUSP2, RGCC, CHI3L1, CCL22, TMEM114, ACAT2
## RGS20, TRIM54, LINC01010, LINC02244, TCTEX1D2, MGST1, GPC4, NUMB, CCND1, SYNGR1
## CLU, POLE4, SERTAD2, AC092353.2, SLC20A1, IL1RN, ZNF366, AC005280.2, OCSTAMP, GCLC
## Negative: MARCKS, CD14, TGFBI, BTG1, MS4A6A, HLA-DQA1, FPR3, CTSC, CD163, MPEG1
## MEF2C, AIF1, TIMP1, HLA-DPB1, C1QC, HLA-DQB1, SELENOP, RNASE1, NAMPT, HLA-DQA2
## TCF4, EPB41L3, ARL4C, NUPR1, ALDH2, SAMSN1, ZNF331, HIF1A, MS4A7, C1QA
## PC_ 4
## Positive: MT-ATP8, MTRNR2L1, XIST, PDE4D, AL035446.2, AC073359.2, CCDC26, AL365295.1, CLVS1, PRKCE
## EPHB1, RAD51B, C5orf17, LINC01135, PLEKHA5, LINC02320, LY86-AS1, STX18-AS1, FHIT, LINC01473
## DMGDH, FTX, AC079142.1, PKD1L1, SLC22A15, AF117829.1, GUCY2F, SH3RF3, STOX2, MIR646HG
## Negative: HLA-DRB1, CD74, CTSB, HSP90B1, TUBA1B, GAPDH, HLA-DPA1, HLA-DRA, ACTG1, IFI30
## RGCC, AIF1, HSPA5, CYP1B1, HLA-DPB1, C1QA, C15orf48, LYZ, CA2, S100A4
## LGMN, TUBB, FBP1, TUBA1A, MMP9, RNASE6, HLA-DQB1, CHI3L1, CCL3, TXN
## PC_ 5
## Positive: PTGDS, LINC02244, MMP9, MGST1, LYZ, RGCC, IFI30, CHI3L1, CLU, GCLC
## GCHFR, CRABP2, FABP3, VAMP5, TMEM176B, HLA-DRB1, SYNGR1, CTSB, NCAPH, ISG15
## CD74, NCF1, SOD2, SAT1, FGL2, HSPA5, CDKN2A, MX1, CLEC12A, S100A10
## Negative: TYMS, PCLAF, TK1, MKI67, MYBL2, CENPM, CLSPN, RRM2, BIRC5, SHCBP1
## CEP55, DIAPH3, CDK1, CENPK, PRC1, CENPF, ESCO2, HELLS, NUSAP1, TOP2A
## NCAPG, FAM111B, ANLN, TPX2, CCNA2, KIF11, ASPM, CENPU, MAD2L1, HMMR
comb <- RunHarmony(comb,"sample")
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony 4/10
## Harmony converged after 4 iterations
DimHeatmap(comb, dims = 1:6, cells = 500, balanced = TRUE)
ElbowPlot(comb)
comb <- JackStraw(comb, num.replicate = 100)
comb <- FindNeighbors(comb, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
comb <- FindClusters(comb, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 28971
## Number of edges: 896422
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8952
## Number of communities: 12
## Elapsed time: 5 seconds
comb <- RunUMAP(comb, dims = 1:10)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 13:43:51 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:43:51 Read 28971 rows and found 10 numeric columns
## 13:43:51 Using Annoy for neighbor search, n_neighbors = 30
## 13:43:51 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:43:53 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea3bb493d7
## 13:43:53 Searching Annoy index using 1 thread, search_k = 3000
## 13:44:02 Annoy recall = 100%
## 13:44:03 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:44:05 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:44:05 Commencing optimization for 200 epochs, with 1186670 positive edges
## 13:44:13 Optimization finished
DimPlot(comb, reduction = "umap")
ADGRE1, CCR2, CD169, CX3CR1, CD206, CD163, LYVE1, CD9, TREM2
HLA-DP, HLA-DM, HLA-DOA, HLA-DOB, HLA-DQ, and HLA-DR.
message("macrophage markers")
## macrophage markers
FeaturePlot(comb, features = c("ADGRE1", "CCR2", "SIGLEC1", "CX3CR1", "MRC1", "CD163", "LYVE1", "CD9", "TREM2"))
DimPlot(comb, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
That’s pretty useless. Let’s use celldex pkg to annotate cells and get the macs.
ref <- celldex::MonacoImmuneData()
DefaultAssay(comb) <- "RNA"
comb2 <- as.SingleCellExperiment(comb)
lc <- logcounts(comb2)
pred_imm_broad <- SingleR(test=comb2, ref=ref,
labels=ref$label.main)
head(pred_imm_broad)
## DataFrame with 6 rows and 4 columns
## scores labels
## <matrix> <character>
## mdm_mock1|AAACGAATCACATACG 0.307826:0.325650:0.169064:... Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.311008:0.281528:0.189016:... Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.318344:0.277003:0.170895:... Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.214292:0.190767:0.121117:... Monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.290416:0.291088:0.155001:... Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.292140:0.281397:0.184535:... Monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.1786319 Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.0338547 Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1301308 Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.1116322 NA
## mdm_mock1|AAAGGTAAGCCATATC 0.1794308 Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.0952702 Monocytes
table(pred_imm_broad$pruned.labels)
##
## Basophils Dendritic cells Monocytes Neutrophils NK cells
## 4 657 24154 14 4
## Progenitors
## 2
cellmetadata$label <- pred_imm_broad$pruned.labels
pred_imm_fine <- SingleR(test=comb2, ref=ref,
labels=ref$label.fine)
head(pred_imm_fine)
## DataFrame with 6 rows and 4 columns
## scores labels
## <matrix> <character>
## mdm_mock1|AAACGAATCACATACG 0.1836411:0.485683:0.206442:... Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1961234:0.430705:0.226843:... Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1718613:0.442610:0.188544:... Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0943609:0.264157:0.120656:... Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1563980:0.415082:0.167965:... Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1857623:0.432131:0.207144:... Classical monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.0626269 Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1150706 Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.0651352 Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0648711 Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1073274 Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1521533 Classical monocytes
table(pred_imm_fine$pruned.labels)
##
## Classical monocytes Intermediate monocytes
## 21583 3569
## Low-density basophils Low-density neutrophils
## 2 21
## Myeloid dendritic cells Naive B cells
## 461 1
## Natural killer cells Non classical monocytes
## 3 64
## Non-Vd2 gd T cells Plasmacytoid dendritic cells
## 1 35
## Progenitor cells Terminal effector CD4 T cells
## 4 1
cellmetadata$finelabel <- pred_imm_fine$pruned.labels
col_pal <- c('#e31a1c', '#ff7f00', "#999900", '#cc00ff', '#1f78b4', '#fdbf6f',
'#33a02c', '#fb9a99', "#a6cee3", "#cc6699", "#b2df8a", "#99004d", "#66ff99",
"#669999", "#006600", "#9966ff", "#cc9900", "#e6ccff", "#3399ff", "#ff66cc",
"#ffcc66", "#003399")
annot_df <- data.frame(
barcodes = rownames(pred_imm_broad),
monaco_broad_annotation = pred_imm_broad$labels,
monaco_broad_pruned_labels = pred_imm_broad$pruned.labels,
monaco_fine_annotation = pred_imm_fine$labels,
monaco_fine_pruned_labels = pred_imm_fine$pruned.labels
)
meta_inf <- comb@meta.data
meta_inf$cell_barcode <- colnames(comb)
meta_inf <- meta_inf %>% dplyr::left_join(y = annot_df,
by = c("cell_barcode" = "barcodes"))
rownames(meta_inf) <- colnames(lc)
comb@meta.data <- meta_inf
DimPlot(comb, label=TRUE, group.by = "monaco_broad_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
DimPlot(comb, label=TRUE, group.by = "monaco_fine_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
mdmlist <- mylist[grep("mdm",names(mylist))]
comb1 <- do.call(cbind,mdmlist)
sce1 <- SingleCellExperiment(list(counts=comb1))
sce1
## class: SingleCellExperiment
## dim: 36603 12633
## metadata(0):
## assays(1): counts
## rownames(36603): gene-HIV1-1 gene-HIV1-2 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(12633): mdm_mock1|AAACGAATCACATACG mdm_mock1|AAACGCTCATCAGCGC
## ... mdm_bystander3|TTTCGATCACCTAAAC mdm_bystander3|TTTGACTGTATGTGTC
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
cellmetadata1 <- data.frame(colnames(comb1) ,sapply(strsplit(colnames(comb1),"\\|"),"[[",1))
colnames(cellmetadata1) <- c("cell","sample")
comb1 <- CreateSeuratObject(comb1, project = "mac", assay = "RNA", meta.data = cellmetadata1)
comb1 <- NormalizeData(comb1)
## Normalizing layer: counts
comb1 <- FindVariableFeatures(comb1, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
comb1 <- ScaleData(comb1)
## Centering and scaling data matrix
comb1 <- RunPCA(comb1, features = VariableFeatures(object = comb1))
## PC_ 1
## Positive: FABP4, MT-ATP8, FABP3, GAL, TXN, gene-HIV1-2, PRDX6, STMN1, NUPR1, S100A10
## UTS2, CYSTM1, LILRA1, GAPDH, SLAMF9, NMB, S100A9, GSTM3, LDHB, SEL1L2
## S100P, GDF15, APOC2, HAMP, UCHL1, MT1G, SLC35F1, COX5B, MMP12, OCIAD2
## Negative: ARL15, NEAT1, DPYD, EXOC4, FTX, ZEB2, PLXDC2, LRMDA, DOCK3, RASAL2
## ATG7, JMJD1C, TRIO, MALAT1, MAML2, FHIT, VPS13B, TMEM117, MITF, ELMO1
## TCF12, COP1, JARID2, DENND4C, ATXN1, MEF2A, ARHGAP15, DOCK4, FMNL2, MED13L
## PC_ 2
## Positive: TM4SF19, GPC4, CYSTM1, SPP1, ANO5, CD164, TXNRD1, MGST1, FNIP2, TXN
## C15orf48, ACAT2, FABP3, S100A10, CHI3L1, PRDX6, RETREG1, PSD3, RGCC, SCD
## COX5B, GDF15, TGM2, RGS20, CRABP2, MREG, CCL22, CSF1, GAL, HES2
## Negative: TGFBI, HLA-DPB1, HLA-DRA, MS4A6A, HLA-DPA1, HLA-DQB1, HLA-DQA1, CEBPD, MPEG1, C1QA
## ST8SIA4, FPR3, CD74, RNASE1, C1QC, CD163, TCF4, EPB41L3, CD14, LILRB2
## HLA-DRB5, ARL4C, PDE4B, MS4A7, SIPA1L1, AIF1, SEMA4A, FOS, HLA-DQA2, HLA-DOA
## PC_ 3
## Positive: GSN, CYP1B1, GCLC, S100A4, LPL, DUSP2, TIMP3, CALR, LINC02345, NCAPH
## CRIP1, OCSTAMP, CAMK1, ID2, ACTB, LHFPL2, FAIM2, STAC, MRC1, SPRED1
## IGSF6, RGCC, CA2, CCND2, AC020656.1, RCBTB2, MNDA, IL1RN, LINC02408, PLCL1
## Negative: NUPR1, SAT1, XIST, PLEKHA5, PSAT1, CCPG1, BTG1, G0S2, STMN1, GDF15
## CLGN, SUPV3L1, MARCKS, C5orf17, CCDC26, PHGDH, PRKCE, NAMPT, GRB10, S100P
## HES2, RETREG1, AC073359.2, BCAT1, REV3L, BEX2, MAML3, AL035446.2, NIBAN1, DTNB
## PC_ 4
## Positive: MT-ATP8, OSBP2, AC067751.1, TMEM117, TPRG1, AC092353.2, LINC02320, AC008591.1, LINC02408, KCNJ1
## NFATC3, PKD1L1, CT69, CPEB2, NOS1AP, PLCL1, AL136317.2, MICAL3, LY86-AS1, TTC28
## DNM3, ZNF366, XYLT1, DOCK3, ST5, ZFPM2, VWA8, MCU, AC093865.1, PDE4D
## Negative: IFI30, SRGN, CTSL, GAPDH, SAT1, BLVRB, LGMN, MARCKS, S100A10, ANXA1
## ACTG1, FUCA1, TXN, HSP90B1, PLIN2, TUBB, GLIPR1, HLA-C, CD164, ACTB
## CD74, TUBA1B, CD36, LDHB, TUBA1C, PLEK, HLA-DRB1, PRDX6, TMEM176A, HLA-DPA1
## PC_ 5
## Positive: SLC35F1, CCL7, gene-HIV1-1, CSF1, TPM4, CCL3, gene-HIV1-2, ACTB, MSMO1, SQLE
## ACTG1, TCTEX1D2, LINC01091, INSIG1, MADD, PCSK6, PHLDA1, SPOCD1, C3, MACC1
## ATP13A3, TNFRSF9, CCL2, IL1RN, DHCR24, FMN1, B4GALT5, TM4SF19, LIMA1, PALM2-AKAP2
## Negative: PTGDS, CLU, LINC02244, SYNGR1, BX664727.3, HLA-C, COX5B, NCAPH, RNASE6, MGST1
## RAB6B, CSRP2, FOLR2, CCPG1, RARRES1, LINC01010, AKR7A2, CPE, MEIKIN, AL136317.2
## PEBP4, CAMP, LYZ, RCBTB2, AC015660.2, VAMP5, MT1E, SOD2, MT2A, MT1X
comb1 <- RunHarmony(comb1,"sample")
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony 4/10
## Harmony 5/10
## Harmony converged after 5 iterations
DimHeatmap(comb1, dims = 1:6, cells = 500, balanced = TRUE)
ElbowPlot(comb1)
comb1 <- JackStraw(comb1, num.replicate = 100)
comb1 <- FindNeighbors(comb1, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
comb1 <- FindClusters(comb1, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 12633
## Number of edges: 400608
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8898
## Number of communities: 13
## Elapsed time: 1 seconds
comb1 <- RunUMAP(comb1, dims = 1:10)
## 13:49:44 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:49:44 Read 12633 rows and found 10 numeric columns
## 13:49:44 Using Annoy for neighbor search, n_neighbors = 30
## 13:49:44 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:49:45 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea25b5fc8d
## 13:49:45 Searching Annoy index using 1 thread, search_k = 3000
## 13:49:49 Annoy recall = 100%
## 13:49:50 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:49:51 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:49:52 Commencing optimization for 200 epochs, with 511998 positive edges
## 13:49:56 Optimization finished
DimPlot(comb1, reduction = "umap")
message("macrophage markers")
## macrophage markers
FeaturePlot(comb1, features = c("ADGRE1", "CCR2", "SIGLEC1", "CX3CR1", "MRC1", "CD163", "LYVE1", "CD9", "TREM2"))
DimPlot(comb1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ref <- celldex::MonacoImmuneData()
DefaultAssay(comb1) <- "RNA"
comb21 <- as.SingleCellExperiment(comb1)
lc1 <- logcounts(comb21)
pred_imm_broad1 <- SingleR(test=comb21, ref=ref,labels=ref$label.main)
head(pred_imm_broad1)
## DataFrame with 6 rows and 4 columns
## scores labels
## <matrix> <character>
## mdm_mock1|AAACGAATCACATACG 0.307826:0.325650:0.169064:... Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.311008:0.281528:0.189016:... Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.318344:0.277003:0.170895:... Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.214292:0.190767:0.121117:... Monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.290416:0.291088:0.155001:... Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.292140:0.281397:0.184535:... Monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.1786319 Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.0338547 Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1301308 Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.1116322 NA
## mdm_mock1|AAAGGTAAGCCATATC 0.1794308 Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.0952702 Monocytes
table(pred_imm_broad1$pruned.labels)
##
## Basophils Dendritic cells Monocytes Neutrophils
## 2 249 10174 2
cellmetadata1$label <- pred_imm_broad1$pruned.labels
pred_imm_fine1 <- SingleR(test=comb21, ref=ref, labels=ref$label.fine)
head(pred_imm_fine1)
## DataFrame with 6 rows and 4 columns
## scores labels
## <matrix> <character>
## mdm_mock1|AAACGAATCACATACG 0.1836411:0.485683:0.206442:... Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1961234:0.430705:0.226843:... Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1718613:0.442610:0.188544:... Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0943609:0.264157:0.120656:... Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1563980:0.415082:0.167965:... Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1857623:0.432131:0.207144:... Classical monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.0626269 Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1150706 Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.0651352 Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0648711 Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1073274 Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1521533 Classical monocytes
table(pred_imm_fine1$pruned.labels)
##
## Classical monocytes Intermediate monocytes
## 9388 1390
## Low-density neutrophils Myeloid dendritic cells
## 3 159
## Non classical monocytes Plasmacytoid dendritic cells
## 23 15
cellmetadata1$finelabel <- pred_imm_fine1$pruned.labels
col_pal <- c('#e31a1c', '#ff7f00', "#999900", '#cc00ff', '#1f78b4', '#fdbf6f',
'#33a02c', '#fb9a99', "#a6cee3", "#cc6699", "#b2df8a", "#99004d", "#66ff99",
"#669999", "#006600", "#9966ff", "#cc9900", "#e6ccff", "#3399ff", "#ff66cc",
"#ffcc66", "#003399")
annot_df1 <- data.frame(
barcodes = rownames(pred_imm_broad1),
monaco_broad_annotation = pred_imm_broad1$labels,
monaco_broad_pruned_labels = pred_imm_broad1$pruned.labels,
monaco_fine_annotation = pred_imm_fine1$labels,
monaco_fine_pruned_labels = pred_imm_fine1$pruned.labels)
meta_inf1 <- comb1@meta.data
meta_inf1$cell_barcode <- colnames(comb1)
meta_inf1 <- meta_inf1 %>% dplyr::left_join(y = annot_df1, by = c("cell_barcode" = "barcodes"))
rownames(meta_inf1) <- colnames(lc1)
comb1@meta.data <- meta_inf1
DimPlot(comb1, label=TRUE, group.by = "monaco_broad_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
DimPlot(comb1, label=TRUE, group.by = "monaco_fine_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
message("extract mono")
## extract mono
mono <- comb1[,which(meta_inf1$monaco_broad_annotation == "Monocytes")]
mono_metainf1 <- meta_inf1[which(meta_inf1$monaco_broad_annotation == "Monocytes"),]
mono_metainf1 <- mono_metainf1[grep("monocytes",mono_metainf1$monaco_fine_pruned_labels),]
mono <- mono[,which(colnames(mono) %in% rownames(mono_metainf1))]
mono <- FindVariableFeatures(mono, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
mono <- RunPCA(mono, features = VariableFeatures(object = mono))
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: MT-ND5, MT-ND2, MT-ND1, AL358779.1,
## TSPAN8, FAH, PFN1, TCEA1, GLUL, EPB41L1, DDIT3, EPHA6, TUBB4B, TFRC, CNIH3,
## RABGEF1, AL669970.1, ECE2, SNHG5, MCTP1, RHOF, FRMD4A, FARP1, RERE, MTRNR2L1,
## RAP1GDS1, GPAT3, SETD2, HSP90AA1, STRBP, BCL2A1, EDIL3, NPTX1, RCAN1, RREB1,
## CFD, DZIP1, PEAK1, ACP5, BIN2, CSTA, PPARGC1B, ATP6V1D.
## PC_ 1
## Positive: S100A10, GAPDH, TXN, FABP3, PRDX6, CYSTM1, FABP4, LDHB, COX5B, C15orf48
## BLVRB, S100A9, PSME2, SLAMF9, STMN1, GDF15, SPP1, NMB, LILRA1, GAL
## CTSL, TUBA1B, HAMP, UCHL1, CARD16, ANXA1, NUPR1, LILRA2, TUBA1A, GSTM4
## Negative: ARL15, FTX, DPYD, EXOC4, VPS13B, LRMDA, DOCK3, NEAT1, ELMO1, RASAL2
## ATG7, JMJD1C, TRIO, TMEM117, ZEB2, MAML2, FHIT, DOCK4, MALAT1, TCF12
## COP1, SPIDR, ATP9B, ARHGAP15, MED13L, ZSWIM6, MBD5, RAD51B, TPRG1, TBC1D5
## PC_ 2
## Positive: TM4SF19, ANO5, GPC4, FNIP2, SPP1, TXNRD1, PSD3, CYSTM1, RETREG1, RGS20
## CD164, TCTEX1D2, MGST1, ACAT2, CCDC26, CCL22, MREG, SLC28A3, FABP3, SCD
## LINC01135, C15orf48, TXN, GDF15, GAL, HES2, CHI3L1, NIBAN1, COL8A2, RGCC
## Negative: TGFBI, HLA-DPB1, HLA-DRA, HLA-DPA1, CD74, HLA-DQB1, MS4A6A, HLA-DQA1, C1QA, CEBPD
## C1QC, MPEG1, FPR3, AIF1, CD163, MS4A7, CD14, HLA-DRB5, ST8SIA4, LILRB2
## RNASE1, TCF4, HLA-DRB1, FOS, EPB41L3, ARL4C, HLA-DQA2, TIMP1, SEMA4A, PDE4B
## PC_ 3
## Positive: NUPR1, SAT1, BTG1, CCPG1, MARCKS, PSAT1, G0S2, XIST, STMN1, PLEKHA5
## GDF15, SUPV3L1, CLGN, NAMPT, PHGDH, ADAMDEC1, HES2, S100P, BCAT1, CXCR4
## SDS, C5orf17, SLAMF9, REV3L, RETREG1, CCDC26, GRB10, PRKCE, BEX2, CD14
## Negative: GCLC, GSN, CYP1B1, DUSP2, LINC02345, S100A4, TIMP3, LPL, NCAPH, OCSTAMP
## CRIP1, STAC, CAMK1, FAIM2, CALR, SPRED1, ID2, LHFPL2, MRC1, LINC02408
## PLCL1, RCBTB2, HIVEP3, ACTB, CA2, RNF128, AC020656.1, CCND2, IGSF6, RGCC
## PC_ 4
## Positive: MT-ATP8, PTGDS, AC008591.1, LY86-AS1, AL136317.2, RCBTB2, KCNJ1, PKD1L1, BX664727.3, AC067751.1
## EPHB1, LINC01010, MT1G, CLU, TMEM117, PDE4D, OSBP2, LINC02320, LINC02244, AC015660.2
## XYLT1, MT1H, CT69, KLRD1, LINC02408, CC2D2B, MT1X, AP000331.1, SPON2, KCNMA1
## Negative: ACTG1, ACTB, TPM4, TUBA1B, GAPDH, TUBB, GLIPR1, CSF1, CD36, CCL3
## CTSL, PLEK, LGMN, LDHA, MGLL, HSP90B1, TUBA1C, MARCKS, CALR, INSIG1
## MSMO1, SLC35F1, C15orf48, ANXA1, PHLDA1, BLVRB, CD164, TMEM176A, DUSP6, IFI30
## PC_ 5
## Positive: CSF1, ACTB, TCTEX1D2, MSMO1, gene-HIV1-1, CCL3, gene-HIV1-2, TPM4, ACTG1, TM4SF19
## MADD, SPOCD1, CCL7, SLC35F1, PCSK6, UGCG, PHLDA1, PPARG, SPRED1, ATP13A3
## ANK2, IL1RN, CD36, CBLB, SQLE, TM4SF19-AS1, B4GALT5, TNFRSF9, NRP1, TIMP3
## Negative: TYMS, PCLAF, TK1, MKI67, MYBL2, CEP55, BIRC5, CENPM, DIAPH3, CENPK
## RRM2, KIF11, PRC1, CDK1, SHCBP1, CLSPN, CENPF, ANLN, RAD51AP1, NUSAP1
## CENPU, KIF4A, ESCO2, TOP2A, TPX2, ASPM, CIT, MAD2L1, FAM111B, KNL1
DimHeatmap(mono, dims = 1:2, cells = 500, balanced = TRUE)
DimHeatmap(mono, dims = 3:4, cells = 500, balanced = TRUE)
ElbowPlot(mono)
mono <- FindNeighbors(mono, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
mono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)
## 13:50:52 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:50:52 Read 10758 rows and found 4 numeric columns
## 13:50:52 Using Annoy for neighbor search, n_neighbors = 30
## 13:50:52 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:50:53 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea6a14858b
## 13:50:53 Searching Annoy index using 1 thread, search_k = 3000
## 13:50:57 Annoy recall = 100%
## 13:50:58 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:51:01 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:51:01 Commencing optimization for 200 epochs, with 364062 positive edges
## 13:51:04 Optimization finished
DimPlot(mono, reduction = "umap", label=TRUE)
DimPlot(mono, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE)
DimPlot(mono, group.by="sample" , reduction = "umap", label=TRUE)
alvlist <- mylist[grep("alv",names(mylist))]
comb1 <- do.call(cbind,alvlist)
sce1 <- SingleCellExperiment(list(counts=comb1))
sce1
## class: SingleCellExperiment
## dim: 36603 13240
## metadata(0):
## assays(1): counts
## rownames(36603): gene-HIV1-1 gene-HIV1-2 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(13240): alv_mock1|AAACCCAGTGCTGCAC alv_mock1|AAAGAACCATTGAGGG
## ... alv_bystander4|TTTCGATGTGAGCAGT alv_bystander4|TTTGATCTCGGCTTGG
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
cellmetadata1 <- data.frame(colnames(comb1) ,sapply(strsplit(colnames(comb1),"\\|"),"[[",1))
colnames(cellmetadata1) <- c("cell","sample")
comb1 <- CreateSeuratObject(comb1, project = "mac", assay = "RNA", meta.data = cellmetadata1)
comb1 <- NormalizeData(comb1)
## Normalizing layer: counts
comb1 <- FindVariableFeatures(comb1, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
comb1 <- ScaleData(comb1)
## Centering and scaling data matrix
comb1 <- RunPCA(comb1, features = VariableFeatures(object = comb1))
## PC_ 1
## Positive: FABP4, gene-HIV1-2, NUPR1, TIMP1, S100A8, STMN1, VSIG4, CAMP, RBP1, LILRA1
## ORM1, PRG2, HLA-DOA, CD79A, CCL5, PDPN, VMO1, MDK, LGALS2, UTS2
## AC078850.1, AC026369.3, RBP4, MS4A6E, PTGES, HCFC1-AS1, PPBP, CXCL10, FXYD2, LINC01467
## Negative: DOCK3, RASAL2, ARL15, MALAT1, NEAT1, ASAP1, PLXDC2, LRMDA, TMEM117, DPYD
## EXOC4, MITF, ATG7, FTX, FRMD4B, JMJD1C, FMNL2, ZNF438, UBE2E2, FHIT
## MAML2, LPP, CHD9, DENND4C, TPRG1, ZEB2, LRCH1, RAPGEF1, ELMO1, MEF2A
## PC_ 2
## Positive: CYSTM1, CSTB, TM4SF19, GAL, ATP6V1D, FDX1, GM2A, CCL22, TGM2, ACAT2
## LGALS1, SCD, CD164, ELOC, FAH, PRDX6, CSF1, RHOF, GAPLINC, FCMR
## CYTOR, C15orf48, UPP1, GOLGA7B, CIR1, TNFRSF12A, NCAPH, SNHG32, DCSTAMP, BCAT1
## Negative: TRPS1, SAMSN1, SLC8A1, PDGFC, AL356124.1, TGFBI, CAMK1D, SPRED1, SLC9A9, ARHGAP15
## RTN1, RCBTB2, ATP8B4, MRC1, FCHSD2, CD14, XYLT1, FOS, HLA-DPB1, ME1
## C1QC, TCF7L2, AIF1, MS4A6A, MAML3, VAMP5, C1QA, DLEU2, TCF12, DIAPH2
## PC_ 3
## Positive: KCNJ1, C2orf92, LINC02320, AC067751.1, AC015660.2, NCAPH, RGS20, LINC01800, AL136317.2, CCL22
## XIST, BX664727.3, LINC01010, AP000331.1, CT69, MICAL3, GDA, AC012668.3, ANO5, TRIM54
## AC008591.1, NOS1AP, AC092353.2, PDE4D, TCTEX1D2, LY86-AS1, TDRD3, MIR646HG, OSBP2, SLC25A48
## Negative: AIF1, IFI30, CTSZ, MARCKS, CD74, HLA-DRA, HLA-DPA1, LGMN, HLA-DRB1, C1QA
## HLA-DPB1, CTSL, CTSB, ACTB, FPR3, C1QC, BLVRB, ACTG1, HLA-DQB1, FUCA1
## TUBA1B, CTSC, HLA-DQA1, GAPDH, FCGR3A, IL4I1, TUBB, C1QB, FBP1, CD163
## PC_ 4
## Positive: FMN1, MARCO, SLC11A1, XIST, SNCA, HAMP, S100A9, FRMD4A, CCDC26, BCAT1
## ME3, CLMP, SASH1, TDRD9, TTC39B, PHACTR1, SLC16A10, DLEU1, UGCG, AC023282.1
## PLEKHA5, SEL1L2, SLA, AL035446.2, MSR1, OSBP2, GYPC, DNAJC6, AP005262.2, B4GALT5
## Negative: PTGDS, GCLC, CRIP1, CLU, MX1, MMP9, TMEM176A, ALOX5AP, TUBA1A, FGL2
## PLEK, LINC02244, VAMP5, IFIT3, IFI6, MMP7, KCNMA1, DUSP2, RGCC, MGST1
## MX2, RNF128, SYNGR1, ISG15, NCAPH, IFIT1, C1QTNF4, MIF, FAIM2, HLA-DRB1
## PC_ 5
## Positive: CFD, CTSZ, IFI30, SAT1, CTSB, FUCA1, LYZ, GCHFR, CD74, AC020656.1
## HLA-DRB1, CSTB, IFI6, CHI3L1, ATP1B1, C15orf48, GPX3, RARRES1, TFRC, HLA-DQB1
## HLA-DRB5, MS4A7, GM2A, RGCC, HLA-DRA, MMP9, CA2, HLA-DPB1, SELENOP, RASSF4
## Negative: TYMS, PCLAF, TK1, CLSPN, MKI67, CENPM, RRM2, DIAPH3, MYBL2, BIRC5
## CDK1, SHCBP1, ESCO2, CEP55, NCAPG, CENPK, FAM111B, NUSAP1, TOP2A, CENPF
## HELLS, ATAD2, PRC1, CCNA2, HMMR, ASPM, TPX2, KIF11, MCM10, DTL
comb1 <- RunHarmony(comb1,"sample")
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony 4/10
## Harmony converged after 4 iterations
DimHeatmap(comb1, dims = 1:6, cells = 500, balanced = TRUE)
ElbowPlot(comb1)
comb1 <- JackStraw(comb1, num.replicate = 100)
comb1 <- FindNeighbors(comb1, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
comb1 <- FindClusters(comb1, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 13240
## Number of edges: 416149
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8777
## Number of communities: 11
## Elapsed time: 1 seconds
comb1 <- RunUMAP(comb1, dims = 1:10)
## 13:55:41 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:55:41 Read 13240 rows and found 10 numeric columns
## 13:55:41 Using Annoy for neighbor search, n_neighbors = 30
## 13:55:41 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:55:42 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea5e26fe47
## 13:55:42 Searching Annoy index using 1 thread, search_k = 3000
## 13:55:46 Annoy recall = 100%
## 13:55:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:55:50 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:55:50 Commencing optimization for 200 epochs, with 539610 positive edges
## 13:55:54 Optimization finished
DimPlot(comb1, reduction = "umap")
message("macrophage markers")
## macrophage markers
FeaturePlot(comb1, features = c("ADGRE1", "CCR2", "SIGLEC1", "CX3CR1", "MRC1", "CD163", "LYVE1", "CD9", "TREM2"))
DimPlot(comb1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ref <- celldex::MonacoImmuneData()
DefaultAssay(comb1) <- "RNA"
comb21 <- as.SingleCellExperiment(comb1)
lc1 <- logcounts(comb21)
pred_imm_broad1 <- SingleR(test=comb21, ref=ref,labels=ref$label.main)
head(pred_imm_broad1)
## DataFrame with 6 rows and 4 columns
## scores labels
## <matrix> <character>
## alv_mock1|AAACCCAGTGCTGCAC 0.273681:0.2723954:0.1612672:... Monocytes
## alv_mock1|AAAGAACCATTGAGGG 0.125532:0.0955108:0.0727378:... Dendritic cells
## alv_mock1|AAAGAACTCTCATGCC 0.093777:0.0888391:0.0837904:... Dendritic cells
## alv_mock1|AAAGGATAGCATGAAT 0.318637:0.3074635:0.1919116:... Monocytes
## alv_mock1|AAAGGATAGTCAGGGT 0.294499:0.2756627:0.1829148:... Monocytes
## alv_mock1|AAAGGATAGTTCCGGC 0.293979:0.2941562:0.1826853:... Monocytes
## delta.next pruned.labels
## <numeric> <character>
## alv_mock1|AAACCCAGTGCTGCAC 0.13448850 Monocytes
## alv_mock1|AAAGAACCATTGAGGG 0.11133884 Dendritic cells
## alv_mock1|AAAGAACTCTCATGCC 0.00571364 Dendritic cells
## alv_mock1|AAAGGATAGCATGAAT 0.10466784 Monocytes
## alv_mock1|AAAGGATAGTCAGGGT 0.13985529 Monocytes
## alv_mock1|AAAGGATAGTTCCGGC 0.12840710 Monocytes
table(pred_imm_broad1$pruned.labels)
##
## Basophils Dendritic cells Monocytes Neutrophils NK cells
## 2 329 11305 10 3
## Progenitors
## 1
cellmetadata1$label <- pred_imm_broad1$pruned.labels
pred_imm_fine1 <- SingleR(test=comb21, ref=ref, labels=ref$label.fine)
head(pred_imm_fine1)
## DataFrame with 6 rows and 4 columns
## scores
## <matrix>
## alv_mock1|AAACCCAGTGCTGCAC 0.1630294:0.398305:0.1767129:...
## alv_mock1|AAAGAACCATTGAGGG 0.0671221:0.151726:0.0730234:...
## alv_mock1|AAAGAACTCTCATGCC 0.0614206:0.101479:0.0761397:...
## alv_mock1|AAAGGATAGCATGAAT 0.1803723:0.435098:0.2088530:...
## alv_mock1|AAAGGATAGTCAGGGT 0.1703657:0.388786:0.1868922:...
## alv_mock1|AAAGGATAGTTCCGGC 0.1664892:0.422481:0.1894544:...
## labels delta.next
## <character> <numeric>
## alv_mock1|AAACCCAGTGCTGCAC Classical monocytes 1.43376e-01
## alv_mock1|AAAGAACCATTGAGGG Myeloid dendritic ce.. 2.22045e-16
## alv_mock1|AAAGAACTCTCATGCC Myeloid dendritic ce.. 5.33208e-03
## alv_mock1|AAAGGATAGCATGAAT Classical monocytes 1.21257e-01
## alv_mock1|AAAGGATAGTCAGGGT Classical monocytes 5.02055e-02
## alv_mock1|AAAGGATAGTTCCGGC Classical monocytes 9.94518e-02
## pruned.labels
## <character>
## alv_mock1|AAACCCAGTGCTGCAC Classical monocytes
## alv_mock1|AAAGAACCATTGAGGG Myeloid dendritic ce..
## alv_mock1|AAAGAACTCTCATGCC Myeloid dendritic ce..
## alv_mock1|AAAGGATAGCATGAAT Classical monocytes
## alv_mock1|AAAGGATAGTCAGGGT Classical monocytes
## alv_mock1|AAAGGATAGTTCCGGC Classical monocytes
table(pred_imm_fine1$pruned.labels)
##
## Classical monocytes Intermediate monocytes
## 9996 1550
## Low-density basophils Low-density neutrophils
## 1 14
## Myeloid dendritic cells Natural killer cells
## 246 3
## Non classical monocytes Plasmacytoid dendritic cells
## 32 9
## Progenitor cells Terminal effector CD4 T cells
## 4 1
cellmetadata1$finelabel <- pred_imm_fine1$pruned.labels
col_pal <- c('#e31a1c', '#ff7f00', "#999900", '#cc00ff', '#1f78b4', '#fdbf6f',
'#33a02c', '#fb9a99', "#a6cee3", "#cc6699", "#b2df8a", "#99004d", "#66ff99",
"#669999", "#006600", "#9966ff", "#cc9900", "#e6ccff", "#3399ff", "#ff66cc",
"#ffcc66", "#003399")
annot_df1 <- data.frame(
barcodes = rownames(pred_imm_broad1),
monaco_broad_annotation = pred_imm_broad1$labels,
monaco_broad_pruned_labels = pred_imm_broad1$pruned.labels,
monaco_fine_annotation = pred_imm_fine1$labels,
monaco_fine_pruned_labels = pred_imm_fine1$pruned.labels)
meta_inf1 <- comb1@meta.data
meta_inf1$cell_barcode <- colnames(comb1)
meta_inf1 <- meta_inf1 %>% dplyr::left_join(y = annot_df1, by = c("cell_barcode" = "barcodes"))
rownames(meta_inf1) <- colnames(lc1)
comb1@meta.data <- meta_inf1
DimPlot(comb1, label=TRUE, group.by = "monaco_broad_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
DimPlot(comb1, label=TRUE, group.by = "monaco_fine_annotation", reduction = "umap",
cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")
message("extract mono")
## extract mono
mono <- comb1[,which(meta_inf1$monaco_broad_annotation == "Monocytes")]
mono_metainf1 <- meta_inf1[which(meta_inf1$monaco_broad_annotation == "Monocytes"),]
mono_metainf1 <- mono_metainf1[grep("monocytes",mono_metainf1$monaco_fine_pruned_labels),]
mono <- mono[,which(colnames(mono) %in% rownames(mono_metainf1))]
mono <- FindVariableFeatures(mono, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
mono <- RunPCA(mono, features = VariableFeatures(object = mono))
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: IGSF6, RPS20, TRIO, MOSPD1, DOP1B,
## SNX10, AC008443.9, SDC2, SCFD1, MCTP1, CRABP2, QKI, HTRA4, S100A6, ALDH1A1,
## ASPH, NABP1, C11orf45, MGAT5, ARHGEF3, CDKAL1, ZFAND3, PTPRJ, LGALS3, HMOX1,
## AC092422.1, RP1L1, UBASH3B, PPP1R12B, RALGAPA2, AEBP2, C1orf21, DPH6, PLA2G7,
## DOCK2, LPAR1, AFF1, ABCA1, MB21D2, ATP6V1H, PITPNC1, ATP6V0D2, DMXL1, MARCH1,
## RNLS, ZPLD1, SLC22A15, ZNF609, SPIRE1, CALM3, CHAC1, SH3PXD2B, SLC4A1.
## PC_ 1
## Positive: GAPDH, LGALS1, MIF, CSTB, BLVRB, TUBA1B, IFI30, PRDX6, FABP4, TMEM176A
## TXN, FAH, TUBA1A, CD74, EMP3, HLA-DPA1, TIMP1, ACTB, S100A9, NUPR1
## MMP9, HLA-DRB1, CFD, C15orf48, HLA-DRA, ELOC, AIF1, H2AFZ, PTGDS, LILRA1
## Negative: DOCK3, ARL15, MALAT1, TMEM117, RASAL2, FTX, EXOC4, FHIT, LRMDA, DPYD
## TPRG1, VPS13B, JMJD1C, NEAT1, ELMO1, ATG7, MAML2, ZNF438, COP1, MITF
## TTC28, FMNL2, VTI1A, MED13L, ERC1, ATP9B, ASAP1, DENND4C, SPIDR, RAD51B
## PC_ 2
## Positive: GAL, CCL22, TM4SF19, FDX1, CYSTM1, ATP6V1D, TGM2, ACAT2, SCD, CIR1
## GM2A, RHOF, CSF1, NCAPH, IARS, GOLGA7B, DUSP13, BCAT1, CSTB, CSRP2
## FAH, SNHG32, SLC20A1, FCMR, CD164, GAPLINC, TCTEX1D2, EPB41L1, SH3BP5, NMB
## Negative: HLA-DPA1, HLA-DPB1, HLA-DRA, AIF1, CD74, C1QA, TGFBI, SAMSN1, FOS, C1QC
## CD14, MRC1, CTSC, SLC8A1, VAMP5, MS4A6A, HLA-DRB1, TRPS1, SLC9A9, FPR3
## CLEC4A, PDGFC, LYZ, SLCO2B1, FGL2, C1QB, RTN1, CAMK1D, CLEC7A, SPRED1
## PC_ 3
## Positive: MARCKS, LGMN, TPM4, CTSL, AIF1, FPR3, HLA-DQA1, CTSZ, HAMP, CD36
## ACTG1, CD163, C1QC, HLA-DQB1, GYPC, CCL3, IL18BP, BLVRB, STMN1, SMS
## FMN1, MARCO, C1QA, FUCA1, PLAU, IFI30, FABP4, FCGR3A, COLEC12, S100A9
## Negative: PTGDS, LINC02244, CLU, LINC01010, NCAPH, CRIP1, AC015660.2, LINC01800, SYNGR1, AC067751.1
## GCLC, KCNMA1, RCBTB2, C1QTNF4, KCNJ1, DUSP2, SPON2, NRCAM, SERTAD2, C2orf92
## MGST1, FGL2, TRIM54, CT69, RNF128, MX1, NCF1, NFATC3, RGS20, IFIT3
## PC_ 4
## Positive: XIST, JAML, HLA-DRB5, QPCT, SLC11A1, GCHFR, MARCO, PAX8-AS1, S100A9, CCDC26
## AC020656.1, SASH1, MSR1, MS4A7, AL136317.2, GPX3, C22orf34, AL035446.2, TMTC1, FRMD4A
## OSBP2, RARRES1, MIR646HG, TDRD9, ST6GAL1, DLEU1, SLC6A16, NMB, GPRIN3, CLMP
## Negative: TYMS, PCLAF, CLSPN, TK1, CENPM, RRM2, DIAPH3, MYBL2, MKI67, SHCBP1
## CDK1, BIRC5, CEP55, HELLS, ESCO2, FAM111B, NCAPG, NUSAP1, HMMR, TOP2A
## KIF11, ATAD2, CENPF, PRC1, CENPK, TPX2, CENPU, MCM10, CCNA2, ASPM
## PC_ 5
## Positive: TYMS, MKI67, PCLAF, TK1, JAML, RRM2, CDK1, CLSPN, DIAPH3, CENPK
## CENPM, MYBL2, XIST, BIRC5, SHCBP1, ESCO2, RAD51AP1, TOP2A, NCAPG, CENPF
## CCNA2, FAM111B, ATAD2, CEP55, NUSAP1, HLA-DRB5, ASPM, ANLN, PRC1, DTL
## Negative: TMEM176A, PLEK, LINC01091, MIF, IL1RN, ACTG1, TUBA1A, CYTOR, MYL9, ACTB
## LINC00278, MMP19, MGLL, PHLDA1, ADAMDEC1, TTTY14, CSF1, AC006001.2, GLIPR1, TEX14
## PRDX6, RABGEF1, HLA-DPA1, LCP1, PPM1H, CCND1, HDAC9, DPYSL3, TIMP3, HLA-DRA
DimHeatmap(mono, dims = 1:2, cells = 500, balanced = TRUE)
DimHeatmap(mono, dims = 3:4, cells = 500, balanced = TRUE)
ElbowPlot(mono)
mono <- FindNeighbors(mono, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
mono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)
## 13:56:49 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:56:49 Read 11565 rows and found 4 numeric columns
## 13:56:49 Using Annoy for neighbor search, n_neighbors = 30
## 13:56:49 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:56:50 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea29100cdd
## 13:56:50 Searching Annoy index using 1 thread, search_k = 3000
## 13:56:54 Annoy recall = 100%
## 13:56:55 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:56:57 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:56:57 Commencing optimization for 200 epochs, with 384128 positive edges
## 13:57:01 Optimization finished
DimPlot(mono, reduction = "umap", label=TRUE)
DimPlot(mono, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE)
DimPlot(mono, group.by="sample" , reduction = "umap", label=TRUE)
mono <- comb[,which(meta_inf$monaco_broad_annotation == "Monocytes")]
mono_metainf <- meta_inf[which(meta_inf$monaco_broad_annotation == "Monocytes"),]
mono_metainf1 <- mono_metainf[grep("monocytes",mono_metainf$monaco_fine_pruned_labels),]
mono <- mono[,which(colnames(mono) %in% rownames(mono_metainf))]
mono <- FindVariableFeatures(mono, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
mono <- RunPCA(mono, features = VariableFeatures(object = mono))
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: LINC02739, ARHGAP10, AL157702.2, DUSP1,
## RPS6KA2, CEP85L, CAMK1, OTOA, ANKRD44, RSPO3, TMCC1.
## PC_ 1
## Positive: MT-ATP8, MTRNR2L1, FABP4, FABP3, GAPDH, HAMP, S100A9, AC105402.3, C1orf56, PVALB
## AL136097.2, TUBA1B, CARD16, FAH, BIRC7, C1QC, UTS2, gene-HIV1-2, MT2A, NUPR1
## BLVRB, BEX3, LILRA1, GCHFR, AC025154.2, MT1G, TXN, GSTP1, PRDX6, FXYD2
## Negative: ARL15, DOCK3, NEAT1, EXOC4, FTX, MALAT1, RASAL2, LRMDA, DPYD, JMJD1C
## TMEM117, FHIT, PLXDC2, VPS13B, MITF, ZEB2, DENND4C, MAML2, ATG7, ZNF438
## TPRG1, COP1, FRMD4B, ELMO1, JARID2, ATXN1, FMNL2, TCF12, ERC1, ARID1B
## PC_ 2
## Positive: CYSTM1, FAH, GDF15, PSAT1, CD164, FDX1, ATP6V1D, CCPG1, BCAT1, SAT1
## SLAMF9, TXN, PHGDH, HES2, CSTA, HEBP2, NUPR1, STMN1, TCEA1, CTSL
## GARS, MARCO, RETREG1, S100A10, RILPL2, CLGN, SNHG5, SCD, RHOQ, SNHG32
## Negative: SPRED1, RCBTB2, KCNMA1, GCLC, HLA-DRB1, HLA-DRA, FGL2, CD74, HLA-DPA1, MRC1
## SLCO2B1, AC020656.1, HLA-DPB1, C1QA, RTN1, LINC02345, STAC, MERTK, PDGFC, VAMP5
## CCDC102B, DPYD, RNF128, MX1, ATP8B4, TGFBI, CRIP1, NFATC3, HLA-DRB5, ATG7
## PC_ 3
## Positive: NCAPH, CRABP2, TM4SF19, GAL, DUSP2, RGCC, CHI3L1, CCL22, TMEM114, ACAT2
## TRIM54, RGS20, LINC01010, LINC02244, MGST1, TCTEX1D2, GPC4, NUMB, CCND1, SYNGR1
## CLU, SERTAD2, AC092353.2, POLE4, SLC20A1, IL1RN, ZNF366, AC005280.2, OCSTAMP, GCLC
## Negative: MARCKS, CD14, TGFBI, BTG1, MS4A6A, HLA-DQA1, FPR3, CTSC, CD163, MPEG1
## MEF2C, AIF1, TIMP1, HLA-DPB1, C1QC, HLA-DQB1, SELENOP, RNASE1, HLA-DQA2, NAMPT
## TCF4, EPB41L3, ARL4C, NUPR1, ALDH2, SAMSN1, ZNF331, HIF1A, MS4A7, SEMA4A
## PC_ 4
## Positive: MT-ATP8, MTRNR2L1, XIST, PDE4D, AL035446.2, AC073359.2, AL365295.1, CCDC26, EPHB1, LY86-AS1
## CLVS1, RAD51B, PRKCE, C5orf17, LINC02320, FTX, DMGDH, LINC01473, LINC01135, FHIT
## AC008591.1, PLEKHA5, MIR646HG, BCAS3, PKD1L1, AF117829.1, LINC02649, ZBTB20, STX18-AS1, AC079142.1
## Negative: HLA-DRB1, CTSB, TUBA1B, HSP90B1, CD74, GAPDH, ACTG1, HLA-DPA1, HLA-DRA, IFI30
## RGCC, HSPA5, AIF1, C15orf48, C1QA, LGMN, TUBB, CYP1B1, HLA-DPB1, FBP1
## CA2, CCL3, TUBA1A, HLA-DQB1, S100A4, RNASE6, LYZ, TXN, MMP9, TNFSF13
## PC_ 5
## Positive: RGCC, MMP9, IFI30, CTSB, PTGDS, LINC02244, MGST1, HLA-DRB1, CHI3L1, LYZ
## CD74, C15orf48, HSPA5, FABP3, GCLC, GCHFR, CRABP2, HLA-DPA1, TXN, HLA-DRA
## ISG15, S100A10, NCAPH, TMEM176B, CDKN2A, RALA, HSP90B1, SAT1, CLU, VAMP5
## Negative: TYMS, PCLAF, TK1, MKI67, MYBL2, CENPM, RRM2, BIRC5, CLSPN, CEP55
## SHCBP1, DIAPH3, CDK1, CENPK, PRC1, CENPF, ESCO2, NUSAP1, TOP2A, HELLS
## NCAPG, ANLN, FAM111B, TPX2, CCNA2, KIF11, ASPM, CENPU, MAD2L1, HMMR
DimHeatmap(mono, dims = 1:2, cells = 500, balanced = TRUE)
DimHeatmap(mono, dims = 3:4, cells = 500, balanced = TRUE)
ElbowPlot(mono)
mono <- FindNeighbors(mono, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
mono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)
## 13:57:54 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:57:54 Read 28289 rows and found 4 numeric columns
## 13:57:54 Using Annoy for neighbor search, n_neighbors = 30
## 13:57:54 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:57:56 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea74d9bf5e
## 13:57:57 Searching Annoy index using 1 thread, search_k = 3000
## 13:58:05 Annoy recall = 100%
## 13:58:06 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:58:08 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:58:08 Commencing optimization for 200 epochs, with 936166 positive edges
## 13:58:15 Optimization finished
DimPlot(mono, reduction = "umap", label=TRUE)
DimPlot(mono, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE)
DimPlot(mono, group.by="sample" , reduction = "umap", label=TRUE)
Most cells are classified as monocytes.
cc <- table(meta_inf[,c("sample","monaco_broad_annotation")])
cc %>% kbl(caption="cell counts") %>% kable_paper("hover", full_width = F)
Basophils | Dendritic cells | Monocytes | Neutrophils | NK cells | Progenitors | |
---|---|---|---|---|---|---|
alv_active1 | 1 | 34 | 908 | 2 | 0 | 0 |
alv_active2 | 0 | 8 | 613 | 0 | 0 | 0 |
alv_active3 | 1 | 28 | 684 | 0 | 0 | 0 |
alv_bystander1 | 0 | 29 | 1028 | 2 | 0 | 0 |
alv_bystander2 | 0 | 20 | 1100 | 0 | 1 | 1 |
alv_bystander3 | 0 | 12 | 491 | 1 | 0 | 0 |
alv_bystander4 | 0 | 20 | 558 | 3 | 0 | 0 |
alv_latent1 | 0 | 12 | 740 | 0 | 0 | 0 |
alv_latent2 | 0 | 30 | 1905 | 1 | 0 | 0 |
alv_latent3 | 1 | 26 | 1804 | 0 | 0 | 0 |
alv_latent4 | 0 | 52 | 1823 | 0 | 0 | 0 |
alv_mock1 | 0 | 18 | 956 | 0 | 0 | 0 |
alv_mock2 | 0 | 19 | 785 | 0 | 0 | 0 |
alv_mock3 | 0 | 47 | 1304 | 1 | 2 | 0 |
mdm_active1 | 0 | 17 | 863 | 0 | 0 | 0 |
mdm_active2 | 0 | 3 | 476 | 0 | 0 | 0 |
mdm_active3 | 0 | 13 | 408 | 0 | 0 | 0 |
mdm_active4 | 0 | 4 | 432 | 0 | 0 | 0 |
mdm_bystander1 | 0 | 48 | 1446 | 0 | 0 | 0 |
mdm_bystander2 | 0 | 31 | 1303 | 0 | 0 | 0 |
mdm_bystander3 | 0 | 47 | 490 | 3 | 0 | 0 |
mdm_latent1 | 1 | 17 | 1154 | 0 | 0 | 0 |
mdm_latent2 | 0 | 11 | 1240 | 0 | 0 | 0 |
mdm_latent3 | 0 | 11 | 337 | 0 | 0 | 0 |
mdm_mock1 | 0 | 6 | 790 | 0 | 0 | 0 |
mdm_mock2 | 0 | 7 | 652 | 0 | 0 | 0 |
mdm_mock3 | 0 | 4 | 173 | 0 | 0 | 0 |
mdm_mock4 | 0 | 4 | 811 | 0 | 0 | 0 |
react6 | 1 | 79 | 3015 | 1 | 1 | 1 |
tcc <- t(cc)
pctcc <- apply(tcc,2,function(x) { x/sum(x)*100} )
pctcc %>% kbl(caption="cell proportions") %>% kable_paper("hover", full_width = F)
alv_active1 | alv_active2 | alv_active3 | alv_bystander1 | alv_bystander2 | alv_bystander3 | alv_bystander4 | alv_latent1 | alv_latent2 | alv_latent3 | alv_latent4 | alv_mock1 | alv_mock2 | alv_mock3 | mdm_active1 | mdm_active2 | mdm_active3 | mdm_active4 | mdm_bystander1 | mdm_bystander2 | mdm_bystander3 | mdm_latent1 | mdm_latent2 | mdm_latent3 | mdm_mock1 | mdm_mock2 | mdm_mock3 | mdm_mock4 | react6 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Basophils | 0.1058201 | 0.000000 | 0.1402525 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.054615 | 0.000000 | 0.000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0853242 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
Dendritic cells | 3.5978836 | 1.288245 | 3.9270687 | 2.7384325 | 1.7825312 | 2.3809524 | 3.4423408 | 1.595745 | 1.5495868 | 1.419989 | 2.773333 | 1.848049 | 2.363184 | 3.4711965 | 1.931818 | 0.6263048 | 3.087886 | 0.9174312 | 3.212851 | 2.323838 | 8.7037037 | 1.4505119 | 0.8792966 | 3.160919 | 0.7537688 | 1.062215 | 2.259887 | 0.4907975 | 2.5500323 |
Monocytes | 96.0846561 | 98.711755 | 95.9326788 | 97.0727101 | 98.0392157 | 97.4206349 | 96.0413081 | 98.404255 | 98.3987603 | 98.525396 | 97.226667 | 98.151951 | 97.636816 | 96.3072378 | 98.068182 | 99.3736952 | 96.912114 | 99.0825688 | 96.787149 | 97.676162 | 90.7407407 | 98.4641638 | 99.1207034 | 96.839080 | 99.2462312 | 98.937785 | 97.740113 | 99.5092025 | 97.3208522 |
Neutrophils | 0.2116402 | 0.000000 | 0.0000000 | 0.1888574 | 0.0000000 | 0.1984127 | 0.5163511 | 0.000000 | 0.0516529 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0738552 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.5555556 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
NK cells | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0891266 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.1477105 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
Progenitors | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0891266 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
focus <- meta_inf[grep("latent",rownames(meta_inf),invert=TRUE),]
focus <- focus[grep("bystander",rownames(focus),invert=TRUE),]
focus_mdm <- focus[grep("mdm",rownames(focus)),]
focus_alv <- focus[grep("alv",rownames(focus)),]
mono_focus_mdm <- mono[,which(colnames(mono) %in% rownames(focus_mdm))]
mono_focus_alv <- mono[,which(colnames(mono) %in% rownames(focus_alv))]
# mono_focus_mdm
mono_focus_mdm <- FindVariableFeatures(mono_focus_mdm, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
mono_focus_mdm <- RunPCA(mono_focus_mdm, features = VariableFeatures(object = mono_focus_mdm))
## Warning: The following 45 features requested have zero variance; running
## reduction without them: AC073359.1, KCNIP1, AC093916.1, PLK1, NAALADL2, ANK2,
## PIWIL2, MDN1, CCL5, AC013287.1, LINC01937, FANCM, MT-ND1, FGL2, RABGAP1L,
## TNFRSF4, AHNAK2, LINC02221, LRRK2, L3MBTL4, ZC3H12C, LINGO1, RANBP17, BCL2,
## CREB5, CTHRC1, AC023194.3, VIPR1, DIAPH2, TESC, SLC44A5, SFMBT2, RRAS2,
## ARHGAP26, AC024901.1, LINC00910, ARRDC4, RIIAD1, AL121820.1, AC048382.1,
## AC004448.3, SUCLG2-AS1, LINC01605, BDNF-AS, C22orf24
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: HIST1H2BJ, OASL, ENG, CENPA, HIST1H3J,
## NCALD, SSBP3, PCBP3, AL355388.1, LINC00885, AL358944.1, SCAI, IFT80, C5AR2,
## KLF12, AL009177.1, SLC2A14, LINC02739, LINC01320, TVP23A, CCND2, ZFPM2-AS1,
## AC108134.2, TTK, PPP6R2, TNRC6C, STXBP1, AC106793.1, HCAR2, RBM6, AC064805.2,
## PNPLA7, GPR155, ANKRD34B, KLRG1, SETD2, AL021917.1, AC015660.2, ZBTB41, EIF2B3,
## NUP93, AGPAT5, AC025262.1, LINC02853, ESR2, CLPX, ZZEF1, MARCH1, FAM118B,
## CWC22, RASIP1, ATP6V0A2, AL353759.1, SH3GL1, TMEM212-AS1, ASTL, PKIA-AS1,
## LINC02585, PCA3, P2RY12, RPS4Y1, MT-ND5, AL357873.1, OBI1-AS1, AC007364.1,
## SF3B3, CHRNB4, TMEM220-AS1, PARD3, GYG2, AL157911.1, AC034213.1, CCNB2, AURKC,
## RAP1GDS1, NR2F1-AS1, DUSP16, IFT140, AC072022.2, AC078923.1, HLA-C, BTBD2,
## MT-ND2, AP000462.1, FAM184A, ZNF276, AC133550.2, AC093766.1, MAP4K1,
## AC108879.1, SRRM2-AS1, RORA, CFAP52, LPL, DOK5, CLIC5, SIGLEC10, CREM,
## LINC01054, LRRC9, CLEC4C, GRM7, PSME2, STAT5B, AC117383.1, WWP1, TTC21B-AS1,
## NYX, AC099552.1, SNAP25-AS1, FHOD3, NLRP4, ZNF385D-AS1, ACTR3C, LDHA, PDIA4,
## AC100849.2, CU638689.5, RGS18, DAPK2, JAKMIP3, MRC2, CCDC7, AP000446.1,
## MAP3K15, DPEP3, SETBP1, AL592295.3, LINC02015, AC092490.1, PACS1, MASP2, CDK19,
## AC089983.1, ANK3, AC240274.1, NELL2, APCDD1L, P2RY13, AL360015.1, AC105001.1,
## MMP28, FDXACB1, RAB10, SRGN, AC104459.1, ARC, CPLX1, AC007091.1, SAMD11,
## GAS1RR, SF3B1, AC013799.1, AL031658.1, CRISPLD2, AC114977.1, AC007389.1, OR3A2,
## AL589740.1, AHR, IFI6, GFRA2, MIR155HG, RASSF4, AC117473.1, CARD11, AL132775.1,
## TUBB3, KIF21A, TPTE2, AC016074.2, GTSF1, AC019117.1, AL391807.1, TEX15,
## AC005393.1, ADRA2B, MANF, UBASH3B, PLEK, AC005264.1, CEP112, ZDHHC17,
## AC011509.2, CACNA2D3, AC084809.2, CD93, AC007128.2, AC073569.1, NR2F2-AS1,
## PRRT2, AC026333.3, LINC01600, TSPAN8, Z99758.1, GNG4, CFAP69, AC011365.1,
## BCL2A1, P3H2, CEP135, TFG, JAKMIP2-AS1, AL049828.1, AC137770.1, DHX8, U62317.4,
## MT-CO3, IGSF6, DHCR24, TMEM45A, COX5B, SULF2, TRIM37, CPE, AC096570.1, SRGAP3,
## SESN3, LINC01098, TACC3, MPDZ, GAL3ST4, AP000829.1, BMPR1B, LINC00237,
## LINC00649, STPG2, MKX, CORIN, SNHG12, RTKN2, PRTG, BCAR3, OSBPL6, HIST1H2BG,
## LINC01917, LINC01855, MLLT3, LNX1, PAX8-AS1, SIGLEC7, MAPK8, HIST1H2AC,
## AL136418.1, EPHA6, MCF2L-AS1, DPF1, AC006427.2, AC079062.1, DARS, FKBP1B,
## LINC01353, Z94721.1, ABCA13, PODN, MT-ATP6, SOX15, ACMSD, KCNMB4, TASP1,
## AC021231.1, VGLL3, HIST2H2BE, CHRM3, AC087683.3, CREG2, PPFIA2, LINC01891,
## TRIM71, AC124017.1, AC008115.2, AC009315.1, SCAPER, LINC02112, CNDP1,
## AC073475.1, AKR7A2, RAP2C-AS1, HCAR3, TMCC1, AC021086.1, RUFY2, RFX7, TMEM71,
## MCCC2, TUBB4B, LINC01191, MNAT1, MRVI1, SOX21-AS1, CKMT1A, ATF3, KDM6A, CD72,
## AL031710.2, AL031005.2, KIAA0825, MAP1A, AC010834.2, DNAH3, ZBTB46, AC090630.1,
## PPP1R12B, AC129803.1, GSN, PDE11A, AL096794.1, PTPRB, ZHX2, PCP4L1, LINC00958,
## MNDA, PTPRG-AS1, KCNJ1, GAK, GLIPR1, RALGAPA2, COL8A2, LYPLAL1-DT, TNS3, MREG,
## RHOF, ITGA4, C11orf45, AL356379.2, AC011476.3, TUBB4A, LGALSL, SCFD2, ALDH1A1,
## SPEG, BFSP2, AC019330.1, AC063923.2, ECE2, AF131216.1, AC103726.2, PROCR,
## MARCH3, AEBP2, AC246817.1, SPIB, POF1B, ANGPTL4, COL28A1, AL157834.1, ACTB,
## PRKG2, FILIP1L, RFX3, GMDS, SOCS3, SPTLC3, PAN3, OXSR1, LINC00571, AC097654.1,
## AC092811.1, UBE2T, SPAG7, MCM9, MAFB, LINC00407, LILRB2, NEGR1, PLA2G5, MAPK13,
## KIAA1841, TRIO, AC079584.1, LINC02851, NETO2, CYP4F22, CLDN4, TMEM131L, RNF180,
## SLC26A7, AL365295.2, AC098588.3, C3orf79, ADORA2B, SPOCK3, RNF24, GRK3, NPRL3,
## TMEM72-AS1, SGO1, ENPEP, ITM2A, BUB1, AC114763.1, PCLO, XKR9, AL136298.1,
## RAPGEF1, ST3GAL6, GLRX, LRRIQ4, AP002793.1, LINC00987, PIP4K2C, CDT1,
## AL137076.1, AL158839.1, AL589765.6, AC096639.1, AC114810.1, AL133243.1,
## AC009229.3, AC012358.1, AC063944.3, AC092902.6, AC128709.2, AC021220.2,
## ADAMTS3, AC093801.1, AC116351.1, LINC02149, AC113386.1, AC011374.3, LINC02571,
## Z84484.1, MRAP2, AL357992.1, AL078582.1, AC002480.2, ITPRID1, DEFB136,
## AC084026.1, CALB1, AF178030.1, AF235103.1, AL353764.1, TMEM246-AS1, PPP3R2,
## AL353803.1, AL731559.1, AL121748.1, MTRNR2L7, LINC02625, AC016816.1, SYCE1,
## OR10A4, LINC02751, SAA4, AC013714.1, OR8J1, AC024940.1, LINC01479, AC012464.3,
## AC073863.1, C1QTNF9-AS1, SMAD9-IT1, LINC00563, AL161717.1, CLYBL-AS2,
## AL442125.2, KLHL33, CMA1, AL163953.1, LINC02280, AC091167.5, BAIAP3,
## AL133297.1, AC012186.2, AC092378.1, AC129507.2, RAI1-AS1, AC007952.6,
## AC243773.2, MIR4527HG, AC090409.1, AC105094.2, OR7E24, ANGPTL8, AL021396.1,
## NUDT11, LINC00266-4P, CFD, TFRC, TRAF2, PTPN13, CALR, NKAIN1, STAG1, ATP6V0D2,
## ATP13A3, EXD2, ZC4H2, E2F1, HERPUD1, KIF20B, CHM, ADIRF, BACH1, COL27A1,
## ADAMTS10, CDKAL1, LINC01762, AL645929.2, MMD2, NIPAL2, NFKB1, ENTPD1, LY96,
## ZFP36L1, SH3D19, ITGA9, MGAT5, ARRDC3-AS1, CD200R1, AC092134.3, LGALS3, PKP2,
## PFKFB3, SEMA6D, COBL, MYADM, ELANE, TNIK, TIMP4, FAM135A, AHCYL2, FAM110B,
## PRDM1, LINC02398, LINC02798, STUM, MDH1, KL, DUSP5, LINC01572, MTUS1, GADD45G,
## CACNA2D2, RALYL, PRH1, CASC2, GNG7, CASP1, SPATA5, DOCK10, PLBD1, LINC02476,
## AC083836.1, PCNX2, MCF2, LINC00511, SIRPD, MOCOS, AL158068.2, AC117465.1,
## AC007391.2, AC073050.1, DRD5, RPS6KA6, FAM189A2, SLC41A2, CFAP43, ST3GAL3,
## MT-CYB, FANCC, AC007405.3, EDA, PMAIP1, PRORP, ERBB4, TRAF1, UFL1-AS1, DYM,
## MFSD4B, CARMIL3, AC097515.1, NCAPD3, PAFAH1B1, AL445584.2, CALM3, NPTX2,
## LAMC1-AS1, AC099782.2, WWTR1-AS1, AC104785.1, AL357054.5, AC010719.1,
## AC073210.3, AP005436.3, AP002884.1, LUM, AC063943.1, PLUT, AC011939.3, NPIPB8,
## AC092127.1, AC104423.1, L1CAM, SCLT1, MATN1, AC020743.2, AC021546.1, PRSS3,
## LUZP4, LINC02457, MCTP2, CRYM, PKN2, CAMK2D, SLC7A8, IGF2R, AFF1, ITPR2,
## AC024084.1, AC096992.3, STMN2.
## PC_ 1
## Positive: SLC35F1, PRDX6, TUBA1B, CRABP2, TUBA1A, GAL, MT-ATP8, C15orf48, FABP3, DUSP2
## LILRA1, PTGDS, MYL9, HPGDS, MT2A, UCHL1, CRIP1, MMP7, CCNA1, MT1E
## NCAPH, S100A4, LINC02244, CA2, SEL1L2, AC105402.3, CCL7, RGCC, RNF128, TMEM114
## Negative: FMN1, NEAT1, FHIT, RAD51B, ARL15, MALAT1, SNTB1, MBD5, AL035446.2, FTX
## EXOC4, FNDC3B, VTI1A, PDE4D, FRMD4B, GMDS-DT, COP1, JMJD1C, TBC1D8, SLC8A1
## SLC22A15, DENND1A, DENND4C, DOCK3, CBLB, DOCK4, VPS13B, LRMDA, SBF2, COG5
## PC_ 2
## Positive: GDF15, PSAT1, B4GALT5, TM4SF19, PHGDH, FDX1, SLAMF9, SNCA, TCEA1, S100A10
## TXN, BEX2, OCIAD2, PHLDA1, RAB38, UGCG, TPM4, SUPV3L1, DRAXIN, RETREG1
## NMRK2, HES2, NMB, CCPG1, CTSL, SCD, CLGN, FABP4, gene-HIV1-2, MARCO
## Negative: TGFBI, HLA-DRA, CEBPD, HLA-DPB1, CD74, HLA-DPA1, HLA-DRB1, SIPA1L1, MS4A6A, LYZ
## EPB41L3, SSBP2, ST8SIA4, AOAH, PTGDS, HLA-DQB1, KCNMA1, RGS2, HLA-DRB5, SAMSN1
## RNASE1, FOS, PDE4B, C1QA, TCF4, MS4A7, VMO1, RCBTB2, FPR3, RNASE6
## PC_ 3
## Positive: BTG1, G0S2, MARCKS, CD14, ARL4C, CLEC4E, TGFBI, CXCR4, SDS, FUCA1
## SEMA4A, TIMP1, HIF1A, VMO1, STMN1, NUPR1, MEF2C, MS4A6A, RNASE1, P2RY8
## TCF4, OLFML2B, CLEC4A, PDK4, NAMPT, HES4, BLVRB, LGMN, MPEG1, FPR3
## Negative: NCAPH, CRABP2, RGCC, LINC01010, LINC02244, CHI3L1, GCLC, DUSP2, TMEM114, NUMB
## RGS20, LRCH1, PTGDS, AC092353.2, S100A4, MICAL3, ASAP1, CYP1B1, AC005280.2, DOCK3
## SERTAD2, CRIP1, RASAL2, TIMP3, TMEM117, NOS1AP, TCTEX1D2, SLC28A3, MYO1E, GPC4
## PC_ 4
## Positive: AC078850.1, TUBA1B, ACTG1, SLC35F1, INSIG1, HSP90B1, TUBB, TMEM176A, CCL7, FABP4
## TPM4, MGLL, TUBA1A, LGMN, CD36, C3, CD300LB, TIMP3, FBP1, IL18BP
## TUBA1C, PHLDA1, CSF1, CADM1, LINC01091, CTSL, TMEM176B, HLA-DRB1, ALCAM, IL4I1
## Negative: MT-ATP8, HES2, BEX2, PDE4D, CLGN, S100P, PSAT1, CCPG1, NIBAN1, C5orf17
## NUPR1, CPD, PHGDH, XIST, ST6GALNAC3, RETREG1, EPHB1, HES4, MT1G, LINC02154
## SUPV3L1, PECR, CLVS1, CEP70, ALKAL2, TDRD3, ME1, AL035446.2, DMGDH, G0S2
## PC_ 5
## Positive: TYMS, PCLAF, MKI67, CENPF, CEP55, BIRC5, PRC1, NUSAP1, KIF4A, CDKN3
## DLGAP5, CENPM, CDK1, TK1, HMMR, TPX2, CENPK, ASPM, PTTG1, ANLN
## RAD51AP1, CCNA2, MAD2L1, RRM2, CIT, BUB1B, DIAPH3, SHCBP1, MYBL2, NCAPG
## Negative: gene-HIV1-2, gene-HIV1-1, CCL7, AC078850.1, PHLDA1, LINC01091, TNFRSF9, TCTEX1D2, C3, CSF1
## SLC35F1, MADD, CD36, IL1RN, MACC1, SPOCD1, AC005062.1, DPYSL3, INSIG1, LIMA1
## TIMP3, TPM2, TPM4, CCL3, B4GALT5, CD300LB, TM4SF19, PCSK6, SDS, CADM1
DimHeatmap(mono_focus_mdm, dims = 1:2, cells = 500, balanced = TRUE)
DimHeatmap(mono_focus_mdm, dims = 3:4, cells = 500, balanced = TRUE)
ElbowPlot(mono_focus_mdm)
mono_focus_mdm <- FindNeighbors(mono_focus_mdm, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
mono_focus_mdm <- FindClusters(mono_focus_mdm, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono_focus_mdm <- RunUMAP(mono_focus_mdm, dims = 1:4)
## 13:58:28 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:58:28 Read 4605 rows and found 4 numeric columns
## 13:58:28 Using Annoy for neighbor search, n_neighbors = 30
## 13:58:28 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:58:28 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea732d5a26
## 13:58:28 Searching Annoy index using 1 thread, search_k = 3000
## 13:58:30 Annoy recall = 100%
## 13:58:30 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:58:32 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:58:32 Commencing optimization for 500 epochs, with 162990 positive edges
## 13:58:35 Optimization finished
DimPlot(mono_focus_mdm, reduction = "umap", label=TRUE)
DimPlot(mono_focus_mdm, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE)
DimPlot(mono_focus_mdm, group.by="sample" , reduction = "umap", label=TRUE)
# mono_focus_alv
mono_focus_alv <- FindVariableFeatures(mono_focus_alv, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
mono_focus_alv <- RunPCA(mono_focus_alv, features = VariableFeatures(object = mono_focus_alv))
## Warning: The following 40 features requested have zero variance; running
## reduction without them: F13A1, TMEM163, FCGBP, MARCKS, TIMP3, LINC01954, S100P,
## TACSTD2, FAM9B, SLC17A6, LINC00636, MAML2, GAPDH, MND1, SMS, CD226, AC010834.2,
## CCDC57, CD5L, SEMA3D, MFSD11, LINC01320, ELMO1, SLC6A7, AC092813.1, CLSTN2,
## CDH12, AC007381.1, EYA2, SLC23A2, LINC02777, ANKRD7, LINC01198, TSGA13, IGSF6,
## LINC01948, HMSD, FOLR2, CACNB4, GRB10
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: STEAP1B, NUP210L, CPNE8, AC083837.1,
## AP002075.1, HCAR2, CD28, AL096794.1, PPP1R1C, SFMBT1, IQCA1, ACVR2A, SLC6A16,
## SPINK6, TENM4, IFI6, EFCAB7, NPTX2, C11orf49, IL6, DNAH12, COL25A1,
## SLC7A14-AS1, LINC02068, DNAJC9, HIST1H2BG, AC076968.2, GNLY, LINC02269,
## LINC01605, EMP3, MS4A5, DSG2, INO80, ACTN2, AC207130.1, BANK1, BTNL8,
## LINC00350, AC092821.3, AC013391.2, KIF6, ARMC8, SMCHD1, TEK, ADAM22, EXO1,
## AC135050.3, MIF, COL27A1, GCH1, PLAAT3, CCDC85A, MSRB3, ERAP2, CNTNAP2, TRIM2,
## SCFD1, AC084740.1, AHCTF1, LINC02073, MIR4300HG, CTSW, MARCKSL1, CH25H, BRMS1L,
## MNDA, PMAIP1, AL353596.1, RAD51C, TUBB4A, GIGYF2, LINC00299, GPRC5D, SVEP1,
## LINC02789, KIAA1755, NEURL3, RFX3-AS1, OR10G3, USP10, MRPS6, ZC3H8, PTPN2,
## LINC01344, ASAP3, SHANK2, DYNC1H1, RFTN2, XDH, YJEFN3, SPOCK1, AC007785.1,
## MECP2, CR1, LPCAT1, TNFAIP3, AC019131.1, LINC01862, AC104596.1, AC079313.2,
## RPH3AL, NEK4, SH3TC2, URB1, LINC00511, STK32B, CFAP74, AC110992.1, HSF2BP,
## SPACA3, XPO5, RLF, PPEF1, AC084809.2, EML2-AS1, FBXO4, AC009292.2, ZPLD1,
## MBOAT4, NIPAL2, POU6F2, AC005264.1, PTX3, LGALS1, AC125421.1, COL4A5,
## SMAD1-AS2, DZIP1, HNRNPM, LINC02466, RGS8, LINC00842, NDRG1, INKA2-AS1, WDR90,
## TEPP, SOX6, ABCB1, TIMELESS, SYT12, LINC01999, SRL, ZSCAN32, PLEK, STS, NRG1,
## FCRLB, DCSTAMP, SLC44A5, ADCY3, PMPCA, CNR2, ZNF157, AL135926.2, STXBP6, CCNC,
## FO393415.3, IGSF23, DDN-AS1, NCMAP, LMNB1-DT, LINC01414, CSTB, PAX8-AS1, RFX8,
## CFD, GNRHR, FSD1, CLDN18, MID2, AC073569.2, GTF2IRD2, ABRAXAS2, WIPF3, HRH2,
## MIR2052HG, AKAP6, GATA2, AC011346.1, AC107081.2, AC041005.1, NNMT, HLTF-AS1,
## C2orf72, PLXDC1, PTPRB, HCAR3, TMEM37, PHACTR1, CTSZ, IL17RB, IGF2R, ANOS1,
## DACH1, AL662884.3, AC006441.4, RELN, C22orf34, ERI2, NCR3, PLCH1, EAF2, CKAP5,
## AOC1, UNC5D, COBL, PDE1C, AL592295.3, AC068228.3, SLC2A5, SVOP, AC022035.1,
## PSME2, GOLGA7B, AC004949.1, ANGPT1, AC009435.1, RGS7, AC008415.1, TNR,
## AC097518.2, CFAP57, CNTNAP5, SHCBP1L, TNC, AC007100.1, AC005753.2, NBAT1,
## AC109454.3, NIPAL4, AC025430.1, TPSAB1, SPRY4-AS1, ENTPD1, RAB7B, STPG2, MTFMT,
## NDRG2, ABLIM1, ADM, CDK6, GRID2, KIAA0825, UST, HSPA1B, ULBP2, LDHA,
## AL354733.2, MRC2, DHRS9, AC008609.1, PRDM14, AL161646.1, KRTAP10-4, LINC00378,
## Z99758.1, SHROOM4, AC010997.3, CTXND1, AL513166.1, BPIFC, MEI4, HIP1, GALNT14,
## RHOF, ZNF385D-AS1, FUT2, DENND2A, NEIL3, LGALS3, KLB, OASL, NSG2, CYP4F22,
## LINC01358, VIPR1, GPRIN3, GAPLINC, KIF21A, NRIP3, MIR3142HG, AC116563.1, KDR,
## AL049828.1, SLC4A1, AC110491.2, NRCAM, CRIM1, FOXM1, ZSCAN5A, ARHGAP6, GINS2,
## SCIN, LINC00973, AC073091.4, SH3TC2-DT, LINGO2, CHODL, SERINC2, GRM7, CHDH,
## ITGA2, PARD3, TROAP, APOM, RXFP1, GLIPR1, POU2F2, LINC01643, MEP1A, NES,
## PCDH15, MSH4, AL033530.1, LINC01900, AC092957.1, AC015908.2, SFTPD-AS1,
## AC006460.1, LINC00571, PRKCG, FAXC, ATP1B2, AC093307.1, IAPP, IKZF3, GPX3, SLA,
## YLPM1, CLEC7A, CLPB, RBL1, BRCA1, CRIPT, AC108066.1, CDHR3, AL109930.1, ASMT,
## AIM2, AC011893.1, ABCG2, AC117453.1, CHAC1, SLC22A2, NCAM2, HMOX1, TMEM45A,
## AC104041.1, AC093898.1, DPYS, PCLO, ELOC, LIPG, FAF1, SAMD12, MCAM, LIN28B-AS1,
## CHD5, LINC00519, AF274853.1, AC016831.7, MREG, MOSPD1, SIK2, CD1E, AL359220.1,
## ID2, AC008443.9, AMPD1, PRRX2, AC124852.1, SLC23A1, GRIK5, AC099489.1, ZNF365,
## TNFRSF11A, AC002454.1, AC093010.2, UPK1A-AS1, AC092718.1, LCP1, DHCR24, RHOD,
## PARN, AC007529.2, LINC01924, AC007402.1, ARL9, TGFB3, MAS1, SEMA6D, UBE2E2,
## FHAD1, TBC1D24, AC084816.1, GRHL2, AL360015.1, TDRD9, AC024901.1, SRGAP3,
## MEIS1, XKR9, AC079742.1, LINC02698, ING3, GNAI1, AC246817.1, RNF212,
## AL713998.1, ANKH, SLC35F3, MX2, LINC02109, AC130456.2, AHCYL2, AL136456.1,
## KIAA1217, SPSB1, SLC4A8, AC006333.1, ZAR1L, RBPJL, AC026391.1, ERBB4,
## LINC02752, PLTP, SNAI3, OPRD1, LINC02621, TMEM233, AC018618.1, FBLN1,
## AL355981.1, TMEM213, AC087897.2, DEGS2, WDFY4, AC099560.1, TNIK, SLC9A2,
## LINC01572, AL591845.1, GLT1D1, RASL10A, RALGAPA2, GALNTL6, COL8A2, GLUL, MAP1A,
## SEPTIN4-AS1, ST3GAL6, LINC02742, RNF144A, PTP4A2, LGR4, LINC01800, KDM7A,
## CNGA4, BX004807.1, DNAH2, ELOVL5, AC010343.3, LINC02805, GNGT1, FRRS1, MB21D2,
## CAMK1, CPEB2-DT, CEP126, CSMD2, AP001496.1, FBXW7, EXOSC10, ELF5, NANOS1,
## RP1L1, CFI, PROSER2-AS1, AC009107.1, TEKT1, AC113137.1, SOX5, STXBP5L, SPOCK3,
## TMEM132C, INPP5J, CLDN4, MARCH1, IL3RA, TNFRSF12A, S100A6, FBXO43, PPP1R14C,
## AP000424.1, EOGT, AC079163.2, AC087280.2, GLCCI1, PDLIM4, SKA3, STUM, C1orf143,
## CD96, ASTN2, GM2A, GRAMD2A, LINC01933, ACTB, LINC00894, DNAJC1, LIMCH1, HECW1,
## NR6A1, AL645465.1, HIST2H2BE, AL805961.1, AC239860.4, SCG2, AC021134.1,
## AC137810.1, AC145146.1, AL591519.1, AC108860.2, AL157702.2, AC068305.2, OTOGL,
## AC090515.6, LHCGR, ZBED9, AL162411.1, CLSTN3, AMPD3, HIST1H2AC, RASSF4,
## TMEM131L, C11orf45, SH3PXD2B, GALNT18, CDT1, INPP4B, SPOPL, FA2H, SLC30A8,
## FAM92B, AC016587.1, BICD1, ZNF431, FCMR, LPAR1, AL353595.1, CIDEC, STMND1,
## SSPO, NLRP7, GPAT3, WDR54, MYO16-AS1, AC005808.1, NUDT10, AC096531.2,
## LINC01276, AL354811.2, RBM11, PDE3A, TRMT5, JAML, AC019197.1, CORO1A, SERPINA1,
## KCNJ1, EMP1, LMO4, CDC5L, AC055855.2, BTG3-AS1, RBPJ, SYNE1, KCP, GRXCR2, EML4,
## AP001636.3, SCFD2, LINC01169, GABRR3, AC090099.1, ROBO4, AC092131.1, TPPP3,
## ATP13A3, TYW1B, HIST1H1C, AC005280.1, ITSN2, SLC25A23, FABP6, MPDZ, AC011476.3,
## EFNA2, AEBP2, AP000812.1, AC009315.1, HTRA4, ATP1B1, AC021851.2, AC092746.1,
## AHRR, BTBD11, MARCH3, AL122003.2, PKIB, MIR155HG, FRMD6-AS1, AC104809.2, BARD1,
## ADD2, RNF217-AS1, MS4A4A, AC079793.1, AL731563.2, AC092801.1, AC068587.4,
## AC007262.2, MED13L, LAMA3, NME5, LMCD1-AS1, LINC02224, AXL, ZNF543, ZNF385B,
## SH3BP5, SNX10, LINC02733, PLAT, HPN, RFX7, AL592431.2, ERLEC1, CHD9, SRPX2,
## KCNF1, ABCA6, AC097487.1, KCNN2, AGMO, DMRT1.
## PC_ 1
## Positive: NUPR1, FABP4, STMN1, PSAT1, PHGDH, CLGN, G0S2, IFITM3, BTG1, ALDH2
## DDIT4, CAMP, BLVRB, ADAMDEC1, TRIB3, GDF15, BEX2, AC073359.1, BEX3, CD14
## DNAJC5B, CARD16, CLEC4A, RAB6B, TIMP1, FAH, GYPC, FXYD2, RBP1, CES1
## Negative: RASAL2, DOCK3, TMEM117, AC092353.2, DPYD, ASAP1, LINC01374, PLXDC2, LRMDA, CPEB2
## ATG7, FMNL2, LRCH1, NEAT1, MITF, ARL15, TPRG1, ATXN1, EXOC4, NUMB
## DENND1B, MALAT1, PPARG, VWA8, FHIT, MYO1E, RABGAP1L, MEF2A, JARID2, PDE3B
## PC_ 2
## Positive: gene-HIV1-1, CCL22, IL6R-AS1, GAL, gene-HIV1-2, AL157912.1, SLC20A1, DUSP13, IL1RN, TRIM54
## CYTOR, MYL9, CSF1, POLE4, MIR4435-2HG, CIR1, TM4SF19, BIN2, NMRK2, GPC3
## ACAT2, ATP6V1D, C3, SCD, NCAPH, TCTEX1D2, LAYN, UPP1, C15orf48, TNFSF14
## Negative: SLC8A1, MRC1, AL356124.1, TRPS1, SAMSN1, RCBTB2, PDGFC, FCGR3A, FCHSD2, SPRED1
## ME1, ATP8B4, CTSC, NRP1, CAMK1D, XYLT1, LYZ, AIF1, ARHGAP15, CCDC102B
## DOCK4, SLCO2B1, TGFBI, HLA-DPA1, C1QA, FOS, TCF7L2, TMEM236, HLA-DRA, HLA-DPB1
## PC_ 3
## Positive: CTSL, MARCO, LGMN, BCAT1, TPM4, SAT1, HAMP, BLVRB, SLC11A1, CTSB
## TXN, FMN1, CCL3, CTSC, STMN1, MSR1, NUPR1, FABP4, S100A9, SNCA
## FDX1, PLAU, FCGR3A, CD164, UGCG, FPR3, CD36, FAH, GYPC, SLAMF9
## Negative: AC067751.1, CRIP1, CLU, MMP7, DUSP2, KCNMA1, PTGDS, LINC02408, TRIM54, NCAPH
## SERTAD2, C2orf92, GCLC, RNF128, KCNQ5, ZNF366, C1QTNF4, ST5, KCNA2, NFATC3
## IL1RN, AC092944.1, RCBTB2, DPH6-DT, TMEM176B, MSI2, LINC00910, NOS1AP, FGL2, RGS20
## PC_ 4
## Positive: XIST, LINC02320, CCDC26, AC008591.1, OSBP2, MIR646HG, LINC01340, AL136317.2, LINC01708, SKAP1
## ZFPM2, KCNMB2-AS1, LIX1-AS1, AC068413.1, SEL1L2, PKD1L1, LINC01374, AC012668.3, LINC00607, CLMP
## AP005262.2, AL365295.1, UBA6-AS1, LINC01135, LINC00923, SASH1, STOX2, PLXNA2, AC098829.1, AC079142.1
## Negative: TUBA1A, IL1RN, TUBA1B, ALOX5AP, ACTG1, MMP7, SLC35F1, MMP9, TMEM176B, TUBB
## TMEM176A, HLA-DRB1, RGCC, CD74, CTSB, GCLC, HLA-DRA, CRIP1, HLA-DPA1, PTGDS
## LINC00910, RNF128, AC007952.4, TEX14, AIF1, S100A4, DUSP2, ALCAM, H2AFZ, STAC
## PC_ 5
## Positive: gene-HIV1-1, gene-HIV1-2, CTSB, AIF1, IL6R-AS1, SPRED1, MRC1, C1QA, IL1RN, STAC
## AL157912.1, LINC02345, SLCO2B1, FPR3, HLA-DRA, CCL2, LGMN, ISG15, ALCAM, C1QB
## C15orf48, CD74, RTN1, HLA-DRB1, IL4I1, HLA-DPA1, CCL7, MMP19, CTSC, TNFSF13B
## Negative: CLSPN, TYMS, TK1, MKI67, PCLAF, DIAPH3, SHCBP1, RRM2, CENPK, FAM111B
## HELLS, ATAD2, CDK1, MYBL2, NCAPG, CDCA2, ESCO2, CCNA2, TOP2A, BIRC5
## KIF11, CENPF, POLQ, CENPM, HMMR, WDR76, CIT, DTL, KNL1, NUSAP1
DimHeatmap(mono_focus_alv, dims = 1:2, cells = 500, balanced = TRUE)
DimHeatmap(mono_focus_alv, dims = 3:4, cells = 500, balanced = TRUE)
ElbowPlot(mono_focus_alv)
mono_focus_alv <- FindNeighbors(mono_focus_alv, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
mono_focus_alv <- FindClusters(mono_focus_alv, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono_focus_alv <- RunUMAP(mono_focus_alv, dims = 1:4)
## 13:58:44 UMAP embedding parameters a = 0.9922 b = 1.112
## 13:58:44 Read 5250 rows and found 4 numeric columns
## 13:58:44 Using Annoy for neighbor search, n_neighbors = 30
## 13:58:44 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:58:45 Writing NN index file to temp file /tmp/RtmpuopyiY/file2952ea3ead0f57
## 13:58:45 Searching Annoy index using 1 thread, search_k = 3000
## 13:58:46 Annoy recall = 100%
## 13:58:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 13:58:49 Initializing from normalized Laplacian + noise (using RSpectra)
## 13:58:49 Commencing optimization for 500 epochs, with 185554 positive edges
## 13:58:52 Optimization finished
DimPlot(mono_focus_alv, reduction = "umap", label=TRUE)
DimPlot(mono_focus_alv, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE)
DimPlot(mono_focus_alv, group.by="sample" , reduction = "umap", label=TRUE)
We are going to use muscat for pseudobulk analysis. First need to convert seurat obj to singlecellexperiment object. Then summarise counts to pseudobulk level.
sce <- Seurat::as.SingleCellExperiment(comb, assay = "RNA")
head(colData(sce),2)
## DataFrame with 2 rows and 13 columns
## orig.ident nCount_RNA nFeature_RNA
## <factor> <numeric> <integer>
## mdm_mock1|AAACGAATCACATACG mac 72710 7097
## mdm_mock1|AAACGCTCATCAGCGC mac 49092 6276
## cell sample RNA_snn_res.0.5
## <character> <character> <factor>
## mdm_mock1|AAACGAATCACATACG mdm_mock1|AAACGAATCA.. mdm_mock1 3
## mdm_mock1|AAACGCTCATCAGCGC mdm_mock1|AAACGCTCAT.. mdm_mock1 11
## seurat_clusters cell_barcode
## <factor> <character>
## mdm_mock1|AAACGAATCACATACG 3 mdm_mock1|AAACGAATCA..
## mdm_mock1|AAACGCTCATCAGCGC 11 mdm_mock1|AAACGCTCAT..
## monaco_broad_annotation monaco_broad_pruned_labels
## <character> <character>
## mdm_mock1|AAACGAATCACATACG Monocytes Monocytes
## mdm_mock1|AAACGCTCATCAGCGC Monocytes Monocytes
## monaco_fine_annotation monaco_fine_pruned_labels
## <character> <character>
## mdm_mock1|AAACGAATCACATACG Classical monocytes Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC Classical monocytes Classical monocytes
## ident
## <factor>
## mdm_mock1|AAACGAATCACATACG 3
## mdm_mock1|AAACGCTCATCAGCGC 11
colnames(colData(sce))
## [1] "orig.ident" "nCount_RNA"
## [3] "nFeature_RNA" "cell"
## [5] "sample" "RNA_snn_res.0.5"
## [7] "seurat_clusters" "cell_barcode"
## [9] "monaco_broad_annotation" "monaco_broad_pruned_labels"
## [11] "monaco_fine_annotation" "monaco_fine_pruned_labels"
## [13] "ident"
#muscat library
pb <- aggregateData(sce,
assay = "counts", fun = "sum",
by = c("monaco_broad_annotation", "sample"))
# one sheet per subpopulation
assayNames(pb)
## [1] "Basophils" "Dendritic cells" "Monocytes" "Neutrophils"
## [5] "NK cells" "Progenitors"
t(head(assay(pb)))
## gene-HIV1-1 gene-HIV1-2 MIR1302-2HG FAM138A OR4F5 AL627309.1
## alv_active1 2 29 0 0 0 0
## alv_active2 0 0 0 0 0 0
## alv_active3 0 2 0 0 0 0
## alv_bystander1 0 0 0 0 0 0
## alv_bystander2 0 0 0 0 0 0
## alv_bystander3 0 0 0 0 0 0
## alv_bystander4 0 0 0 0 0 0
## alv_latent1 0 0 0 0 0 0
## alv_latent2 0 0 0 0 0 0
## alv_latent3 0 15 0 0 0 0
## alv_latent4 0 0 0 0 0 0
## alv_mock1 0 0 0 0 0 0
## alv_mock2 0 0 0 0 0 0
## alv_mock3 0 0 0 0 0 0
## mdm_active1 0 0 0 0 0 0
## mdm_active2 0 0 0 0 0 0
## mdm_active3 0 0 0 0 0 0
## mdm_active4 0 0 0 0 0 0
## mdm_bystander1 0 0 0 0 0 0
## mdm_bystander2 0 0 0 0 0 0
## mdm_bystander3 0 0 0 0 0 0
## mdm_latent1 0 2 0 0 0 0
## mdm_latent2 0 0 0 0 0 0
## mdm_latent3 0 0 0 0 0 0
## mdm_mock1 0 0 0 0 0 0
## mdm_mock2 0 0 0 0 0 0
## mdm_mock3 0 0 0 0 0 0
## mdm_mock4 0 0 0 0 0 0
## react6 0 0 0 0 0 0
plotMDS(assays(pb)[[3]], main="Monocytes" )
par(mfrow=c(2,3))
dump <- lapply(1:length(assays(pb)) , function(i) {
cellname=names(assays(pb))[i]
plotMDS(assays(pb)[[i]],cex=1,pch=19,main=paste(cellname),labels=colnames(assays(pb)[[1]]))
})
par(mfrow=c(1,1))
pbmono <- assays(pb)[["Monocytes"]]
head(pbmono)
## alv_active1 alv_active2 alv_active3 alv_bystander1 alv_bystander2
## gene-HIV1-1 29173 20670 14188 0 0
## gene-HIV1-2 120841 105818 103519 0 0
## MIR1302-2HG 0 0 0 0 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 16 11 27 20 38
## alv_bystander3 alv_bystander4 alv_latent1 alv_latent2 alv_latent3
## gene-HIV1-1 0 0 380 1777 1654
## gene-HIV1-2 0 0 2126 12238 13801
## MIR1302-2HG 0 0 0 1 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 9 14 6 65 59
## alv_latent4 alv_mock1 alv_mock2 alv_mock3 mdm_active1 mdm_active2
## gene-HIV1-1 2215 92 115 1038 11089 7725
## gene-HIV1-2 26684 727 1653 6685 65635 42374
## MIR1302-2HG 0 0 0 0 0 0
## FAM138A 0 0 0 0 0 0
## OR4F5 0 0 0 0 0 0
## AL627309.1 76 46 26 52 54 37
## mdm_active3 mdm_active4 mdm_bystander1 mdm_bystander2
## gene-HIV1-1 4907 11852 0 0
## gene-HIV1-2 73030 35142 0 0
## MIR1302-2HG 0 0 0 0
## FAM138A 0 0 0 0
## OR4F5 0 0 0 0
## AL627309.1 30 5 65 51
## mdm_bystander3 mdm_latent1 mdm_latent2 mdm_latent3 mdm_mock1
## gene-HIV1-1 0 649 799 156 37
## gene-HIV1-2 6 6540 7091 3290 542
## MIR1302-2HG 0 0 0 0 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 19 72 78 28 61
## mdm_mock2 mdm_mock3 mdm_mock4 react6
## gene-HIV1-1 67 3 141 895
## gene-HIV1-2 517 100 1078 7528
## MIR1302-2HG 0 0 0 0
## FAM138A 0 0 0 0
## OR4F5 0 0 0 0
## AL627309.1 30 6 15 37
dim(pbmono)
## [1] 36603 29
hiv <- pbmono[1:2,]
pbmono <- pbmono[3:nrow(pbmono),]
Gene ontology.
#go <- gmt_import("c5.go.v2023.2.Hs.symbols.gmt")
#names(go) <- gsub("_"," ",names(go))
#wget https://ziemann-lab.net/public/tmp/go_2024-11.gmt
go <- gmt_import("go_2024-11.gmt")
MDM group.
pbmdm <- pbmono[,grep("mdm",colnames(pbmono))]
pb1m <- pbmdm[,c(grep("active",colnames(pbmdm)),grep("latent",colnames(pbmdm)))]
head(pb1m)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_latent1
## MIR1302-2HG 0 0 0 0 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 54 37 30 5 72
## AL627309.3 0 2 0 0 0
## AL627309.2 0 0 0 0 0
## mdm_latent2 mdm_latent3
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 78 28
## AL627309.3 2 2
## AL627309.2 0 0
pb1mf <- pb1m[which(rowMeans(pb1m)>=10),]
head(pb1mf)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_latent1
## AL627309.1 54 37 30 5 72
## AL627309.5 15 14 15 6 35
## LINC01409 129 114 104 27 153
## LINC01128 262 164 167 94 215
## LINC00115 22 6 7 4 23
## FAM41C 49 39 63 68 62
## mdm_latent2 mdm_latent3
## AL627309.1 78 28
## AL627309.5 49 8
## LINC01409 221 60
## LINC01128 202 87
## LINC00115 10 4
## FAM41C 54 21
colSums(pb1mf)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_latent1 mdm_latent2
## 30843590 22855559 26196614 13879736 35238765 39384235
## mdm_latent3
## 15058506
des1m <- as.data.frame(grepl("active",colnames(pb1mf)))
colnames(des1m) <- "case"
plot(cmdscale(dist(t(pb1mf))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb1mf))), labels=colnames(pb1mf) )
dds <- DESeqDataSetFromMatrix(countData = pb1mf , colData = des1m, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de1mf <- de
write.table(de1mf,"de1mf.tsv",sep="\t")
nrow(subset(de1mf,padj<0.05 & log2FoldChange>0))
## [1] 111
nrow(subset(de1mf,padj<0.05 & log2FoldChange<0))
## [1] 232
head(subset(de1mf,padj<0.05 & log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active MDM cells") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
LRRC25 | 389.10200 | 1.2420866 | 0.2192739 | 5.664542 | 0.0e+00 | 0.0000173 |
NMRK2 | 1221.29447 | 1.8111399 | 0.3328557 | 5.441216 | 1.0e-07 | 0.0000385 |
PHACTR1 | 363.65716 | 1.4735762 | 0.2807122 | 5.249421 | 2.0e-07 | 0.0000797 |
IL7R | 18.77671 | 4.1477458 | 0.8171411 | 5.075923 | 4.0e-07 | 0.0001699 |
DRAXIN | 478.82311 | 1.0070479 | 0.2036374 | 4.945299 | 8.0e-07 | 0.0002897 |
IL6R-AS1 | 79.77999 | 2.5938871 | 0.5445190 | 4.763630 | 1.9e-06 | 0.0005362 |
UBE2J1 | 2818.57491 | 0.4811650 | 0.1023425 | 4.701518 | 2.6e-06 | 0.0006742 |
HBEGF | 408.02932 | 0.8614193 | 0.1855476 | 4.642578 | 3.4e-06 | 0.0008299 |
NDFIP2 | 270.63116 | 0.9089767 | 0.1982830 | 4.584240 | 4.6e-06 | 0.0010361 |
DIRC3 | 378.78002 | 1.1228802 | 0.2496282 | 4.498211 | 6.9e-06 | 0.0014864 |
head(subset(de1mf,padj<0.05 & log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in active MDM cells") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
RCBTB2 | 819.53562 | -1.473879 | 0.1941941 | -7.589725 | 0 | 0.00e+00 |
TGFBI | 461.99975 | -2.449123 | 0.3481134 | -7.035418 | 0 | 0.00e+00 |
AL031123.1 | 525.20294 | -0.995166 | 0.1419224 | -7.012041 | 0 | 0.00e+00 |
HIST1H1C | 223.23751 | -2.154687 | 0.3283867 | -6.561431 | 0 | 2.00e-07 |
PTGDS | 28566.30351 | -2.275352 | 0.3586803 | -6.343678 | 0 | 6.00e-07 |
EVI2A | 562.82781 | -1.145883 | 0.1812752 | -6.321233 | 0 | 6.00e-07 |
MAP1A | 95.16198 | -2.073547 | 0.3294914 | -6.293175 | 0 | 6.00e-07 |
FLRT2 | 303.44509 | -1.893853 | 0.3124055 | -6.062163 | 0 | 2.40e-06 |
IGF1R | 384.50752 | -1.003939 | 0.1694912 | -5.923255 | 0 | 4.90e-06 |
CFD | 3010.50393 | -1.551909 | 0.2685750 | -5.778309 | 0 | 1.06e-05 |
m1m <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 15002
## Note: no. genes in output = 15002
## Note: estimated proportion of input genes in output = 1
mres1m <- mitch_calc(m1m,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres1m$enrichment_result
mitchtbl <- mres1m$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de1mf_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("MDM_latent_vs_active.html") ) {
mitch_report(mres1m,outfile="MDM_latent_vs_active.html")
}
Alv cells.
pbalv <- pbmono[,grep("alv",colnames(pbmono))]
pb1a <- pbalv[,c(grep("active",colnames(pbalv)),grep("latent",colnames(pbalv)))]
head(pb1a)
## alv_active1 alv_active2 alv_active3 alv_latent1 alv_latent2
## MIR1302-2HG 0 0 0 0 1
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 16 11 27 6 65
## AL627309.3 0 0 0 0 1
## AL627309.2 0 0 0 0 0
## alv_latent3 alv_latent4
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 59 76
## AL627309.3 1 0
## AL627309.2 0 1
pb1af <- pb1a[which(rowMeans(pb1a)>=10),]
head(pb1af)
## alv_active1 alv_active2 alv_active3 alv_latent1 alv_latent2
## AL627309.1 16 11 27 6 65
## AL627309.5 26 23 15 18 108
## LINC01409 87 151 92 34 133
## LINC01128 293 284 219 102 299
## LINC00115 10 7 10 5 14
## FAM41C 136 117 84 80 123
## alv_latent3 alv_latent4
## AL627309.1 59 76
## AL627309.5 90 67
## LINC01409 306 198
## LINC01128 413 348
## LINC00115 27 21
## FAM41C 134 150
colSums(pb1af)
## alv_active1 alv_active2 alv_active3 alv_latent1 alv_latent2 alv_latent3
## 30059559 28510626 24033240 16645696 43686559 56858514
## alv_latent4
## 52271268
des1a <- as.data.frame(grepl("active",colnames(pb1af)))
colnames(des1a) <- "case"
plot(cmdscale(dist(t(pb1af))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb1af))), labels=colnames(pb1af) )
dds <- DESeqDataSetFromMatrix(countData = pb1af , colData = des1a, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de1af <- de
write.table(de1af,"de1af.tsv",sep="\t")
nrow(subset(de1af,padj<0.05 & log2FoldChange>0))
## [1] 360
nrow(subset(de1af,padj<0.05 & log2FoldChange<0))
## [1] 259
head(subset(de1af,padj<0.05 & log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active Alv cells") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
LAYN | 646.06412 | 1.949699 | 0.1837933 | 10.608112 | 0 | 0 |
IL6R-AS1 | 496.13994 | 3.348015 | 0.3195846 | 10.476147 | 0 | 0 |
AL157912.1 | 876.82015 | 2.266161 | 0.2172491 | 10.431163 | 0 | 0 |
CDKN1A | 6332.25123 | 1.338772 | 0.1368473 | 9.782967 | 0 | 0 |
DNAJA4 | 337.39117 | 1.595379 | 0.1654746 | 9.641232 | 0 | 0 |
TNFSF15 | 504.64619 | 2.273073 | 0.2526606 | 8.996549 | 0 | 0 |
TCTEX1D2 | 5924.67714 | 1.301338 | 0.1452066 | 8.961980 | 0 | 0 |
IGFL2 | 135.86673 | 4.178248 | 0.4907575 | 8.513876 | 0 | 0 |
MID2 | 153.58862 | 2.106117 | 0.2485700 | 8.472936 | 0 | 0 |
CLIC6 | 76.68185 | 2.633839 | 0.3172527 | 8.302021 | 0 | 0 |
head(subset(de1af,padj<0.05 & log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active Alv cells") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
CEBPD | 942.22837 | -1.5655586 | 0.1606184 | -9.747070 | 0 | 0e+00 |
LYZ | 58220.50593 | -1.1385910 | 0.1231063 | -9.248847 | 0 | 0e+00 |
H1FX | 1153.38317 | -1.2336404 | 0.1591139 | -7.753189 | 0 | 0e+00 |
P2RY11 | 202.29568 | -1.2342577 | 0.1696748 | -7.274256 | 0 | 0e+00 |
JAKMIP2 | 75.00681 | -3.9911471 | 0.5561287 | -7.176662 | 0 | 0e+00 |
IFNGR2 | 2572.40169 | -0.6681364 | 0.0973720 | -6.861687 | 0 | 0e+00 |
CSF1R | 652.75459 | -1.7013376 | 0.2517604 | -6.757765 | 0 | 0e+00 |
SLCO3A1 | 245.66623 | -1.7535308 | 0.2756341 | -6.361806 | 0 | 1e-07 |
CD44 | 5417.83243 | -0.8186934 | 0.1309255 | -6.253126 | 0 | 2e-07 |
GCA | 1488.45196 | -0.9167194 | 0.1478308 | -6.201139 | 0 | 3e-07 |
m1a <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 16484
## Note: no. genes in output = 16484
## Note: estimated proportion of input genes in output = 1
mres1a <- mitch_calc(m1a,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres1a$enrichment_result
mitchtbl <- mres1a$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de1af_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("Alv_latent_vs_active.html") ) {
mitch_report(mres1a,outfile="Alv_latent_vs_active.html")
}
Combined.
mm1 <- merge(m1a,m1m,by=0)
head(mm1)
## Row.names x.x x.y
## 1 A1BG 1.0451580 0.1247177
## 2 A1BG-AS1 -0.6726339 -1.1348606
## 3 A2M 0.5824820 -0.5012939
## 4 A2M-AS1 -0.3548217 -1.0457070
## 5 A2ML1-AS1 -1.6525887 0.7798403
## 6 A4GALT 3.8981565 -0.7614755
rownames(mm1) <- mm1[,1]
mm1[,1]=NULL
colnames(mm1) <- c("Alv","MDM")
plot(mm1)
mylm <- lm(mm1)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0314 -0.8127 -0.0676 0.7340 10.1779
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.032175 0.010989 -2.928 0.00342 **
## MDM 0.442823 0.008061 54.936 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.334 on 14743 degrees of freedom
## Multiple R-squared: 0.1699, Adjusted R-squared: 0.1699
## F-statistic: 3018 on 1 and 14743 DF, p-value: < 2.2e-16
cor.test(mm1$Alv,mm1$MDM)
##
## Pearson's product-moment correlation
##
## data: mm1$Alv and mm1$MDM
## t = 54.936, df = 14743, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3987295 0.4255273
## sample estimates:
## cor
## 0.4122176
mm1r <- as.data.frame(apply(mm1,2,rank))
plot(mm1r,cex=0.3)
mylm <- lm(mm1r)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm1r)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9874.1 -3236.5 -26.6 3193.4 9902.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.479e+03 6.449e+01 69.45 <2e-16 ***
## MDM 3.925e-01 7.575e-03 51.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3915 on 14743 degrees of freedom
## Multiple R-squared: 0.1541, Adjusted R-squared: 0.154
## F-statistic: 2685 on 1 and 14743 DF, p-value: < 2.2e-16
cor.test(mm1r$Alv,mm1r$MDM)
##
## Pearson's product-moment correlation
##
## data: mm1r$Alv and mm1r$MDM
## t = 51.821, df = 14743, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3787906 0.4060997
## sample estimates:
## cor
## 0.3925316
MDM group.
pb2m <- pbmdm[,c(grep("bystander",colnames(pbmdm)),grep("latent",colnames(pbmdm)))]
head(pb2m)
## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1
## MIR1302-2HG 0 0 0 0
## FAM138A 0 0 0 0
## OR4F5 0 0 0 0
## AL627309.1 65 51 19 72
## AL627309.3 2 4 0 0
## AL627309.2 0 0 0 0
## mdm_latent2 mdm_latent3
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 78 28
## AL627309.3 2 2
## AL627309.2 0 0
pb2mf <- pb2m[which(rowMeans(pb2m)>=10),]
head(pb2mf)
## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1 mdm_latent2
## AL627309.1 65 51 19 72 78
## AL627309.5 33 35 8 35 49
## LINC01409 136 246 63 153 221
## LINC01128 216 208 111 215 202
## LINC00115 18 24 7 23 10
## FAM41C 43 69 19 62 54
## mdm_latent3
## AL627309.1 28
## AL627309.5 8
## LINC01409 60
## LINC01128 87
## LINC00115 4
## FAM41C 21
colSums(pb2mf)
## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1 mdm_latent2
## 40549303 36895352 16235047 35242448 39389780
## mdm_latent3
## 15060675
des2m <- as.data.frame(grepl("latent",colnames(pb2mf)))
colnames(des2m) <- "case"
plot(cmdscale(dist(t(pb2mf))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb2mf))), labels=colnames(pb2mf) )
dds <- DESeqDataSetFromMatrix(countData = pb2mf , colData = des2m, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de2mf <- de
#write.table(de2mf,"de2mf.tsv",sep="\t")
nrow(subset(de2mf,padj<0.05 & log2FoldChange>0))
## [1] 0
nrow(subset(de2mf,padj<0.05 & log2FoldChange<0))
## [1] 0
head(subset(de2mf,log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in latent compared to bystander (MDM unpaired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
ZFPM2 | 815.724490 | 1.0134336 | 0.3505885 | 2.890664 | 0.0038443 | 0.9997657 |
AC015660.1 | 13.134808 | 1.7822471 | 0.6580113 | 2.708536 | 0.0067581 | 0.9997657 |
PAEP | 8.352154 | 3.2437035 | 1.2087877 | 2.683435 | 0.0072870 | 0.9997657 |
MAP1A | 161.014052 | 0.5807312 | 0.2246662 | 2.584862 | 0.0097418 | 0.9997657 |
GRIP1 | 22.096324 | 1.3653641 | 0.5493578 | 2.485382 | 0.0129412 | 0.9997657 |
AL603840.1 | 64.653599 | 0.6972954 | 0.2851908 | 2.445013 | 0.0144847 | 0.9997657 |
LINC02257 | 59.540802 | 1.0232779 | 0.4203183 | 2.434531 | 0.0149111 | 0.9997657 |
FAM83G | 16.640778 | 1.2031135 | 0.5112349 | 2.353348 | 0.0186052 | 0.9997657 |
LINC02479 | 31.481351 | 1.0513244 | 0.4468126 | 2.352943 | 0.0186255 | 0.9997657 |
PTH2R | 9.387033 | 1.9614365 | 0.8560226 | 2.291337 | 0.0219439 | 0.9997657 |
head(subset(de2mf,log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in latent compared to bystander (MDM unpaired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
HSPA5 | 3129.492289 | -0.5464083 | 0.2136635 | -2.557330 | 0.0105479 | 0.9997657 |
AKR1C3 | 49.148938 | -0.8102721 | 0.3304780 | -2.451818 | 0.0142136 | 0.9997657 |
PLCB1 | 87.139407 | -2.1553080 | 0.8985866 | -2.398553 | 0.0164600 | 0.9997657 |
ZNF93 | 235.615217 | -0.4680771 | 0.2074674 | -2.256148 | 0.0240614 | 0.9997657 |
AC092142.1 | 9.750122 | -1.6296531 | 0.8050444 | -2.024302 | 0.0429391 | 0.9997657 |
ACTG1 | 17039.224704 | -0.2770873 | 0.1382883 | -2.003692 | 0.0451031 | 0.9997657 |
ZNF821 | 226.318656 | -0.3641906 | 0.1830925 | -1.989107 | 0.0466894 | 0.9997657 |
SMIM15-AS1 | 38.595801 | -0.7115038 | 0.3610366 | -1.970725 | 0.0487553 | 0.9997657 |
SYN2 | 12.534778 | -1.2709728 | 0.6477603 | -1.962104 | 0.0497504 | 0.9997657 |
PARTICL | 19.357134 | -1.0880757 | 0.5595985 | -1.944386 | 0.0518489 | 0.9997657 |
des2m$sample <- rep(1:3,2)
dds <- DESeqDataSetFromMatrix(countData = pb2mf , colData = des2m, design = ~ sample + case)
## 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))
de <- as.data.frame(zz[order(zz$pvalue),])
de2mf <- de
write.table(de2mf,"de2mf.tsv",sep="\t")
nrow(subset(de2mf,padj<0.05 & log2FoldChange>0))
## [1] 0
nrow(subset(de2mf,padj<0.05 & log2FoldChange<0))
## [1] 0
head(subset(de2mf,log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in latent compared to bystander (MDM paired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
ZFPM2 | 815.724490 | 1.0325654 | 0.2720485 | 3.795519 | 0.0001473 | 0.9999834 |
IL1RN | 254.031816 | 0.7031503 | 0.1985602 | 3.541246 | 0.0003982 | 0.9999834 |
PAEP | 8.352154 | 3.4470829 | 1.0303317 | 3.345605 | 0.0008210 | 0.9999834 |
MAP1A | 161.014052 | 0.5736414 | 0.1906963 | 3.008141 | 0.0026285 | 0.9999834 |
IGFBP6 | 155.859355 | 0.6138135 | 0.2083494 | 2.946078 | 0.0032183 | 0.9999834 |
AC015660.1 | 13.134808 | 1.7526778 | 0.6688347 | 2.620495 | 0.0087802 | 0.9999834 |
GRIP1 | 22.096324 | 1.3575888 | 0.5660267 | 2.398454 | 0.0164645 | 0.9999834 |
PRKN | 1022.174743 | 0.2546727 | 0.1065789 | 2.389522 | 0.0168703 | 0.9999834 |
LINC02257 | 59.540802 | 0.9626412 | 0.4034090 | 2.386266 | 0.0170204 | 0.9999834 |
NEO1 | 23.796659 | 1.1376615 | 0.4798766 | 2.370738 | 0.0177526 | 0.9999834 |
head(subset(de2mf,log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in latent compared to bystander (MDM paired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
HSPA5 | 3129.492289 | -0.5271919 | 0.1925378 | -2.738121 | 0.0061791 | 0.9999834 |
VSIG4 | 137.115837 | -0.5384120 | 0.2116862 | -2.543444 | 0.0109766 | 0.9999834 |
ZNF93 | 235.615217 | -0.4639859 | 0.1972344 | -2.352459 | 0.0186497 | 0.9999834 |
AKR1C3 | 49.148938 | -0.8055583 | 0.3468676 | -2.322380 | 0.0202125 | 0.9999834 |
ACTG1 | 17039.224704 | -0.2776672 | 0.1205940 | -2.302496 | 0.0213072 | 0.9999834 |
SOWAHD | 243.662376 | -0.4059598 | 0.1797295 | -2.258726 | 0.0239004 | 0.9999834 |
PLCB1 | 87.139407 | -1.5516291 | 0.7104441 | -2.184027 | 0.0289603 | 0.9999834 |
AC092142.1 | 9.750122 | -1.7293643 | 0.8532186 | -2.026871 | 0.0426756 | 0.9999834 |
ZNF821 | 226.318656 | -0.3630509 | 0.1824288 | -1.990097 | 0.0465803 | 0.9999834 |
CYTL1 | 716.355028 | -0.4856587 | 0.2455353 | -1.977959 | 0.0479334 | 0.9999834 |
m2m <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 15300
## Note: no. genes in output = 15300
## Note: estimated proportion of input genes in output = 1
mres2m <- mitch_calc(m2m,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres2m$enrichment_result
mitchtbl <- mres2m$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de2mf_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("MDM_bystander_vs_latent.html") ) {
mitch_report(mres2m,outfile="MDM_bystander_vs_latent.html")
}
Alv cells.
pb2a <- pbalv[,c(grep("bystander",colnames(pbalv)),grep("latent",colnames(pbalv)))]
head(pb2a)
## alv_bystander1 alv_bystander2 alv_bystander3 alv_bystander4
## MIR1302-2HG 0 0 0 0
## FAM138A 0 0 0 0
## OR4F5 0 0 0 0
## AL627309.1 20 38 9 14
## AL627309.3 0 1 0 0
## AL627309.2 0 0 0 0
## alv_latent1 alv_latent2 alv_latent3 alv_latent4
## MIR1302-2HG 0 1 0 0
## FAM138A 0 0 0 0
## OR4F5 0 0 0 0
## AL627309.1 6 65 59 76
## AL627309.3 0 1 1 0
## AL627309.2 0 0 0 1
pb2af <- pb2a[which(rowMeans(pb2a)>=10),]
head(pb2af)
## alv_bystander1 alv_bystander2 alv_bystander3 alv_bystander4
## AL627309.1 20 38 9 14
## AL627309.5 27 51 24 12
## LINC01409 58 83 75 73
## LINC01128 141 162 109 105
## LINC00115 1 7 5 2
## FAM41C 82 60 27 37
## alv_latent1 alv_latent2 alv_latent3 alv_latent4
## AL627309.1 6 65 59 76
## AL627309.5 18 108 90 67
## LINC01409 34 133 306 198
## LINC01128 102 299 413 348
## LINC00115 5 14 27 21
## FAM41C 80 123 134 150
colSums(pb2af)
## alv_bystander1 alv_bystander2 alv_bystander3 alv_bystander4 alv_latent1
## 20821170 22496343 14124595 13510489 16643623
## alv_latent2 alv_latent3 alv_latent4
## 43678024 56851234 52263981
des2a <- as.data.frame(grepl("latent",colnames(pb2af)))
colnames(des2a) <- "case"
plot(cmdscale(dist(t(pb2af))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb2af))), labels=colnames(pb2af) )
dds <- DESeqDataSetFromMatrix(countData = pb2af , colData = des2a, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de2af <- de
#write.table(de2af,"de2af.tsv",sep="\t")
nrow(subset(de2af,padj<0.05 & log2FoldChange>0))
## [1] 0
nrow(subset(de2af,padj<0.05 & log2FoldChange<0))
## [1] 0
head(subset(de2af, log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in latent compared to bystander (Alv unpaired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
IGFL2 | 6.831276 | 4.1473419 | 1.1419919 | 3.631674 | 0.0002816 | 0.9999275 |
IL6R-AS1 | 59.471366 | 0.9936364 | 0.3155831 | 3.148573 | 0.0016407 | 0.9999275 |
GUCY1A2 | 57.478671 | 1.1652601 | 0.3867694 | 3.012803 | 0.0025885 | 0.9999275 |
CLIC6 | 12.514507 | 1.3804460 | 0.5392894 | 2.559750 | 0.0104748 | 0.9999275 |
IL7R | 45.847087 | 1.0647822 | 0.4599362 | 2.315065 | 0.0206094 | 0.9999275 |
AC131254.1 | 20.634480 | 1.5538947 | 0.6839264 | 2.272020 | 0.0230853 | 0.9999275 |
CCL17 | 34.574222 | 1.5455779 | 0.6967303 | 2.218330 | 0.0265323 | 0.9999275 |
TLCD5 | 102.898040 | 0.5295161 | 0.2426812 | 2.181941 | 0.0291139 | 0.9999275 |
CXCL10 | 28.584967 | 4.3629570 | 2.0025608 | 2.178689 | 0.0293548 | 0.9999275 |
AC113404.1 | 13.483464 | 1.8436836 | 0.8522499 | 2.163313 | 0.0305171 | 0.9999275 |
head(subset(de2af, log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in latent compared to bystander (Alv unpaired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
PBK | 15.68054 | -1.5165527 | 0.5250931 | -2.888160 | 0.0038750 | 0.9999275 |
AL162586.1 | 23.37148 | -0.8275908 | 0.3459676 | -2.392105 | 0.0167521 | 0.9999275 |
CENPF | 110.08618 | -1.0023498 | 0.4539709 | -2.207960 | 0.0272470 | 0.9999275 |
TOP2A | 100.15582 | -0.9579134 | 0.4374134 | -2.189950 | 0.0285279 | 0.9999275 |
GTSE1 | 39.83766 | -0.9404875 | 0.4642714 | -2.025728 | 0.0427927 | 0.9999275 |
CCNA2 | 54.74930 | -0.7658193 | 0.3795303 | -2.017808 | 0.0436113 | 0.9999275 |
CDCA3 | 13.65163 | -1.0069050 | 0.5029028 | -2.002186 | 0.0452647 | 0.9999275 |
ASPM | 52.29393 | -1.0870301 | 0.5442275 | -1.997382 | 0.0457837 | 0.9999275 |
HJURP | 12.85119 | -1.4352668 | 0.7272403 | -1.973580 | 0.0484295 | 0.9999275 |
AC016590.1 | 10.67140 | -0.9703872 | 0.5041304 | -1.924873 | 0.0542452 | 0.9999275 |
des2a$sample <- rep(1:4,2)
dds <- DESeqDataSetFromMatrix(countData = pb2af , colData = des2a, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de2af <- de
write.table(de2af,"de2af.tsv",sep="\t")
nrow(subset(de2af,padj<0.05 & log2FoldChange>0))
## [1] 0
nrow(subset(de2af,padj<0.05 & log2FoldChange<0))
## [1] 0
head(subset(de2af, log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in latent compared to bystander (Alv paired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
IGFL2 | 6.831276 | 4.1473419 | 1.1419919 | 3.631674 | 0.0002816 | 0.9999275 |
IL6R-AS1 | 59.471366 | 0.9936364 | 0.3155831 | 3.148573 | 0.0016407 | 0.9999275 |
GUCY1A2 | 57.478671 | 1.1652601 | 0.3867694 | 3.012803 | 0.0025885 | 0.9999275 |
CLIC6 | 12.514507 | 1.3804460 | 0.5392894 | 2.559750 | 0.0104748 | 0.9999275 |
IL7R | 45.847087 | 1.0647822 | 0.4599362 | 2.315065 | 0.0206094 | 0.9999275 |
AC131254.1 | 20.634480 | 1.5538947 | 0.6839264 | 2.272020 | 0.0230853 | 0.9999275 |
CCL17 | 34.574222 | 1.5455779 | 0.6967303 | 2.218330 | 0.0265323 | 0.9999275 |
TLCD5 | 102.898040 | 0.5295161 | 0.2426812 | 2.181941 | 0.0291139 | 0.9999275 |
CXCL10 | 28.584967 | 4.3629570 | 2.0025608 | 2.178689 | 0.0293548 | 0.9999275 |
AC113404.1 | 13.483464 | 1.8436836 | 0.8522499 | 2.163313 | 0.0305171 | 0.9999275 |
head(subset(de2af, log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in latent compared to bystander (Alv paired)") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
PBK | 15.68054 | -1.5165527 | 0.5250931 | -2.888160 | 0.0038750 | 0.9999275 |
AL162586.1 | 23.37148 | -0.8275908 | 0.3459676 | -2.392105 | 0.0167521 | 0.9999275 |
CENPF | 110.08618 | -1.0023498 | 0.4539709 | -2.207960 | 0.0272470 | 0.9999275 |
TOP2A | 100.15582 | -0.9579134 | 0.4374134 | -2.189950 | 0.0285279 | 0.9999275 |
GTSE1 | 39.83766 | -0.9404875 | 0.4642714 | -2.025728 | 0.0427927 | 0.9999275 |
CCNA2 | 54.74930 | -0.7658193 | 0.3795303 | -2.017808 | 0.0436113 | 0.9999275 |
CDCA3 | 13.65163 | -1.0069050 | 0.5029028 | -2.002186 | 0.0452647 | 0.9999275 |
ASPM | 52.29393 | -1.0870301 | 0.5442275 | -1.997382 | 0.0457837 | 0.9999275 |
HJURP | 12.85119 | -1.4352668 | 0.7272403 | -1.973580 | 0.0484295 | 0.9999275 |
AC016590.1 | 10.67140 | -0.9703872 | 0.5041304 | -1.924873 | 0.0542452 | 0.9999275 |
m2a <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 15956
## Note: no. genes in output = 15956
## Note: estimated proportion of input genes in output = 1
mres2a <- mitch_calc(m2a,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres2a$enrichment_result
mitchtbl <- mres2a$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de2af_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("Alv_bystander_vs_latent.html") ) {
mitch_report(mres2a,outfile="Alv_bystander_vs_latent.html")
}
Combined.
mm2 <- merge(m2a,m2m,by=0)
head(mm2)
## Row.names x.x x.y
## 1 A1BG 0.144282200 -0.4751495
## 2 A1BG-AS1 0.428994559 0.6099293
## 3 A2M -0.252894005 0.4934169
## 4 A2M-AS1 0.161704527 -0.2630109
## 5 A2ML1-AS1 -0.647095488 -0.3808451
## 6 AAAS -0.009550646 -0.4414420
rownames(mm2) <- mm2[,1]
mm2[,1]=NULL
colnames(mm2) <- c("Alv","MDM")
plot(mm2)
mylm <- lm(mm2)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3098 -0.2762 -0.0051 0.2754 3.6223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004641 0.003933 -1.18 0.238
## MDM 0.082903 0.006446 12.86 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4781 on 14792 degrees of freedom
## Multiple R-squared: 0.01106, Adjusted R-squared: 0.01099
## F-statistic: 165.4 on 1 and 14792 DF, p-value: < 2.2e-16
cor.test(mm2$Alv,mm2$MDM)
##
## Pearson's product-moment correlation
##
## data: mm2$Alv and mm2$MDM
## t = 12.861, df = 14792, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08919709 0.12106940
## sample estimates:
## cor
## 0.1051602
mm2r <- as.data.frame(apply(mm2,2,rank))
plot(mm2r,cex=0.3)
mylm <- lm(mm2r)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm2r)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8172.8 -3619.7 46.7 3659.7 8139.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.600e+03 6.982e+01 94.53 <2e-16 ***
## MDM 1.078e-01 8.174e-03 13.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4246 on 14792 degrees of freedom
## Multiple R-squared: 0.01162, Adjusted R-squared: 0.01155
## F-statistic: 173.9 on 1 and 14792 DF, p-value: < 2.2e-16
cor.test(mm2r$Alv,mm2r$MDM)
##
## Pearson's product-moment correlation
##
## data: mm2r$Alv and mm2r$MDM
## t = 13.188, df = 14792, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.09184355 0.12369775
## sample estimates:
## cor
## 0.1077983
MDM group.
pb3m <- pbmdm[,c(grep("active",colnames(pbmdm)),grep("mock",colnames(pbmdm)))]
head(pb3m)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_mock1 mdm_mock2
## MIR1302-2HG 0 0 0 0 0 0
## FAM138A 0 0 0 0 0 0
## OR4F5 0 0 0 0 0 0
## AL627309.1 54 37 30 5 61 30
## AL627309.3 0 2 0 0 0 0
## AL627309.2 0 0 0 0 0 0
## mdm_mock3 mdm_mock4
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 6 15
## AL627309.3 0 0
## AL627309.2 0 0
pb3mf <- pb3m[which(rowMeans(pb3m)>=10),]
head(pb3mf)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_mock1 mdm_mock2
## AL627309.1 54 37 30 5 61 30
## AL627309.5 15 14 15 6 17 27
## LINC01409 129 114 104 27 119 116
## LINC01128 262 164 167 94 173 118
## FAM41C 49 39 63 68 61 25
## NOC2L 183 126 246 155 185 153
## mdm_mock3 mdm_mock4
## AL627309.1 6 15
## AL627309.5 1 16
## LINC01409 33 57
## LINC01128 57 159
## FAM41C 14 85
## NOC2L 75 221
colSums(pb3mf)
## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_mock1 mdm_mock2
## 30838055 22850851 26192488 13877828 29112212 21615267
## mdm_mock3 mdm_mock4
## 7549870 20729914
des3m <- as.data.frame(grepl("active",colnames(pb3mf)))
colnames(des3m) <- "case"
plot(cmdscale(dist(t(pb3mf))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb3mf))), labels=colnames(pb3mf) )
dds <- DESeqDataSetFromMatrix(countData = pb3mf , colData = des3m, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de3mf <- de
write.table(de3mf,"de3mf.tsv",sep="\t")
nrow(subset(de3mf,padj<0.05 & log2FoldChange>0))
## [1] 102
nrow(subset(de3mf,padj<0.05 & log2FoldChange<0))
## [1] 220
head(subset(de3mf,padj<0.05 & log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active MDM cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
PLAAT3 | 35.54458 | 2.6018024 | 0.4350913 | 5.979899 | 0.0e+00 | 0.0000031 |
LYRM1 | 207.51046 | 0.9010836 | 0.1550569 | 5.811308 | 0.0e+00 | 0.0000060 |
FAM241B | 98.31075 | 1.4123422 | 0.2609457 | 5.412399 | 1.0e-07 | 0.0000431 |
CDKN1A | 2611.64337 | 1.1118551 | 0.2071428 | 5.367578 | 1.0e-07 | 0.0000484 |
GRIP1 | 33.91692 | 2.0231034 | 0.3825670 | 5.288233 | 1.0e-07 | 0.0000641 |
DDB2 | 470.91041 | 0.6805458 | 0.1292762 | 5.264279 | 1.0e-07 | 0.0000706 |
KCNMB2-AS1 | 215.01389 | 1.5174273 | 0.2887412 | 5.255319 | 1.0e-07 | 0.0000716 |
FAS | 177.02499 | 0.9989147 | 0.1966111 | 5.080661 | 4.0e-07 | 0.0001505 |
MDM2 | 3009.65137 | 1.0736852 | 0.2190317 | 4.901962 | 9.0e-07 | 0.0003209 |
SFXN1 | 174.71064 | 0.8474226 | 0.1751359 | 4.838657 | 1.3e-06 | 0.0004321 |
head(subset(de1mf,padj<0.05 & log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in active MDM cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
RCBTB2 | 819.53562 | -1.473879 | 0.1941941 | -7.589725 | 0 | 0.00e+00 |
TGFBI | 461.99975 | -2.449123 | 0.3481134 | -7.035418 | 0 | 0.00e+00 |
AL031123.1 | 525.20294 | -0.995166 | 0.1419224 | -7.012041 | 0 | 0.00e+00 |
HIST1H1C | 223.23751 | -2.154687 | 0.3283867 | -6.561431 | 0 | 2.00e-07 |
PTGDS | 28566.30351 | -2.275352 | 0.3586803 | -6.343678 | 0 | 6.00e-07 |
EVI2A | 562.82781 | -1.145883 | 0.1812752 | -6.321233 | 0 | 6.00e-07 |
MAP1A | 95.16198 | -2.073547 | 0.3294914 | -6.293175 | 0 | 6.00e-07 |
FLRT2 | 303.44509 | -1.893853 | 0.3124055 | -6.062163 | 0 | 2.40e-06 |
IGF1R | 384.50752 | -1.003939 | 0.1694912 | -5.923255 | 0 | 4.90e-06 |
CFD | 3010.50393 | -1.551909 | 0.2685750 | -5.778309 | 0 | 1.06e-05 |
m3m <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 14592
## Note: no. genes in output = 14592
## Note: estimated proportion of input genes in output = 1
mres3m <- mitch_calc(m3m,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres3m$enrichment_result
mitchtbl <- mres3m$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de3mf_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("MDM_mock_vs_active.html") ) {
mitch_report(mres3m,outfile="MDM_mock_vs_active.html")
}
pb3a <- pbalv[,c(grep("active",colnames(pbalv)),grep("mock",colnames(pbalv)))]
head(pb3a)
## alv_active1 alv_active2 alv_active3 alv_mock1 alv_mock2 alv_mock3
## MIR1302-2HG 0 0 0 0 0 0
## FAM138A 0 0 0 0 0 0
## OR4F5 0 0 0 0 0 0
## AL627309.1 16 11 27 46 26 52
## AL627309.3 0 0 0 0 0 1
## AL627309.2 0 0 0 0 0 0
pb3af <- pb3a[which(rowMeans(pb3a)>=10),]
head(pb3af)
## alv_active1 alv_active2 alv_active3 alv_mock1 alv_mock2 alv_mock3
## AL627309.1 16 11 27 46 26 52
## AL627309.5 26 23 15 63 29 50
## LINC01409 87 151 92 59 103 140
## LINC01128 293 284 219 176 200 250
## LINC00115 10 7 10 10 12 15
## FAM41C 136 117 84 74 53 98
colSums(pb3af)
## alv_active1 alv_active2 alv_active3 alv_mock1 alv_mock2 alv_mock3
## 30051410 28505235 24028409 20446755 25545005 34873898
des3a <- as.data.frame(grepl("active",colnames(pb3af)))
colnames(des3a) <- "case"
plot(cmdscale(dist(t(pb3af))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb3af))), labels=colnames(pb3af) )
dds <- DESeqDataSetFromMatrix(countData = pb3af , colData = des3a, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de3af <- de
write.table(de3af,"de3af.tsv",sep="\t")
nrow(subset(de3af,padj<0.05 & log2FoldChange>0))
## [1] 891
nrow(subset(de3af,padj<0.05 & log2FoldChange<0))
## [1] 653
head(subset(de3af,padj<0.05 & log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active Alv cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
AL157912.1 | 754.26670 | 2.731970 | 0.2091846 | 13.060089 | 0 | 0 |
LAYN | 553.09496 | 2.261020 | 0.1871843 | 12.079108 | 0 | 0 |
HES4 | 382.20667 | 2.359190 | 0.2058619 | 11.460064 | 0 | 0 |
CDKN1A | 5303.58884 | 1.555197 | 0.1361511 | 11.422583 | 0 | 0 |
OCIAD2 | 353.43173 | 2.483107 | 0.2219754 | 11.186406 | 0 | 0 |
SDS | 4763.87880 | 2.107153 | 0.1953024 | 10.789184 | 0 | 0 |
IL6R-AS1 | 455.12951 | 3.520018 | 0.3307596 | 10.642225 | 0 | 0 |
CCL7 | 1444.11412 | 1.947067 | 0.1938562 | 10.043870 | 0 | 0 |
TNFRSF9 | 169.65948 | 2.509770 | 0.2514923 | 9.979509 | 0 | 0 |
CLIC6 | 64.19776 | 4.054497 | 0.4350349 | 9.319935 | 0 | 0 |
head(subset(de3af,padj<0.05 & log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active Alv cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
NDRG2 | 317.0287 | -1.8033036 | 0.1656273 | -10.887716 | 0 | 0 |
HIST1H1C | 791.4543 | -1.1579180 | 0.1105796 | -10.471351 | 0 | 0 |
CEBPD | 686.4608 | -1.4769083 | 0.1637149 | -9.021221 | 0 | 0 |
ADA2 | 1791.4889 | -1.0408494 | 0.1158943 | -8.981020 | 0 | 0 |
CRIM1 | 388.2713 | -1.6778340 | 0.1948415 | -8.611277 | 0 | 0 |
LYZ | 46011.7884 | -1.1761399 | 0.1385676 | -8.487841 | 0 | 0 |
RPL10A | 8121.1434 | -0.8692328 | 0.1046647 | -8.304927 | 0 | 0 |
CHST13 | 459.7071 | -1.6239334 | 0.1970679 | -8.240475 | 0 | 0 |
TGFBI | 395.6418 | -2.3302926 | 0.2874197 | -8.107631 | 0 | 0 |
PIK3CD | 583.8467 | -1.0404514 | 0.1382467 | -7.526050 | 0 | 0 |
m3a <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 15748
## Note: no. genes in output = 15748
## Note: estimated proportion of input genes in output = 1
mres3a <- mitch_calc(m3a,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres3a$enrichment_result
mitchtbl <- mres3a$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de3af_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("Alv_mock_vs_active.html") ) {
mitch_report(mres3a,outfile="Alv_mock_vs_active.html")
}
Combined.
mm3 <- merge(m3a,m3m,by=0)
head(mm3)
## Row.names x.x x.y
## 1 A1BG 1.2688786 0.37303646
## 2 A1BG-AS1 -0.7086023 -1.24450352
## 3 A2M -1.1312069 -0.45366896
## 4 A2ML1-AS1 -0.5102606 0.25614343
## 5 A4GALT 3.2131861 0.03284305
## 6 AAAS 0.8200706 -0.35943212
rownames(mm3) <- mm3[,1]
mm3[,1]=NULL
colnames(mm3) <- c("Alv","MDM")
plot(mm3)
mylm <- lm(mm3)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3766 -0.8804 -0.0974 0.7740 11.4780
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.093161 0.011884 7.839 4.85e-15 ***
## MDM 0.839716 0.009027 93.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.423 on 14332 degrees of freedom
## Multiple R-squared: 0.3765, Adjusted R-squared: 0.3764
## F-statistic: 8653 on 1 and 14332 DF, p-value: < 2.2e-16
cor.test(mm3$Alv,mm3$MDM)
##
## Pearson's product-moment correlation
##
## data: mm3$Alv and mm3$MDM
## t = 93.021, df = 14332, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6032502 0.6236681
## sample estimates:
## cor
## 0.6135617
mm3r <- as.data.frame(apply(mm3,2,rank))
plot(mm3r,cex=0.3)
mylm <- lm(mm3r)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm3r)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9964.7 -2495.7 -94.3 2415.6 11009.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.776e+03 5.464e+01 50.81 <2e-16 ***
## MDM 6.126e-01 6.602e-03 92.80 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3271 on 14332 degrees of freedom
## Multiple R-squared: 0.3753, Adjusted R-squared: 0.3753
## F-statistic: 8611 on 1 and 14332 DF, p-value: < 2.2e-16
cor.test(mm3r$Alv,mm3r$MDM)
##
## Pearson's product-moment correlation
##
## data: mm3r$Alv and mm3r$MDM
## t = 92.797, df = 14332, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6023124 0.6227672
## sample estimates:
## cor
## 0.6126424
MDM group.
pb4m <- pbmdm[,c(grep("mock",colnames(pbmdm)),grep("bystander",colnames(pbmdm)))]
head(pb4m)
## mdm_mock1 mdm_mock2 mdm_mock3 mdm_mock4 mdm_bystander1
## MIR1302-2HG 0 0 0 0 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 61 30 6 15 65
## AL627309.3 0 0 0 0 2
## AL627309.2 0 0 0 0 0
## mdm_bystander2 mdm_bystander3
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 51 19
## AL627309.3 4 0
## AL627309.2 0 0
pb4mf <- pb4m[which(rowMeans(pb4m)>=10),]
head(pb4mf)
## mdm_mock1 mdm_mock2 mdm_mock3 mdm_mock4 mdm_bystander1
## AL627309.1 61 30 6 15 65
## AL627309.5 17 27 1 16 33
## LINC01409 119 116 33 57 136
## LINC01128 173 118 57 159 216
## LINC00115 19 6 3 1 18
## FAM41C 61 25 14 85 43
## mdm_bystander2 mdm_bystander3
## AL627309.1 51 19
## AL627309.5 35 8
## LINC01409 246 63
## LINC01128 208 111
## LINC00115 24 7
## FAM41C 69 19
colSums(pb4mf)
## mdm_mock1 mdm_mock2 mdm_mock3 mdm_mock4 mdm_bystander1
## 29115660 21618680 7551139 20731699 40542425
## mdm_bystander2 mdm_bystander3
## 36888573 16232197
des4m <- as.data.frame(grepl("bystander",colnames(pb4mf)))
colnames(des4m) <- "case"
plot(cmdscale(dist(t(pb4mf))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb4mf))), labels=colnames(pb4mf) )
dds <- DESeqDataSetFromMatrix(countData = pb4mf , colData = des4m, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de4mf <- de
write.table(de4mf,"de4mf.tsv",sep="\t")
nrow(subset(de4mf,padj<0.05 & log2FoldChange>0))
## [1] 1
nrow(subset(de4mf,padj<0.05 & log2FoldChange<0))
## [1] 27
head(subset(de4mf, log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in bystander MDM cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
RBP7 | 3.371428e+02 | 0.9028624 | 0.2108153 | 4.282718 | 0.0000185 | 0.0182163 |
IFI6 | 1.236179e+04 | 1.0642894 | 0.2742746 | 3.880380 | 0.0001043 | 0.0514515 |
CCL8 | 9.256934e+00 | 3.2369960 | 0.8663313 | 3.736441 | 0.0001866 | 0.0863224 |
PLAAT3 | 2.037398e+01 | 1.5551917 | 0.4436038 | 3.505813 | 0.0004552 | 0.1727485 |
ISG15 | 9.101177e+02 | 0.7042685 | 0.2014310 | 3.496326 | 0.0004717 | 0.1745331 |
TNFAIP6 | 1.180878e+02 | 1.3809310 | 0.4155475 | 3.323161 | 0.0008900 | 0.2533181 |
PSME2 | 3.473007e+03 | 0.3948661 | 0.1219777 | 3.237198 | 0.0012071 | 0.3370754 |
HEXIM2 | 9.105406e+01 | 0.7473266 | 0.2461812 | 3.035676 | 0.0024000 | 0.5919922 |
HES1 | 7.898897e+01 | 1.8750551 | 0.6238906 | 3.005423 | 0.0026521 | 0.6330866 |
APOE | 2.885174e+05 | 0.3003852 | 0.1018483 | 2.949339 | 0.0031845 | 0.7141099 |
head(subset(de4mf, log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top downregulated genes in bystander MDM cells compared to mock") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
CDK1 | 60.01187 | -3.338466 | 0.5884689 | -5.673138 | 0.0e+00 | 0.0002075 |
CENPF | 126.75646 | -2.876893 | 0.5474866 | -5.254728 | 1.0e-07 | 0.0010970 |
CIT | 20.94293 | -2.835207 | 0.5566413 | -5.093419 | 4.0e-07 | 0.0017349 |
CEP55 | 46.62027 | -2.159231 | 0.4297160 | -5.024785 | 5.0e-07 | 0.0018648 |
CCL4L2 | 38.47946 | -4.011528 | 0.8118221 | -4.941388 | 8.0e-07 | 0.0022960 |
TYMS | 130.85145 | -2.198421 | 0.4486049 | -4.900574 | 1.0e-06 | 0.0023571 |
UBE2C | 70.95719 | -3.619243 | 0.7488856 | -4.832838 | 1.3e-06 | 0.0027504 |
RRM2 | 48.14325 | -2.669279 | 0.5545961 | -4.813015 | 1.5e-06 | 0.0027504 |
MKI67 | 150.04105 | -3.068534 | 0.6451095 | -4.756610 | 2.0e-06 | 0.0032374 |
TOP2A | 77.58402 | -2.914842 | 0.6184587 | -4.713074 | 2.4e-06 | 0.0036113 |
m4m <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 14825
## Note: no. genes in output = 14825
## Note: estimated proportion of input genes in output = 1
mres4m <- mitch_calc(m4m,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres4m$enrichment_result
mitchtbl <- mres4m$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de4mf_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("MDM_mock_vs_bystander.html") ) {
mitch_report(mres4m,outfile="MDM_mock_vs_bystander.html")
}
Alv cells.
pb4a <- pbalv[,c(grep("mock",colnames(pbalv)),grep("bystander",colnames(pbalv)))]
head(pb4a)
## alv_mock1 alv_mock2 alv_mock3 alv_bystander1 alv_bystander2
## MIR1302-2HG 0 0 0 0 0
## FAM138A 0 0 0 0 0
## OR4F5 0 0 0 0 0
## AL627309.1 46 26 52 20 38
## AL627309.3 0 0 1 0 1
## AL627309.2 0 0 0 0 0
## alv_bystander3 alv_bystander4
## MIR1302-2HG 0 0
## FAM138A 0 0
## OR4F5 0 0
## AL627309.1 9 14
## AL627309.3 0 0
## AL627309.2 0 0
pb4af <- pb4a[which(rowMeans(pb4a)>=10),]
head(pb4af)
## alv_mock1 alv_mock2 alv_mock3 alv_bystander1 alv_bystander2
## AL627309.1 46 26 52 20 38
## AL627309.5 63 29 50 27 51
## LINC01409 59 103 140 58 83
## LINC01128 176 200 250 141 162
## FAM41C 74 53 98 82 60
## SAMD11 19 44 37 1 29
## alv_bystander3 alv_bystander4
## AL627309.1 9 14
## AL627309.5 24 12
## LINC01409 75 73
## LINC01128 109 105
## FAM41C 27 37
## SAMD11 16 24
colSums(pb4af)
## alv_mock1 alv_mock2 alv_mock3 alv_bystander1 alv_bystander2
## 20441114 25539265 34865693 20814782 22486763
## alv_bystander3 alv_bystander4
## 14120665 13506372
des4a <- as.data.frame(grepl("bystander",colnames(pb4af)))
colnames(des4a) <- "case"
plot(cmdscale(dist(t(pb4af))), xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",pch=19,col="gray",cex=2)
text(cmdscale(dist(t(pb4af))), labels=colnames(pb4af) )
dds <- DESeqDataSetFromMatrix(countData = pb4af , colData = des4a, design = ~ case)
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
de <- as.data.frame(zz[order(zz$pvalue),])
de4af <- de
write.table(de4af,"de4af.tsv",sep="\t")
nrow(subset(de4af,padj<0.05 & log2FoldChange>0))
## [1] 4
nrow(subset(de4af,padj<0.05 & log2FoldChange<0))
## [1] 0
head(subset(de4af,padj<0.05 & log2FoldChange>0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in latent Alv cells compared to ") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
IFITM3 | 396.1914 | 1.8164917 | 0.3245822 | 5.596399 | 0.0e+00 | 0.0001983 |
EPSTI1 | 421.0338 | 2.9145326 | 0.5241180 | 5.560833 | 0.0e+00 | 0.0001983 |
IFI6 | 13043.5823 | 2.2745088 | 0.4392235 | 5.178477 | 2.0e-07 | 0.0011013 |
STAT1 | 2314.6957 | 0.8931956 | 0.1992453 | 4.482893 | 7.4e-06 | 0.0271889 |
head(subset(de4af, log2FoldChange<0),10)[,1:6] %>%
kbl(caption="Top upregulated genes in active Alv cells") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
SPP1 | 40458.66137 | -0.4122804 | 0.0965604 | -4.269664 | 0.0000196 | 0.0578258 |
PPBP | 21.33519 | -6.4579302 | 1.5515157 | -4.162337 | 0.0000315 | 0.0775391 |
F13A1 | 30.49474 | -5.3300819 | 1.2950413 | -4.115762 | 0.0000386 | 0.0814198 |
FCGBP | 20.90544 | -4.3156170 | 1.2150996 | -3.551657 | 0.0003828 | 0.3140987 |
MTSS1 | 195.93714 | -0.7031042 | 0.2021971 | -3.477320 | 0.0005065 | 0.3561808 |
SCN9A | 18.86558 | -1.2563392 | 0.4509299 | -2.786108 | 0.0053345 | 0.9998790 |
CTSV | 120.56641 | -0.9427725 | 0.3458440 | -2.726005 | 0.0064106 | 0.9998790 |
SLC9A3R2 | 32.00124 | -0.9907193 | 0.3969273 | -2.495972 | 0.0125613 | 0.9998790 |
PLD4 | 55.32855 | -0.8244766 | 0.3314254 | -2.487669 | 0.0128583 | 0.9998790 |
BIN2 | 739.68199 | -0.4605689 | 0.1929669 | -2.386776 | 0.0169968 | 0.9998790 |
m4a <- mitch_import(de,DEtype="deseq2",joinType="full")
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 15089
## Note: no. genes in output = 15089
## Note: estimated proportion of input genes in output = 1
mres4a <- mitch_calc(m4a,genesets=go,minsetsize=5,cores=4,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res <- mres4a$enrichment_result
mitchtbl <- mres4a$enrichment_result
goid <- sapply(strsplit(mitchtbl$set," "),"[[",1)
mysplit <- strsplit(mitchtbl$set," ")
mysplit <- lapply(mysplit, function(x) { x[2:length(x)] } )
godescription <- unlist(lapply(mysplit, function(x) { paste(x,collapse=" ") } ))
em <- data.frame(goid,godescription,mitchtbl$pANOVA,mitchtbl$p.adjustANOVA,sign(mitchtbl$s.dist),mitchtbl$s.dist)
colnames(em) <- c("GO.ID","Description","p.Val","FDR","Phenotype","ES")
write.table(em,"de4af_em.tsv",row.names = FALSE, quote=FALSE,sep="\t")
res <- subset(res,p.adjustANOVA<0.05)
resup <- subset(res,s.dist>0)
resdn <- subset(res,s.dist<0)
s <- c(head(resup$s.dist,10), head(resdn$s.dist,10))
names(s) <- c(head(resup$set,10),head(resdn$set,10))
s <- s[order(s)]
cols <- gsub("1","red",gsub("-1","blue",as.character(sign(s))))
par(mar=c(5,27,3,1))
barplot(abs(s),las=1,horiz=TRUE,col=cols,xlab="ES",cex.names=0.8,main="")
if (! file.exists("Alv_mock_vs_bystander.html") ) {
mitch_report(mres4a,outfile="Alv_mock_vs_bystander.html")
}
Combined.
mm4 <- merge(m4a,m4m,by=0)
head(mm4)
## Row.names x.x x.y
## 1 A1BG -0.07091926 0.31515557
## 2 A1BG-AS1 -0.45060890 -0.64197143
## 3 A2M -1.45649601 -0.27561959
## 4 A2ML1-AS1 1.53565784 -0.01805401
## 5 AAAS 0.58832021 -0.11873620
## 6 AACS 0.13939707 -0.80459468
rownames(mm4) <- mm4[,1]
mm4[,1]=NULL
colnames(mm4) <- c("Alv","MDM")
plot(mm4)
mylm <- lm(mm4)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4069 -0.4651 -0.0270 0.4554 5.4426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.099212 0.005777 17.175 < 2e-16 ***
## MDM 0.041525 0.006321 6.569 5.24e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6877 on 14276 degrees of freedom
## Multiple R-squared: 0.003014, Adjusted R-squared: 0.002944
## F-statistic: 43.15 on 1 and 14276 DF, p-value: 5.236e-11
cor.test(mm4$Alv,mm4$MDM)
##
## Pearson's product-moment correlation
##
## data: mm4$Alv and mm4$MDM
## t = 6.5691, df = 14276, p-value = 5.236e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03852895 0.07123594
## sample estimates:
## cor
## 0.05489717
mm4r <- as.data.frame(apply(mm4,2,rank))
plot(mm4r,cex=0.3)
mylm <- lm(mm4r)
abline(mylm,col="red",lty=2,lwd=3)
summary(mylm)
##
## Call:
## lm(formula = mm4r)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7325.0 -3551.3 12.5 3565.9 7363.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.886e+03 6.895e+01 99.86 < 2e-16 ***
## MDM 3.554e-02 8.364e-03 4.25 2.16e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4119 on 14276 degrees of freedom
## Multiple R-squared: 0.001263, Adjusted R-squared: 0.001193
## F-statistic: 18.06 on 1 and 14276 DF, p-value: 2.156e-05
cor.test(mm4r$Alv,mm4r$MDM)
##
## Pearson's product-moment correlation
##
## data: mm4r$Alv and mm4r$MDM
## t = 4.2495, df = 14276, p-value = 2.156e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01915191 0.05191631
## sample estimates:
## cor
## 0.03554366
l1 <- list("de1a"=de1af,"de1m"=de1mf,"de2a"=de2af,"de2m"=de2mf,
"de3a"=de3af,"de3m"=de3mf,"de4a"=de4af,"de4m"=de4mf)
lm <- mitch_import(x=l1,DEtype="deseq2",joinType="inner")
## Note: Mean no. genes in input = 15374.5
## Note: no. genes in output = 13992
## Note: estimated proportion of input genes in output = 0.91
lmres <- mitch_calc(x=lm,genesets=go,minsetsize=5,cores=8,priority="effect")
## Note: Enrichments with large effect sizes may not be
## statistically significant.
top <- head(lmres$enrichment_result,50)
top %>%
kbl(caption="Top pathways across all contrasts") %>%
kable_paper("hover", full_width = F)
set | setSize | pMANOVA | s.de1a | s.de1m | s.de2a | s.de2m | s.de3a | s.de3m | s.de4a | s.de4m | p.de1a | p.de1m | p.de2a | p.de2m | p.de3a | p.de3m | p.de4a | p.de4m | s.dist | SD | p.adjustMANOVA | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
2039 | GO:0019773 proteasome core complex, alpha-subunit complex | 7 | 0.0000005 | -0.0805455 | 0.3509168 | -0.3864447 | -0.8354972 | 0.2738342 | 0.7930231 | 0.8250983 | 0.8143317 | 0.7121416 | 0.1079160 | 0.0766589 | 0.0001290 | 0.2096693 | 0.0002794 | 0.0001564 | 0.0001904 | 1.739201 | 0.6141127 | 0.0000219 |
2051 | GO:0019864 IgG binding | 6 | 0.0003750 | -0.8217027 | -0.5944039 | -0.5871586 | -0.4254493 | -0.7562563 | -0.7475094 | 0.5545069 | 0.1266981 | 0.0004905 | 0.0116901 | 0.0127520 | 0.0711396 | 0.0013360 | 0.0015190 | 0.0186682 | 0.5910092 | 1.734502 | 0.4909392 | 0.0064502 |
4609 | GO:0070508 cholesterol import | 5 | 0.0017572 | 0.7392150 | 0.7706728 | -0.2950025 | 0.3536570 | 0.7544005 | 0.6239937 | 0.5328519 | -0.5624508 | 0.0042016 | 0.0028406 | 0.2533459 | 0.1708785 | 0.0034841 | 0.0156782 | 0.0390812 | 0.0294088 | 1.706237 | 0.5137242 | 0.0223181 |
2587 | GO:0032395 MHC class II receptor activity | 6 | 0.0000000 | -0.8632919 | -0.3760904 | 0.3394108 | -0.2494399 | -0.3413175 | -0.7069689 | 0.9591020 | -0.4419658 | 0.0002499 | 0.1106666 | 0.1499766 | 0.2900629 | 0.1477004 | 0.0027088 | 0.0000472 | 0.0608409 | 1.672020 | 0.5907155 | 0.0000008 |
3683 | GO:0045656 negative regulation of monocyte differentiation | 5 | 0.0080738 | -0.7485951 | -0.3788804 | 0.3954672 | 0.5263602 | -0.8770859 | -0.6009151 | -0.5030242 | -0.4858368 | 0.0037441 | 0.1423596 | 0.1256973 | 0.0415310 | 0.0006820 | 0.0199687 | 0.0514372 | 0.0599376 | 1.660375 | 0.5160364 | 0.0676118 |
3932 | GO:0048245 eosinophil chemotaxis | 7 | 0.0000511 | 0.6365085 | 0.7221104 | 0.5989172 | 0.4531488 | 0.6475407 | 0.7461770 | -0.1213034 | -0.4145564 | 0.0035421 | 0.0009375 | 0.0060692 | 0.0378899 | 0.0030085 | 0.0006286 | 0.5784183 | 0.0575364 | 1.628893 | 0.4338978 | 0.0012065 |
3271 | GO:0042613 MHC class II protein complex | 12 | 0.0000000 | -0.8022651 | -0.5341440 | 0.2022890 | -0.4797926 | -0.3812470 | -0.7281354 | 0.8564139 | -0.1764187 | 0.0000015 | 0.0013563 | 0.2250775 | 0.0040060 | 0.0222244 | 0.0000125 | 0.0000003 | 0.2900624 | 1.624857 | 0.5501007 | 0.0000000 |
374 | GO:0002503 peptide antigen assembly with MHC class II protein complex | 11 | 0.0000000 | -0.7876209 | -0.5446808 | 0.1782874 | -0.4826745 | -0.3406766 | -0.7168105 | 0.8480145 | -0.1272181 | 0.0000061 | 0.0017601 | 0.3059731 | 0.0055748 | 0.0504383 | 0.0000384 | 0.0000011 | 0.4651049 | 1.595925 | 0.5425071 | 0.0000000 |
2244 | GO:0030292 protein tyrosine kinase inhibitor activity | 5 | 0.0057415 | -0.8815186 | -0.6753271 | -0.2132409 | -0.3907772 | -0.8098806 | -0.5730035 | 0.2059770 | 0.1067992 | 0.0006404 | 0.0089188 | 0.4089930 | 0.1302467 | 0.0017106 | 0.0264981 | 0.4251371 | 0.6792243 | 1.571430 | 0.4078682 | 0.0524191 |
3335 | GO:0043009 chordate embryonic development | 5 | 0.0101461 | -0.4091657 | -0.5252163 | -0.3357260 | 0.7828555 | -0.6367198 | -0.6649746 | -0.4130836 | -0.4469722 | 0.1131168 | 0.0419759 | 0.1936174 | 0.0024320 | 0.0136779 | 0.0100225 | 0.1097073 | 0.0834973 | 1.544949 | 0.4644070 | 0.0784181 |
1457 | GO:0008541 proteasome regulatory particle, lid subcomplex | 8 | 0.0000735 | 0.0808066 | 0.5838995 | -0.5467499 | -0.6032966 | 0.2541655 | 0.7735090 | 0.6165082 | 0.5232408 | 0.6923097 | 0.0042383 | 0.0074095 | 0.0031278 | 0.2132352 | 0.0001513 | 0.0025309 | 0.0103866 | 1.525500 | 0.5309661 | 0.0016901 |
922 | GO:0006271 DNA strand elongation involved in DNA replication | 8 | 0.0000000 | -0.5730478 | -0.1114667 | -0.7970717 | 0.5430134 | -0.6506543 | -0.4310462 | -0.1739846 | -0.5926952 | 0.0050053 | 0.5851546 | 0.0000944 | 0.0078245 | 0.0014379 | 0.0347670 | 0.3941930 | 0.0036967 | 1.503897 | 0.4294211 | 0.0000013 |
3139 | GO:0038094 Fc-gamma receptor signaling pathway | 8 | 0.0010790 | -0.7986091 | -0.6071582 | -0.2996818 | -0.2346253 | -0.6835491 | -0.6791154 | 0.3622890 | 0.1619887 | 0.0000915 | 0.0029413 | 0.1422019 | 0.2505476 | 0.0008137 | 0.0008798 | 0.0760168 | 0.4276086 | 1.495729 | 0.4263174 | 0.0152976 |
2619 | GO:0032489 regulation of Cdc42 protein signal transduction | 5 | 0.0709729 | -0.7237149 | -0.5813827 | 0.2108672 | 0.0616144 | -0.7318653 | -0.7415600 | -0.4250947 | -0.0678773 | 0.0050696 | 0.0243674 | 0.4142279 | 0.8114408 | 0.0045948 | 0.0040827 | 0.0997590 | 0.7926912 | 1.476781 | 0.3885175 | 0.2684563 |
2508 | GO:0031726 CCR1 chemokine receptor binding | 6 | 0.0000043 | 0.5503599 | 0.5440917 | 0.6158063 | 0.4796225 | 0.5921159 | 0.4842462 | 0.2049430 | -0.5855856 | 0.0195697 | 0.0210045 | 0.0089970 | 0.0419115 | 0.0120166 | 0.0399744 | 0.3847119 | 0.0129935 | 1.475871 | 0.4030879 | 0.0001377 |
5422 | GO:1903238 positive regulation of leukocyte tethering or rolling | 6 | 0.0002540 | -0.8001335 | -0.6353020 | 0.3086896 | -0.1596835 | -0.6827065 | -0.6750798 | 0.0389199 | -0.2793508 | 0.0006879 | 0.0070412 | 0.1904322 | 0.4982287 | 0.0037791 | 0.0041876 | 0.8688866 | 0.2360749 | 1.471728 | 0.4010405 | 0.0047086 |
3095 | GO:0036150 phosphatidylserine acyl-chain remodeling | 7 | 0.0127078 | -0.6422698 | -0.6399408 | 0.2954901 | 0.1527044 | -0.6631697 | -0.7263599 | -0.4122682 | -0.2904643 | 0.0032536 | 0.0033675 | 0.1758352 | 0.4842095 | 0.0023776 | 0.0008743 | 0.0589281 | 0.1832994 | 1.467786 | 0.3935202 | 0.0907612 |
3281 | GO:0042719 mitochondrial intermembrane space protein transporter complex | 6 | 0.0051529 | 0.5746222 | 0.6736737 | 0.1501263 | 0.2599266 | 0.7952715 | 0.4638448 | 0.5239525 | -0.3720387 | 0.0147919 | 0.0042672 | 0.5242922 | 0.2702599 | 0.0007416 | 0.0491296 | 0.0262520 | 0.1145616 | 1.461053 | 0.3697564 | 0.0488673 |
510 | GO:0004185 serine-type carboxypeptidase activity | 5 | 0.0009746 | -0.4496032 | -0.7268607 | -0.0723672 | -0.6067205 | -0.2650604 | -0.2766712 | 0.5152070 | 0.7822264 | 0.0816953 | 0.0048813 | 0.7793224 | 0.0188027 | 0.3047414 | 0.2840446 | 0.0460434 | 0.0024517 | 1.458745 | 0.5314035 | 0.0142406 |
4686 | GO:0071139 resolution of DNA recombination intermediates | 5 | 0.0681896 | -0.7636663 | -0.2820476 | -0.2642883 | -0.0390791 | -0.7269465 | -0.7304640 | -0.0474870 | -0.5688282 | 0.0031031 | 0.2747908 | 0.3061524 | 0.8797293 | 0.0048763 | 0.0046734 | 0.8541178 | 0.0276185 | 1.456695 | 0.3064802 | 0.2619653 |
3655 | GO:0045569 TRAIL binding | 5 | 0.0008343 | 0.5551870 | 0.2297133 | 0.4938729 | 0.1994852 | 0.8454851 | 0.5646243 | 0.1587903 | 0.6388075 | 0.0315686 | 0.3737613 | 0.0558283 | 0.4398751 | 0.0010592 | 0.0287877 | 0.5386660 | 0.0133723 | 1.453148 | 0.2430014 | 0.0124773 |
4545 | GO:0070106 interleukin-27-mediated signaling pathway | 8 | 0.0000000 | -0.4510512 | -0.4420767 | 0.2545051 | -0.3206701 | 0.6387300 | 0.0834883 | 0.8810605 | 0.5833631 | 0.0271709 | 0.0303805 | 0.2126241 | 0.1163138 | 0.0017570 | 0.6826460 | 0.0000159 | 0.0042736 | 1.448450 | 0.5223175 | 0.0000000 |
811 | GO:0005839 proteasome core complex | 18 | 0.0000000 | -0.1652036 | 0.1315379 | -0.3506989 | -0.7031948 | 0.2096274 | 0.4878823 | 0.7592911 | 0.7359064 | 0.2250778 | 0.3340891 | 0.0100090 | 0.0000002 | 0.1237123 | 0.0003392 | 0.0000000 | 0.0000001 | 1.436009 | 0.5222824 | 0.0000000 |
4210 | GO:0051383 kinetochore organization | 5 | 0.0104385 | -0.4150854 | -0.3544577 | -0.5133481 | 0.4800887 | -0.4371059 | -0.6966612 | 0.3763352 | -0.6336312 | 0.1079968 | 0.1699119 | 0.0468342 | 0.0630283 | 0.0905437 | 0.0069803 | 0.1450599 | 0.0141415 | 1.418260 | 0.4487812 | 0.0795739 |
3113 | GO:0036402 proteasome-activating activity | 6 | 0.0023304 | 0.0461175 | 0.5064588 | -0.4429668 | -0.5022165 | 0.1682158 | 0.8227513 | 0.5707612 | 0.5154917 | 0.8449226 | 0.0316941 | 0.0602586 | 0.0331510 | 0.4755556 | 0.0004824 | 0.0154756 | 0.0287740 | 1.415551 | 0.4853636 | 0.0276216 |
3253 | GO:0042555 MCM complex | 11 | 0.0000003 | -0.5265913 | -0.1052922 | -0.2900495 | 0.1136933 | -0.6778875 | -0.6640636 | -0.4832077 | -0.6751565 | 0.0024940 | 0.5454647 | 0.0958158 | 0.5138801 | 0.0000989 | 0.0001368 | 0.0055226 | 0.0001055 | 1.405446 | 0.2944633 | 0.0000152 |
4955 | GO:0097100 supercoiled DNA binding | 6 | 0.0066141 | -0.7991325 | -0.4206588 | -0.3615520 | -0.1812050 | -0.7859288 | -0.4602936 | 0.3706325 | 0.1649268 | 0.0006987 | 0.0743809 | 0.1251444 | 0.4421532 | 0.0008559 | 0.0508907 | 0.1159384 | 0.4842278 | 1.404723 | 0.4155287 | 0.0579291 |
2620 | GO:0032494 response to peptidoglycan | 5 | 0.0699122 | -0.6597698 | -0.1965968 | -0.5409738 | 0.0590405 | -0.6087224 | -0.6217917 | 0.2737828 | -0.5850719 | 0.0106221 | 0.4465256 | 0.0361901 | 0.8191783 | 0.0184146 | 0.0160493 | 0.2891021 | 0.0234775 | 1.394460 | 0.3600882 | 0.2662736 |
848 | GO:0005942 phosphatidylinositol 3-kinase complex | 6 | 0.0252338 | -0.6262930 | -0.5649221 | 0.0452119 | 0.5048620 | -0.5719291 | -0.7129272 | 0.0726679 | -0.3550932 | 0.0078918 | 0.0165625 | 0.8479305 | 0.0322359 | 0.0152659 | 0.0024923 | 0.7579220 | 0.1320322 | 1.391061 | 0.4351368 | 0.1470388 |
5171 | GO:0140367 antibacterial innate immune response | 5 | 0.0037706 | -0.2229356 | -0.3465933 | 0.5975406 | 0.6729249 | -0.0896404 | -0.6593122 | -0.1035104 | -0.7079288 | 0.3880243 | 0.1795852 | 0.0206751 | 0.0091648 | 0.7285287 | 0.0106763 | 0.6885764 | 0.0061172 | 1.390797 | 0.5129721 | 0.0397019 |
2333 | GO:0030836 positive regulation of actin filament depolymerization | 5 | 0.0493005 | 0.2536784 | 0.8464860 | -0.2788732 | -0.0843498 | 0.4178308 | 0.6617716 | 0.5427468 | -0.3829699 | 0.3259787 | 0.0010448 | 0.2802293 | 0.7439718 | 0.1056855 | 0.0103878 | 0.0355838 | 0.1381018 | 1.385475 | 0.4521868 | 0.2183580 |
495 | GO:0004045 aminoacyl-tRNA hydrolase activity | 5 | 0.0358800 | 0.3955244 | 0.5036820 | -0.2125545 | -0.4395653 | 0.6315150 | 0.7064417 | 0.5412597 | 0.2227926 | 0.1256426 | 0.0511330 | 0.4105026 | 0.0887439 | 0.0144670 | 0.0062254 | 0.0360917 | 0.3883287 | 1.374410 | 0.4139049 | 0.1829010 |
2386 | GO:0031048 regulatory ncRNA-mediated heterochromatin formation | 7 | 0.0035132 | -0.5724603 | -0.2635375 | -0.3095255 | 0.4347209 | -0.6489504 | -0.5691098 | 0.1552173 | -0.6528117 | 0.0087217 | 0.2273159 | 0.1561927 | 0.0464143 | 0.0029459 | 0.0091227 | 0.4770470 | 0.0027805 | 1.370120 | 0.4037799 | 0.0376852 |
5517 | GO:1904948 midbrain dopaminergic neuron differentiation | 5 | 0.0625096 | -0.6873096 | -0.2812469 | -0.5601058 | 0.0144277 | -0.5597340 | -0.6172732 | 0.4095088 | -0.3505684 | 0.0077779 | 0.2761558 | 0.0300917 | 0.9554508 | 0.0302012 | 0.0168350 | 0.1128149 | 0.1746458 | 1.360274 | 0.3749702 | 0.2510905 |
903 | GO:0006172 ADP biosynthetic process | 5 | 0.0201749 | 0.6002288 | 0.3193966 | 0.5374848 | 0.5289054 | 0.6607135 | 0.6062344 | -0.0520340 | -0.1110317 | 0.0201106 | 0.2161911 | 0.0374089 | 0.0405555 | 0.0105111 | 0.0188980 | 0.8403287 | 0.6672611 | 1.360154 | 0.3062642 | 0.1260674 |
5238 | GO:0140948 histone H3K9 monomethyltransferase activity | 5 | 0.0412565 | -0.6338600 | -0.6695789 | -0.4182026 | -0.0031315 | -0.7124187 | -0.5421177 | -0.1074283 | -0.0231787 | 0.0141067 | 0.0095175 | 0.1053755 | 0.9903260 | 0.0058010 | 0.0357979 | 0.6774408 | 0.9284886 | 1.355923 | 0.2999000 | 0.1953033 |
1726 | GO:0014909 smooth muscle cell migration | 5 | 0.0978603 | 0.6563952 | 0.4643598 | 0.0810610 | -0.2764138 | 0.7266605 | 0.5528991 | 0.3417602 | 0.3892615 | 0.0110277 | 0.0721644 | 0.7536241 | 0.2844928 | 0.0048931 | 0.0322764 | 0.1857296 | 0.1317440 | 1.353336 | 0.3282024 | 0.3206772 |
5368 | GO:1902254 negative regulation of intrinsic apoptotic signaling pathway by p53 class mediator | 6 | 0.0042211 | 0.5770532 | 0.1011011 | 0.5082702 | 0.1512465 | 0.6712903 | 0.5179465 | 0.0758854 | 0.6892845 | 0.0143752 | 0.6680639 | 0.0310890 | 0.5212016 | 0.0044051 | 0.0280217 | 0.7475608 | 0.0034563 | 1.350628 | 0.2589746 | 0.0427940 |
5577 | GO:1990414 replication-born double-strand break repair via sister chromatid exchange | 6 | 0.0119313 | -0.3756375 | -0.4652986 | -0.6112303 | 0.5056247 | -0.6144954 | -0.5279089 | -0.2128080 | -0.3483722 | 0.1110966 | 0.0484235 | 0.0095211 | 0.0319762 | 0.0091444 | 0.0251398 | 0.3667345 | 0.1395077 | 1.345291 | 0.3648702 | 0.0868390 |
4446 | GO:0060992 response to fungicide | 6 | 0.0497383 | -0.6623767 | -0.6846370 | -0.0978598 | 0.2016064 | -0.5906383 | -0.6910005 | -0.1310358 | 0.0008103 | 0.0049577 | 0.0036816 | 0.6780987 | 0.3924987 | 0.0122318 | 0.0033763 | 0.5783667 | 0.9972577 | 1.342087 | 0.3625275 | 0.2189863 |
3270 | GO:0042612 MHC class I protein complex | 9 | 0.0000000 | -0.2064332 | -0.5217526 | -0.1754909 | -0.4302923 | 0.3301231 | -0.1629677 | 0.8378030 | 0.6505121 | 0.2836124 | 0.0067213 | 0.3620235 | 0.0254078 | 0.0863902 | 0.3972917 | 0.0000134 | 0.0007263 | 1.338435 | 0.5040533 | 0.0000001 |
3149 | GO:0038156 interleukin-3-mediated signaling pathway | 6 | 0.0150804 | -0.8707279 | -0.3100243 | 0.0131322 | -0.1334668 | -0.8502312 | -0.4150579 | 0.1197150 | -0.0780304 | 0.0002208 | 0.1885224 | 0.9555822 | 0.5713374 | 0.0003098 | 0.0783224 | 0.6116248 | 0.7406784 | 1.337107 | 0.3762780 | 0.1024677 |
491 | GO:0004017 adenylate kinase activity | 7 | 0.0024552 | 0.6206752 | 0.2577558 | 0.5443690 | 0.3783952 | 0.7255631 | 0.5853721 | 0.1314572 | 0.0739466 | 0.0044588 | 0.2376767 | 0.0126297 | 0.0830017 | 0.0008858 | 0.0073195 | 0.5470331 | 0.7347970 | 1.335308 | 0.2412153 | 0.0287593 |
3245 | GO:0042492 gamma-delta T cell differentiation | 5 | 0.0542480 | -0.5320798 | -0.4154286 | -0.6826196 | -0.2984915 | -0.5092014 | -0.5967112 | 0.2953457 | -0.1973976 | 0.0393660 | 0.1077057 | 0.0082079 | 0.2477752 | 0.0486396 | 0.0208520 | 0.2527941 | 0.4446762 | 1.323745 | 0.3103825 | 0.2331044 |
4689 | GO:0071162 CMG complex | 10 | 0.0000000 | -0.4103276 | 0.1139322 | -0.3473609 | 0.0933200 | -0.5041482 | -0.6597196 | -0.3704763 | -0.7840652 | 0.0246621 | 0.5327866 | 0.0571912 | 0.6094150 | 0.0057717 | 0.0003031 | 0.0425169 | 0.0000175 | 1.323678 | 0.3214564 | 0.0000001 |
2919 | GO:0035005 1-phosphatidylinositol-4-phosphate 3-kinase activity | 6 | 0.0243224 | -0.0530769 | -0.7187664 | 0.2353067 | 0.6745078 | -0.2369751 | -0.6149959 | -0.4994757 | -0.1865675 | 0.8218875 | 0.0022957 | 0.3182617 | 0.0042198 | 0.3148425 | 0.0090879 | 0.0341224 | 0.4287646 | 1.322283 | 0.4634371 | 0.1428913 |
2040 | GO:0019774 proteasome core complex, beta-subunit complex | 10 | 0.0000003 | -0.2968674 | 0.0192676 | -0.2923330 | -0.7058790 | 0.1167072 | 0.2804320 | 0.7055214 | 0.6892862 | 0.1040859 | 0.9159902 | 0.1094807 | 0.0001108 | 0.5228550 | 0.1246935 | 0.0001117 | 0.0001602 | 1.318091 | 0.4933941 | 0.0000143 |
5282 | GO:1900028 negative regulation of ruffle assembly | 5 | 0.0579430 | -0.6502181 | -0.4106813 | -0.0623865 | -0.3339530 | -0.5656252 | -0.7264603 | 0.2996926 | -0.2793594 | 0.0118059 | 0.1117880 | 0.8091232 | 0.1959812 | 0.0285056 | 0.0049049 | 0.2458774 | 0.2793916 | 1.312400 | 0.3362544 | 0.2406799 |
731 | GO:0005658 alpha DNA polymerase:primase complex | 5 | 0.0263270 | -0.4710803 | -0.1615643 | -0.6857368 | 0.5147208 | -0.6494173 | -0.3316937 | -0.0013012 | -0.4492314 | 0.0681361 | 0.5315953 | 0.0079198 | 0.0462491 | 0.0119102 | 0.1990238 | 0.9959801 | 0.0819480 | 1.310254 | 0.3950033 | 0.1511033 |
5375 | GO:1902425 positive regulation of attachment of mitotic spindle microtubules to kinetochore | 8 | 0.0000755 | -0.2920302 | -0.0800915 | -0.6364416 | 0.2548269 | -0.3103368 | -0.4766161 | 0.5589960 | -0.7132616 | 0.1526696 | 0.6948952 | 0.0018252 | 0.2120463 | 0.1285547 | 0.0195809 | 0.0061842 | 0.0004764 | 1.306285 | 0.4387112 | 0.0017228 |
colfunc <- colorRampPalette(c("blue", "white", "red"))
mx <- top[,grep("^s\\.",colnames(top))]
mx <- mx[,-ncol(mx)]
rownames(mx) <- top$set
heatmap.2(as.matrix(mx),scale="none",trace="none",margins=c(6,25),
col=colfunc(25),cexRow=0.6,cexCol=0.8)
For reproducibility.
save.image("scanalysis.Rdata")
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 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] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gplots_3.2.0 limma_3.60.6
## [3] SingleR_2.6.0 celldex_1.14.0
## [5] harmony_1.2.1 Rcpp_1.0.13-1
## [7] ggplot2_3.5.1 kableExtra_1.4.0
## [9] vioplot_0.5.0 zoo_1.8-12
## [11] sm_2.2-6.0 mitch_1.19.3
## [13] DESeq2_1.44.0 muscat_1.18.0
## [15] beeswarm_0.4.0 stringi_1.8.4
## [17] SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
## [19] Biobase_2.64.0 GenomicRanges_1.56.2
## [21] GenomeInfoDb_1.40.1 IRanges_2.38.1
## [23] S4Vectors_0.42.1 BiocGenerics_0.50.0
## [25] MatrixGenerics_1.16.0 matrixStats_1.4.1
## [27] hdf5r_1.3.11 Seurat_5.1.0
## [29] SeuratObject_5.0.2 sp_2.1-4
## [31] plyr_1.8.9
##
## loaded via a namespace (and not attached):
## [1] R.methodsS3_1.8.2 progress_1.2.3
## [3] goftest_1.2-3 HDF5Array_1.32.1
## [5] Biostrings_2.72.1 vctrs_0.6.5
## [7] spatstat.random_3.3-2 digest_0.6.37
## [9] png_0.1-8 corpcor_1.6.10
## [11] shape_1.4.6.1 gypsum_1.0.1
## [13] ggrepel_0.9.6 echarts4r_0.4.5
## [15] deldir_2.0-4 parallelly_1.39.0
## [17] MASS_7.3-61 reshape2_1.4.4
## [19] httpuv_1.6.15 foreach_1.5.2
## [21] withr_3.0.2 xfun_0.49
## [23] survival_3.7-0 memoise_2.0.1.9000
## [25] ggbeeswarm_0.7.2 systemfonts_1.1.0
## [27] GlobalOptions_0.1.2 gtools_3.9.5
## [29] pbapply_1.7-2 R.oo_1.27.0
## [31] prettyunits_1.2.0 GGally_2.2.1
## [33] KEGGREST_1.44.1 promises_1.3.1
## [35] httr_1.4.7 globals_0.16.3
## [37] fitdistrplus_1.2-1 rhdf5filters_1.16.0
## [39] rhdf5_2.48.0 rstudioapi_0.17.1
## [41] UCSC.utils_1.0.0 miniUI_0.1.1.1
## [43] generics_0.1.3 curl_6.0.1
## [45] zlibbioc_1.50.0 ScaledMatrix_1.12.0
## [47] polyclip_1.10-7 GenomeInfoDbData_1.2.12
## [49] ExperimentHub_2.12.0 SparseArray_1.4.8
## [51] xtable_1.8-4 stringr_1.5.1
## [53] doParallel_1.0.17 evaluate_1.0.1
## [55] S4Arrays_1.4.1 BiocFileCache_2.12.0
## [57] hms_1.1.3 irlba_2.3.5.1
## [59] colorspace_2.1-1 filelock_1.0.3
## [61] ROCR_1.0-11 reticulate_1.40.0
## [63] spatstat.data_3.1-4 magrittr_2.0.3
## [65] lmtest_0.9-40 later_1.4.0
## [67] viridis_0.6.5 lattice_0.22-6
## [69] spatstat.geom_3.3-4 future.apply_1.11.3
## [71] scattermore_1.2 scuttle_1.14.0
## [73] cowplot_1.1.3 RcppAnnoy_0.0.22
## [75] pillar_1.9.0 nlme_3.1-166
## [77] iterators_1.0.14 caTools_1.18.3
## [79] compiler_4.4.1 beachmat_2.20.0
## [81] RSpectra_0.16-2 tensor_1.5
## [83] minqa_1.2.8 crayon_1.5.3
## [85] abind_1.4-8 scater_1.32.1
## [87] blme_1.0-6 locfit_1.5-9.10
## [89] bit_4.5.0 dplyr_1.1.4
## [91] codetools_0.2-20 BiocSingular_1.20.0
## [93] bslib_0.8.0 alabaster.ranges_1.4.2
## [95] GetoptLong_1.0.5 plotly_4.10.4
## [97] remaCor_0.0.18 mime_0.12
## [99] splines_4.4.1 circlize_0.4.16
## [101] fastDummies_1.7.4 dbplyr_2.5.0
## [103] sparseMatrixStats_1.16.0 knitr_1.49
## [105] blob_1.2.4 utf8_1.2.4
## [107] clue_0.3-66 BiocVersion_3.19.1
## [109] lme4_1.1-35.5 listenv_0.9.1
## [111] DelayedMatrixStats_1.26.0 Rdpack_2.6.2
## [113] tibble_3.2.1 Matrix_1.7-1
## [115] statmod_1.5.0 svglite_2.1.3
## [117] fANCOVA_0.6-1 pkgconfig_2.0.3
## [119] tools_4.4.1 cachem_1.1.0
## [121] RhpcBLASctl_0.23-42 rbibutils_2.3
## [123] RSQLite_2.3.8 viridisLite_0.4.2
## [125] DBI_1.2.3 numDeriv_2016.8-1.1
## [127] fastmap_1.2.0 rmarkdown_2.29
## [129] scales_1.3.0 grid_4.4.1
## [131] ica_1.0-3 broom_1.0.7
## [133] AnnotationHub_3.12.0 sass_0.4.9
## [135] patchwork_1.3.0 BiocManager_1.30.25
## [137] ggstats_0.7.0 dotCall64_1.2
## [139] RANN_2.6.2 alabaster.schemas_1.4.0
## [141] farver_2.1.2 reformulas_0.4.0
## [143] aod_1.3.3 mgcv_1.9-1
## [145] yaml_2.3.10 cli_3.6.3
## [147] purrr_1.0.2 leiden_0.4.3.1
## [149] lifecycle_1.0.4 uwot_0.2.2
## [151] glmmTMB_1.1.10 mvtnorm_1.3-2
## [153] backports_1.5.0 BiocParallel_1.38.0
## [155] gtable_0.3.6 rjson_0.2.23
## [157] ggridges_0.5.6 progressr_0.15.1
## [159] jsonlite_1.8.9 edgeR_4.2.2
## [161] RcppHNSW_0.6.0 bitops_1.0-9
## [163] bit64_4.5.2 Rtsne_0.17
## [165] alabaster.matrix_1.4.2 spatstat.utils_3.1-1
## [167] BiocNeighbors_1.22.0 alabaster.se_1.4.1
## [169] jquerylib_0.1.4 spatstat.univar_3.1-1
## [171] R.utils_2.12.3 pbkrtest_0.5.3
## [173] lazyeval_0.2.2 alabaster.base_1.4.2
## [175] shiny_1.9.1 htmltools_0.5.8.1
## [177] sctransform_0.4.1 rappdirs_0.3.3
## [179] glue_1.8.0 spam_2.11-0
## [181] httr2_1.0.7 XVector_0.44.0
## [183] gridExtra_2.3 EnvStats_3.0.0
## [185] boot_1.3-31 igraph_2.1.1
## [187] variancePartition_1.34.0 TMB_1.9.15
## [189] R6_2.5.1 tidyr_1.3.1
## [191] labeling_0.4.3 cluster_2.1.6
## [193] Rhdf5lib_1.26.0 nloptr_2.1.1
## [195] DelayedArray_0.30.1 tidyselect_1.2.1
## [197] vipor_0.4.7 xml2_1.3.6
## [199] AnnotationDbi_1.66.0 future_1.34.0
## [201] rsvd_1.0.5 munsell_0.5.1
## [203] KernSmooth_2.23-24 data.table_1.16.2
## [205] htmlwidgets_1.6.4 ComplexHeatmap_2.20.0
## [207] RColorBrewer_1.1-3 rlang_1.1.4
## [209] spatstat.sparse_3.1-0 spatstat.explore_3.3-3
## [211] lmerTest_3.1-3 fansi_1.0.6