Single cell transcriptomes of macrophages with HIV 2024 - comparative analysis
Burnet Bioinformatics group
2025-01-06
Introduction
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")
})Load data
| 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")Make single cell experiment object
comb <- do.call(cbind,mylist)
sce <- SingleCellExperiment(list(counts=comb))
sce## class: SingleCellExperiment
## dim: 36622 28959
## metadata(0):
## assays(1): counts
## rownames(36622): HIV_LTRR HIV_LTRU5 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(28959): mdm_mock1|AAACGAATCACATACG mdm_mock1|AAACGCTCATCAGCGC
## ... react6|TTTGTTGTCTGAACGT react6|TTTGTTGTCTTGGCTC
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):Normalise data
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)## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')comb <- NormalizeData(comb)## Normalizing layer: countscomb <- FindVariableFeatures(comb, selection.method = "vst", nfeatures = 2000)## Finding variable features for layer countscomb <- ScaleData(comb)## Centering and scaling data matrixPCA and Cluster
comb <- RunPCA(comb, features = VariableFeatures(object = comb))## PC_ 1
## Positive: ARL15, NEAT1, DOCK3, EXOC4, MALAT1, FTX, LRMDA, RASAL2, DPYD, PLXDC2
## JMJD1C, TMEM117, FHIT, ZEB2, MITF, DENND4C, FRMD4B, ATG7, VPS13B, ZNF438
## MAML2, COP1, TPRG1, FMNL2, ATXN1, JARID2, ASAP1, ELMO1, FNDC3B, TCF12
## Negative: MT-ATP8, MTRNR2L1, FABP4, FABP3, HAMP, C1orf56, PVALB, GAPDH, AL136097.2, BIRC7
## S100A9, UTS2, MT1A, C1QC, MT1G, AC105402.3, AC025154.2, MT1H, CARD16, COX7A1
## BEX3, AL450311.1, PPP1R1A, LILRA1, MT2A, NUPR1, CCL23, FXYD2, HIV-Nef, CCL15
## PC_ 2
## Positive: CYSTM1, FAH, CD164, FDX1, GDF15, PSAT1, ATP6V1D, SAT1, TXN, BCAT1
## SLAMF9, CCPG1, HES2, HEBP2, CSTA, CTSL, S100A10, PHGDH, GARS, NUPR1
## TCEA1, STMN1, SCD, MARCO, RILPL2, SNHG32, SNHG5, NMB, TM4SF19, RETREG1
## Negative: SPRED1, RCBTB2, KCNMA1, GCLC, HLA-DRA, FGL2, MRC1, HLA-DRB1, HLA-DPA1, CD74
## SLCO2B1, HLA-DPB1, AC020656.1, RTN1, C1QA, PDGFC, MERTK, DPYD, LINC02345, ATP8B4
## CCDC102B, TGFBI, STAC, VAMP5, NFATC3, RNF128, MX1, ATG7, HLA-DRB5, XYLT1
## PC_ 3
## Positive: NCAPH, CRABP2, TM4SF19, GAL, DUSP2, RGCC, CHI3L1, CCL22, TRIM54, TMEM114
## ACAT2, RGS20, TCTEX1D2, NUMB, LINC02244, CCND1, IL1RN, MGST1, SERTAD2, AC092353.2
## GPC4, POLE4, CLU, ZNF366, SYNGR1, SLC20A1, GCLC, AC067751.1, OCSTAMP, MICAL3
## Negative: MARCKS, CD14, BTG1, TGFBI, MS4A6A, HLA-DQA1, FPR3, CTSC, CD163, MPEG1
## MEF2C, AIF1, TIMP1, HLA-DPB1, C1QC, HLA-DQB1, SELENOP, NUPR1, NAMPT, RNASE1
## HLA-DQA2, ALDH2, TCF4, EPB41L3, ARL4C, HIF1A, ZNF331, MS4A7, SAMSN1, CLEC4A
## PC_ 4
## Positive: MT-ATP8, MTRNR2L1, XIST, PDE4D, AL035446.2, HIV-Polp15p31, AC073359.2, AL365295.1, CCDC26, PRKCE
## RAD51B, CLVS1, EPHB1, C5orf17, HIV-LTRU5, LINC02320, FHIT, LINC01473, LY86-AS1, LINC01135
## FTX, STX18-AS1, DMGDH, PLEKHA5, AC079142.1, AF117829.1, KCNMB2-AS1, MIR646HG, AC008591.1, SLC22A15
## Negative: HLA-DRB1, CTSB, CD74, HSP90B1, TUBA1B, GAPDH, HLA-DPA1, ACTG1, HLA-DRA, IFI30
## RGCC, HSPA5, CYP1B1, AIF1, HLA-DPB1, C1QA, C15orf48, LYZ, CA2, FBP1
## TUBB, LGMN, S100A4, MMP9, TUBA1A, RNASE6, HLA-DQB1, CHI3L1, TXN, CCL3
## PC_ 5
## Positive: TYMS, PCLAF, TK1, MKI67, CENPM, MYBL2, RRM2, CLSPN, BIRC5, SHCBP1
## DIAPH3, CEP55, CDK1, CENPK, PRC1, ESCO2, CENPF, HELLS, TOP2A, NUSAP1
## NCAPG, TPX2, FAM111B, ANLN, CCNA2, ASPM, KIF11, CENPU, MAD2L1, HMMR
## Negative: PTGDS, LINC02244, MGST1, MMP9, LYZ, CLU, SYNGR1, IFI30, GCHFR, VAMP5
## CHI3L1, CRABP2, GCLC, RGCC, FABP3, CLEC12A, CSTA, NCAPH, SOD2, TMEM176B
## NCF1, FGL2, PLAUR, RARRES1, SPP1, MX1, ALDH2, S100A10, HLA-DRB1, RALAcomb <- RunHarmony(comb,"sample")## Transposing data matrix## Initializing state using k-means centroids initialization## Warning: Quick-TRANSfer stage steps exceeded maximum (= 1447950)## Harmony 1/10## Harmony 2/10## Harmony 3/10## Harmony 4/10## Harmony converged after 4 iterationsDimHeatmap(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 SNNcomb <- FindClusters(comb, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 28959
## Number of edges: 899987
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9036
## Number of communities: 14
## Elapsed time: 7 secondsUMAP
comb <- RunUMAP(comb, dims = 1:10)## 11:42:54 UMAP embedding parameters a = 0.9922 b = 1.112## 11:42:54 Read 28959 rows and found 10 numeric columns## 11:42:54 Using Annoy for neighbor search, n_neighbors = 30## 11:42:54 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:42:56 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd1ddf71ab
## 11:42:56 Searching Annoy index using 1 thread, search_k = 3000
## 11:43:06 Annoy recall = 100%
## 11:43:08 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:43:11 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:43:11 Commencing optimization for 200 epochs, with 1175456 positive edges
## 11:43:20 Optimization finishedDimPlot(comb, reduction = "umap")
Assign names to clusters
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 markersFeaturePlot(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.306472:0.325537:0.166567:... Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.311527:0.281167:0.189514:... Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.316271:0.276472:0.170736:... Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.214303:0.190773:0.121133:... Monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.290547:0.291036:0.154976:... Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.290907:0.279384:0.181857:... Monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.1823977 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 Monocytestable(pred_imm_broad$pruned.labels)##
## Basophils Dendritic cells Monocytes Neutrophils NK cells
## 4 655 24152 14 3
## Progenitors
## 2cellmetadata$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.1800566:0.485292:0.202974:... Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1959728:0.430960:0.226764:... Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1705944:0.441313:0.186890:... Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0958445:0.266677:0.121020:... Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1562428:0.415082:0.167816:... Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1850059:0.431679:0.205883:... Classical monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.0675290 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.1076301 Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1521533 Classical monocytestable(pred_imm_fine$pruned.labels)##
## Classical monocytes Intermediate monocytes
## 21584 3537
## Low-density basophils Low-density neutrophils
## 2 22
## Myeloid dendritic cells Naive B cells
## 466 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 1cellmetadata$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")
Make MDM object
mdmlist <- mylist[grep("mdm",names(mylist))]
comb1 <- do.call(cbind,mdmlist)
sce1 <- SingleCellExperiment(list(counts=comb1))
sce1## class: SingleCellExperiment
## dim: 36622 12416
## metadata(0):
## assays(1): counts
## rownames(36622): HIV_LTRR HIV_LTRU5 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(12416): 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)## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')comb1 <- NormalizeData(comb1)## Normalizing layer: countscomb1 <- FindVariableFeatures(comb1, selection.method = "vst", nfeatures = 2000)## Finding variable features for layer countscomb1 <- ScaleData(comb1)## Centering and scaling data matrixcomb1 <- RunPCA(comb1, features = VariableFeatures(object = comb1))## PC_ 1
## Positive: ARL15, NEAT1, DPYD, EXOC4, FTX, ZEB2, PLXDC2, LRMDA, DOCK3, RASAL2
## ATG7, JMJD1C, FHIT, TRIO, MALAT1, MAML2, VPS13B, TMEM117, MITF, ELMO1
## COP1, TCF12, DENND4C, JARID2, ATXN1, ARHGAP15, DOCK4, MEF2A, FNDC3B, FMNL2
## Negative: FABP4, MT-ATP8, FABP3, GAL, TXN, PRDX6, STMN1, S100A10, LILRA1, UTS2
## NUPR1, HIV-LTRU5, S100A9, NMB, SLAMF9, GSTM3, SEL1L2, S100P, GDF15, HAMP
## UCHL1, APOC2, MT1G, HIV-Nef, SLC35F1, SPP1, OCIAD2, MMP12, COX5B, MT1H
## PC_ 2
## Positive: TM4SF19, GPC4, ANO5, SPP1, TXNRD1, MGST1, C15orf48, ACAT2, TXN, FABP3
## PSD3, RGCC, CHI3L1, S100A10, RETREG1, RGS20, CRABP2, PRDX6, MREG, SCD
## COX5B, CCL22, TCTEX1D2, TGM2, CSF1, GAL, GDF15, COL8A2, SLC28A3, HES2
## Negative: TGFBI, HLA-DPB1, HLA-DRA, MS4A6A, HLA-DPA1, CEBPD, HLA-DQA1, HLA-DQB1, MPEG1, C1QA
## FPR3, CD74, RNASE1, ST8SIA4, CD163, C1QC, TCF4, CD14, LILRB2, EPB41L3
## ARL4C, PDE4B, MS4A7, AIF1, SEMA4A, HLA-DQA2, FOS, SIPA1L1, HLA-DRB5, HLA-DOA
## PC_ 3
## Positive: GSN, CYP1B1, CALR, S100A4, LPL, TIMP3, GCLC, DUSP2, ACTB, ID2
## NCAPH, LINC02345, CRIP1, CAMK1, IGSF6, OCSTAMP, LHFPL2, MRC1, FAIM2, RGCC
## STAC, MNDA, SPRED1, HLA-DRB1, CA2, CCND2, HSP90B1, AC020656.1, IL1RN, PDIA4
## Negative: NUPR1, XIST, PLEKHA5, PSAT1, G0S2, CCPG1, C5orf17, PRKCE, CCDC26, CLGN
## GRB10, BTG1, SUPV3L1, GDF15, AC073359.2, RETREG1, S100P, PHGDH, HES2, PDE4D
## STMN1, AL035446.2, BEX2, NIBAN1, DMGDH, ST6GALNAC3, MAML3, REV3L, EPHB1, SLC22A15
## PC_ 4
## Positive: IFI30, CTSL, BLVRB, MARCKS, SRGN, S100A10, LGMN, TXN, ANXA1, PLIN2
## TUBB, ACTG1, FUCA1, CD36, TUBA1C, GLIPR1, TPM4, PRDX6, PLEK, HSP90B1
## HLA-C, ADAMDEC1, TUBA1B, C15orf48, BTG1, GDF15, TMEM176A, TIMP1, SOD2, SLAMF9
## Negative: MT-ATP8, AC067751.1, LINC02408, OSBP2, TMEM117, KCNJ1, PLCL1, NFATC3, TPRG1, AC092353.2
## CT69, RCBTB2, LINC02345, ZNF366, AC008591.1, CPEB2, XYLT1, DNM3, STAC, ST5
## AL136317.2, LINC02320, MICAL3, VWA8, PKD1L1, NOS1AP, ZFPM2, AC093865.1, GCLC, CCDC102B
## PC_ 5
## Positive: PTGDS, CLU, LINC02244, BX664727.3, SYNGR1, RAB6B, COX5B, CSRP2, CCPG1, FOLR2
## HLA-C, NCAPH, RARRES1, RNASE6, LINC01010, MGST1, LYZ, SOD2, CPE, AL136317.2
## MT2A, AKR7A2, CAMP, PLEKHA5, MEIKIN, AC015660.2, HES2, CLEC12A, MT1E, PEBP4
## Negative: HIV-LTRU5, HIV-BaLEnv, HIV-Gagp17, HIV-Polprot, HIV-Polp15p31, HIV-Nef, HIV-TatEx1, HIV-Vif, HIV-Gagp1Pol, HIV-LTRR
## HIV-Gagp2p7, HIV-TatEx2Rev, HIV-EGFP, HIV-Vpu, HIV-EnvStart, CCL3, HIV-Vpr, CCL7, CSF1, SLC35F1
## SQLE, MSMO1, IL1RN, ACTB, TNFRSF9, TCTEX1D2, INSIG1, TPM4, LINC01091, C3comb1 <- 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 iterationsDimHeatmap(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 SNNcomb1 <- FindClusters(comb1, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 12416
## Number of edges: 395562
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8918
## Number of communities: 12
## Elapsed time: 1 secondscomb1 <- RunUMAP(comb1, dims = 1:10)## 11:48:17 UMAP embedding parameters a = 0.9922 b = 1.112## 11:48:17 Read 12416 rows and found 10 numeric columns## 11:48:17 Using Annoy for neighbor search, n_neighbors = 30## 11:48:17 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:48:18 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd168b857
## 11:48:18 Searching Annoy index using 1 thread, search_k = 3000
## 11:48:21 Annoy recall = 100%
## 11:48:22 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:48:24 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:48:24 Commencing optimization for 200 epochs, with 500376 positive edges
## 11:48:28 Optimization finishedDimPlot(comb1, reduction = "umap")
message("macrophage markers")## macrophage markersFeaturePlot(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.306472:0.325537:0.166567:... Monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.311527:0.281167:0.189514:... Monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.316271:0.276472:0.170736:... Monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.214303:0.190773:0.121133:... Monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.290547:0.291036:0.154976:... Monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.290907:0.279384:0.181857:... Monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.1823977 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 Monocytestable(pred_imm_broad1$pruned.labels)##
## Basophils Dendritic cells Monocytes Neutrophils
## 2 239 10015 2cellmetadata1$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.1800566:0.485292:0.202974:... Classical monocytes
## mdm_mock1|AAACGCTCATCAGCGC 0.1959728:0.430960:0.226764:... Classical monocytes
## mdm_mock1|AAACGCTGTCGAGTGA 0.1705944:0.441313:0.186890:... Classical monocytes
## mdm_mock1|AAAGAACGTGCAACAG 0.0958445:0.266677:0.121020:... Classical monocytes
## mdm_mock1|AAAGGTAAGCCATATC 0.1562428:0.415082:0.167816:... Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1850059:0.431679:0.205883:... Classical monocytes
## delta.next pruned.labels
## <numeric> <character>
## mdm_mock1|AAACGAATCACATACG 0.0675290 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.1076301 Classical monocytes
## mdm_mock1|AAAGGTAGTTTCCAAG 0.1521533 Classical monocytestable(pred_imm_fine1$pruned.labels)##
## Classical monocytes Intermediate monocytes
## 9249 1360
## Low-density neutrophils Myeloid dendritic cells
## 3 153
## Non classical monocytes Plasmacytoid dendritic cells
## 23 15cellmetadata1$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 monomono <- 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 countsmono <- 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, ZNF25-DT, NYX,
## GAPDH, LDHB, AL669970.1, TUBB4B, CYSTM1, PFN1, BCAT1, EPB41L1, GRM7, RERE,
## SAT1, RAP1GDS1, TFRC, CNIH3, FNIP2, FRMD4A, DDIT3, MRAP2, ME3, CD164, GPAT3,
## RHOF, JAZF1, MCTP1, CFD, FAH, SETD2, TCEA1, RABGEF1, SLC7A8, FARP1, STRBP, ENG,
## HIST1H2AC, NPTX1.## PC_ 1
## Positive: TXN, S100A10, FABP3, PRDX6, FABP4, COX5B, C15orf48, S100A9, BLVRB, PSME2
## SLAMF9, STMN1, GDF15, SPP1, LILRA1, NMB, GAL, TUBA1B, CTSL, HAMP
## UCHL1, CARD16, ANXA1, LILRA2, GSTM4, TUBA1A, NUPR1, IFI30, ACAT2, HES2
## Negative: ARL15, FTX, DPYD, EXOC4, VPS13B, NEAT1, LRMDA, DOCK3, ELMO1, ATG7
## FHIT, RASAL2, JMJD1C, TRIO, DOCK4, ZEB2, TMEM117, MAML2, MALAT1, COP1
## MBD5, ARHGAP15, TCF12, SPIDR, ATP9B, MED13L, ZSWIM6, RAD51B, TPRG1, SBF2
## PC_ 2
## Positive: TGFBI, HLA-DPB1, HLA-DRA, HLA-DPA1, CD74, MS4A6A, HLA-DQA1, HLA-DQB1, CEBPD, C1QA
## C1QC, MPEG1, FPR3, AIF1, CD14, CD163, MS4A7, ST8SIA4, RNASE1, LILRB2
## TCF4, ARL4C, HLA-DRB1, EPB41L3, FOS, TIMP1, HLA-DQA2, HLA-DRB5, SEMA4A, PDE4B
## Negative: TM4SF19, ANO5, GPC4, SPP1, PSD3, TXNRD1, RGS20, TCTEX1D2, RETREG1, MGST1
## ACAT2, CCL22, SLC28A3, MREG, CCDC26, FABP3, C15orf48, RGCC, CHI3L1, CRABP2
## LINC01135, SCD, GAL, CSF1, COL8A2, TXN, NUMB, NIBAN1, RAI14, GDF15
## PC_ 3
## Positive: NUPR1, CCPG1, BTG1, PSAT1, G0S2, XIST, GDF15, PLEKHA5, MARCKS, STMN1
## CLGN, SUPV3L1, PHGDH, HES2, NAMPT, S100P, C5orf17, RETREG1, GRB10, BEX2
## CCDC26, PRKCE, SLAMF9, ADAMDEC1, REV3L, SDS, NIBAN1, CXCR4, RAB38, B4GALT5
## Negative: GSN, GCLC, CYP1B1, DUSP2, LINC02345, TIMP3, S100A4, LPL, NCAPH, OCSTAMP
## CRIP1, STAC, CALR, FAIM2, CAMK1, ID2, SPRED1, MRC1, LHFPL2, LINC02408
## RCBTB2, ACTB, PLCL1, HIVEP3, CA2, AC020656.1, IGSF6, RNF128, CCND2, RGCC
## PC_ 4
## Positive: PTGDS, MT-ATP8, CLU, AL136317.2, BX664727.3, LY86-AS1, RCBTB2, AC008591.1, LINC01010, LINC02244
## KCNJ1, PKD1L1, EPHB1, AC015660.2, MT1G, AC067751.1, SPON2, TMEM117, AP000331.1, CT69
## OSBP2, NCAPH, FAIM2, MT1X, FLRT2, XYLT1, KCNMA1, KLRD1, MT1H, LINC02320
## Negative: ACTB, ACTG1, TPM4, CSF1, TUBA1B, CCL3, GLIPR1, CD36, TUBB, MGLL
## INSIG1, MSMO1, LDHA, SLC35F1, LGMN, PHLDA1, MARCKS, PLEK, CALR, CCL7
## BLVRB, CTSL, DHCR24, C15orf48, SQLE, HIV-Gagp17, DUSP6, TUBA1C, C3, LIMA1
## PC_ 5
## Positive: HIV-LTRU5, HIV-Polp15p31, HIV-BaLEnv, HIV-Gagp17, HIV-Polprot, HIV-Nef, HIV-TatEx1, HIV-Vif, HIV-Gagp1Pol, HIV-LTRR
## HIV-Gagp2p7, HIV-TatEx2Rev, HIV-Vpu, HIV-EGFP, HIV-EnvStart, HIV-Vpr, TCTEX1D2, SPRED1, MT-ATP8, MADD
## IL1RN, TNFRSF9, CCL3, ANK2, CCL7, PPARG, DTNA, PLCL1, PRKCA, LINC01091
## Negative: PCLAF, TYMS, CLU, TK1, PTGDS, SYNGR1, STMN1, RAD51AP1, COX5B, CENPM
## MKI67, BX664727.3, BIRC5, MYBL2, SRGN, LINC02244, CSRP2, RAB6B, CEP55, HLA-C
## CENPK, CENPF, KIF11, PRC1, DIAPH3, RRM2, FOLR2, SHCBP1, CLSPN, CCPG1DimHeatmap(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 SNNmono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)## 11:49:15 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:49:15 Read 10590 rows and found 4 numeric columns
## 11:49:15 Using Annoy for neighbor search, n_neighbors = 30
## 11:49:15 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:49:16 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd2ff29ac
## 11:49:16 Searching Annoy index using 1 thread, search_k = 3000
## 11:49:19 Annoy recall = 100%
## 11:49:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:49:22 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:49:22 Commencing optimization for 200 epochs, with 356498 positive edges
## 11:49:26 Optimization finishedDimPlot(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)
Make Alv object
alvlist <- mylist[grep("alv",names(mylist))]
comb1 <- do.call(cbind,alvlist)
sce1 <- SingleCellExperiment(list(counts=comb1))
sce1## class: SingleCellExperiment
## dim: 36622 13445
## metadata(0):
## assays(1): counts
## rownames(36622): HIV_LTRR HIV_LTRU5 ... AC007325.4 AC007325.2
## rowData names(0):
## colnames(13445): alv_mock1|AAACCCAGTGCTGCAC alv_mock1|AAAGAACCATTGAGGG
## ... alv_bystander4|TTTGATCTCGGCTTGG alv_bystander4|TTTGGTTAGACTGTTC
## 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)## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')comb1 <- NormalizeData(comb1)## Normalizing layer: countscomb1 <- FindVariableFeatures(comb1, selection.method = "vst", nfeatures = 2000)## Finding variable features for layer countscomb1 <- ScaleData(comb1)## Centering and scaling data matrixcomb1 <- RunPCA(comb1, features = VariableFeatures(object = comb1))## PC_ 1
## Positive: FABP4, NUPR1, HIV-Nef, HIV-LTRU5, TIMP1, S100A8, STMN1, AC105402.3, LILRA1, VSIG4
## CAMP, RBP1, ORM1, PRG2, CD79A, CCL5, MDK, UTS2, PDPN, AC078850.1
## RBP4, VMO1, AC026369.3, LGALS2, HLA-DOA, PCSK1N, NCF1, MS4A6E, PTGES, HCFC1-AS1
## Negative: DOCK3, RASAL2, ARL15, MALAT1, NEAT1, PLXDC2, ASAP1, LRMDA, TMEM117, DPYD
## EXOC4, MITF, ATG7, FTX, JMJD1C, FRMD4B, FMNL2, UBE2E2, ZNF438, MAML2
## FHIT, LPP, TPRG1, CHD9, DENND4C, ZEB2, RAPGEF1, ELMO1, LRCH1, VPS13B
## PC_ 2
## Positive: GAL, CYSTM1, TM4SF19, CSTB, ATP6V1D, CCL22, FDX1, GM2A, TGM2, SCD
## ACAT2, CD164, LGALS1, FAH, CSF1, RHOF, ELOC, PRDX6, CIR1, CYTOR
## FCMR, GAPLINC, C15orf48, UPP1, GOLGA7B, SNHG32, NCAPH, BCAT1, DCSTAMP, TNFRSF12A
## Negative: SAMSN1, TRPS1, SLC8A1, TGFBI, PDGFC, HLA-DPB1, MRC1, CAMK1D, SPRED1, AL356124.1
## SLC9A9, RCBTB2, HLA-DPA1, RTN1, FOS, CD14, HLA-DRA, ARHGAP15, AIF1, TET2
## C1QC, ATP8B4, FCHSD2, C1QA, ME1, XYLT1, VAMP5, MS4A6A, SLCO2B1, FGL2
## PC_ 3
## Positive: IFI30, CTSZ, AIF1, MARCKS, LGMN, CD74, HLA-DRA, CTSL, HLA-DPA1, CTSB
## HLA-DRB1, C1QA, ACTG1, HLA-DPB1, BLVRB, FUCA1, GAPDH, FPR3, HLA-DQB1, TUBA1B
## C1QC, CTSC, HLA-DQA1, TUBB, FBP1, FCGR3A, IL4I1, C1QB, MIF, TPM4
## Negative: KCNJ1, AC067751.1, C2orf92, LINC02320, AC015660.2, RGS20, AL136317.2, CT69, XIST, AC008591.1
## NOS1AP, NCAPH, AP000331.1, PDE4D, GDA, MICAL3, LINC01010, LY86-AS1, AC012668.3, CCL22
## AC092353.2, TRIM54, ANO5, BX664727.3, MIR646HG, TCTEX1D2, AC092944.1, SERTAD2, OSBP2, TDRD3
## PC_ 4
## Positive: FMN1, SNCA, HIV-Polp15p31, MARCO, SLC11A1, HIV-Polprot, HIV-BaLEnv, XIST, HIV-LTRU5, HAMP
## HIV-Vif, HIV-Nef, HIV-Gagp17, FRMD4A, CLMP, BCAT1, HIV-TatEx1, SASH1, ME3, PHACTR1
## CCDC26, SLA, AC023282.1, S100A9, DLEU1, HIV-LTRR, HIV-Gagp1Pol, TTC39B, SLC16A10, GYPC
## Negative: PTGDS, CRIP1, GCLC, LINC02244, MMP9, CLU, ALOX5AP, MX1, TMEM176A, PLEK
## TUBA1A, FGL2, SYNGR1, DUSP2, RGCC, IFIT3, VAMP5, KCNMA1, NCAPH, MMP7
## MX2, RNF128, MT2A, PRDX6, C1QTNF4, LYZ, FCMR, FAIM2, ISG15, IFIT1
## PC_ 5
## Positive: TYMS, PCLAF, TK1, CLSPN, MKI67, RRM2, CENPM, DIAPH3, MYBL2, SHCBP1
## CDK1, BIRC5, ESCO2, CEP55, CENPK, NCAPG, FAM111B, NUSAP1, HELLS, TOP2A
## CENPF, ATAD2, PRC1, DTL, TPX2, CCNA2, HMMR, KIF11, MCM10, ASPM
## Negative: CTSB, CD74, SAT1, IFI30, CTSZ, HLA-DRB1, RGCC, CFD, FUCA1, HLA-DRA
## HLA-DPA1, HLA-DPB1, C15orf48, HIV-BaLEnv, HSPA5, MMP9, ISG15, HIV-Gagp17, HIV-Polprot, ATP1B1
## LYZ, HIV-LTRU5, HIV-Nef, PLEK, GCLC, VAMP5, TNFSF13B, HIV-TatEx1, HIV-Polp15p31, FGL2comb1 <- 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 iterationsDimHeatmap(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 SNNcomb1 <- FindClusters(comb1, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 13445
## Number of edges: 425364
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8875
## Number of communities: 12
## Elapsed time: 1 secondscomb1 <- RunUMAP(comb1, dims = 1:10)## 11:53:25 UMAP embedding parameters a = 0.9922 b = 1.112## 11:53:25 Read 13445 rows and found 10 numeric columns## 11:53:25 Using Annoy for neighbor search, n_neighbors = 30## 11:53:25 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:53:26 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd6bbe7d2a
## 11:53:26 Searching Annoy index using 1 thread, search_k = 3000
## 11:53:30 Annoy recall = 100%
## 11:53:30 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:53:33 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:53:33 Commencing optimization for 200 epochs, with 544452 positive edges
## 11:53:37 Optimization finishedDimPlot(comb1, reduction = "umap")
message("macrophage markers")## macrophage markersFeaturePlot(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.273560:0.2723439:0.1612425:... Monocytes
## alv_mock1|AAAGAACCATTGAGGG 0.125010:0.0949686:0.0761100:... Dendritic cells
## alv_mock1|AAAGAACTCTCATGCC 0.093777:0.0888391:0.0837904:... Dendritic cells
## alv_mock1|AAAGGATAGCATGAAT 0.318711:0.3073772:0.1916631:... Monocytes
## alv_mock1|AAAGGATAGTCAGGGT 0.294485:0.2756727:0.1829141:... Monocytes
## alv_mock1|AAAGGATAGTTCCGGC 0.294678:0.2947249:0.1839451:... 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.10464548 Monocytes
## alv_mock1|AAAGGATAGTCAGGGT 0.14001755 Monocytes
## alv_mock1|AAAGGATAGTTCCGGC 0.12840710 Monocytestable(pred_imm_broad1$pruned.labels)##
## Basophils Dendritic cells Monocytes Neutrophils NK cells
## 2 337 11488 10 2
## Progenitors
## 1cellmetadata1$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.1630174:0.398232:0.1766975:...
## alv_mock1|AAAGAACCATTGAGGG 0.0698896:0.151804:0.0746815:...
## alv_mock1|AAAGAACTCTCATGCC 0.0614206:0.101479:0.0761397:...
## alv_mock1|AAAGGATAGCATGAAT 0.1802343:0.435146:0.2087306:...
## alv_mock1|AAAGGATAGTCAGGGT 0.1699671:0.389207:0.1877512:...
## alv_mock1|AAAGGATAGTTCCGGC 0.1664617:0.422480:0.1894660:...
## 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.21392e-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 monocytestable(pred_imm_fine1$pruned.labels)##
## Classical monocytes Intermediate monocytes
## 10164 1558
## Low-density basophils Low-density neutrophils
## 1 14
## Myeloid dendritic cells Naive B cells
## 257 1
## Natural killer cells Non classical monocytes
## 3 31
## Plasmacytoid dendritic cells Progenitor cells
## 8 4
## Terminal effector CD4 T cells
## 1cellmetadata1$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 monomono <- 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 countsmono <- RunPCA(mono, features = VariableFeatures(object = mono))## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: MGST1, NR1H3, GLIPR1, ACTB, IGSF6, IFI6,
## SSH2, DOCK10, EMP1, RPS20, AP003718.1, MOSPD1, TRIO, DOP1B, C11orf45, SNX10,
## PTPRJ, AL162411.1, SCFD1, SDC2, S100A6, NABP1, JAML, CRABP2, QKI, AC008443.9,
## MGAT5, ABCA1, UBASH3B, ZFAND3, CDKAL1, ARHGEF3, ASPH, C1orf21, DENND1A,
## PPP1R12B, ALDH1A1, DOCK2, LGALS3, LPAR1, AFF1, HMOX1, LINC01800, FAM89A,
## COL8A2, SEPTIN14, RALGAPA2, NREP-AS1, PITPNC1.## PC_ 1
## Positive: DOCK3, ARL15, MALAT1, TMEM117, FTX, EXOC4, RASAL2, LRMDA, FHIT, DPYD
## TPRG1, VPS13B, JMJD1C, NEAT1, ELMO1, ATG7, MAML2, ZNF438, COP1, MITF
## TTC28, FMNL2, VTI1A, MED13L, ERC1, ATP9B, ASAP1, SPIDR, DENND4C, ZSWIM6
## Negative: GAPDH, LGALS1, MIF, CSTB, BLVRB, TUBA1B, PRDX6, IFI30, FABP4, TXN
## TMEM176A, FAH, TUBA1A, EMP3, CD74, TIMP1, S100A9, HLA-DPA1, MMP9, NUPR1
## C15orf48, ELOC, HLA-DRB1, CFD, H2AFZ, LILRA1, PTGDS, CYSTM1, HLA-DRA, VAMP5
## PC_ 2
## Positive: GAL, CCL22, TM4SF19, CYSTM1, ATP6V1D, FDX1, SCD, TGM2, ACAT2, CIR1
## RHOF, GM2A, CSF1, IARS, DUSP13, GOLGA7B, NCAPH, BCAT1, CSTB, CSRP2
## FAH, SLC20A1, TCTEX1D2, SNHG32, NMB, EPB41L1, CD164, CYTOR, DCSTAMP, FCMR
## Negative: HLA-DPA1, HLA-DRA, HLA-DPB1, AIF1, CD74, TGFBI, SAMSN1, C1QA, FOS, C1QC
## MRC1, CD14, CTSC, SLC8A1, VAMP5, HLA-DRB1, MS4A6A, LYZ, TRPS1, PDGFC
## SLC9A9, SLCO2B1, FPR3, CLEC4A, FGL2, SELENOP, CLEC7A, C1QB, RTN1, CAMK1D
## PC_ 3
## Positive: MARCKS, LGMN, TPM4, CTSL, AIF1, FPR3, HLA-DQA1, CTSZ, HAMP, C1QC
## HLA-DQB1, CD36, CD163, CCL3, GYPC, IL18BP, ACTG1, BLVRB, STMN1, C1QA
## SMS, FMN1, MARCO, PLAU, FUCA1, FCGR3A, IFI30, COLEC12, MS4A7, MS4A6A
## Negative: PTGDS, LINC02244, CLU, NCAPH, LINC01010, CRIP1, AC015660.2, SYNGR1, AC067751.1, GCLC
## KCNMA1, RCBTB2, C1QTNF4, KCNJ1, DUSP2, SERTAD2, SPON2, NRCAM, C2orf92, TRIM54
## RGS20, CT69, FGL2, GPC4, RNF128, NFATC3, MX1, NCF1, ST5, IFIT3
## PC_ 4
## Positive: XIST, HLA-DRB5, GCHFR, SLC11A1, QPCT, MS4A7, PAX8-AS1, MSR1, SASH1, GPX3
## C22orf34, MARCO, AL035446.2, AC020656.1, ST6GAL1, FRMD4A, CCDC26, AL136317.2, S100A9, MIR646HG
## DLEU1, TMTC1, GPRIN3, CLMP, NMB, RARRES1, SERINC2, AC008591.1, OSBP2, SLC6A16
## Negative: PCLAF, TYMS, CLSPN, TK1, CENPM, RRM2, MYBL2, DIAPH3, MKI67, SHCBP1
## CDK1, BIRC5, HELLS, CEP55, ESCO2, FAM111B, NCAPG, NUSAP1, CENPK, ATAD2
## CENPU, TOP2A, MCM7, HMMR, PRC1, DTL, KIF11, TPX2, CENPF, GINS2
## PC_ 5
## Positive: XIST, HLA-DRB5, MKI67, TYMS, TK1, PCLAF, GCHFR, AC020656.1, RAD51AP1, CLSPN
## TMTC1, RRM2, CENPK, CDK1, QPCT, MYBL2, SHCBP1, BIRC5, DIAPH3, AL136317.2
## NCAPG, ESCO2, TDRD9, CENPM, ATAD2, TOP2A, CCNA2, DTL, CEP55, CENPF
## Negative: TMEM176A, PLEK, IL1RN, MIF, LINC01091, TUBA1A, MYL9, CYTOR, LINC00278, MMP19
## TTTY14, MGLL, HIV-Gagp17, PHLDA1, ADAMDEC1, HIV-BaLEnv, HIV-LTRU5, HDAC9, ACTG1, MMP7
## SERPINE1, TEX14, PPM1H, HIV-TatEx1, AC006001.2, HIV-Polprot, HIV-LTRR, RABGEF1, HLA-DPA1, UTYDimHeatmap(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 SNNmono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)## 11:54:30 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:54:30 Read 11741 rows and found 4 numeric columns
## 11:54:30 Using Annoy for neighbor search, n_neighbors = 30
## 11:54:30 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:54:31 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd7f1c49
## 11:54:31 Searching Annoy index using 1 thread, search_k = 3000
## 11:54:35 Annoy recall = 100%
## 11:54:36 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:54:38 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:54:38 Commencing optimization for 200 epochs, with 392702 positive edges
## 11:54:42 Optimization finishedDimPlot(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)
Extract monocytes MDM only
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 countsmono <- RunPCA(mono, features = VariableFeatures(object = mono))## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: PCAT5, ADGRF1, GSTP1, LINC01010,
## RPS6KA2, ANKRD44, CAMK1, CEP85L, RHOQ, ARHGAP10, DUSP1, ENTPD1, STPG2, KCNA2,
## SIGLEC15, TMCC1, TREM1.## PC_ 1
## Positive: ARL15, DOCK3, NEAT1, EXOC4, FTX, MALAT1, LRMDA, RASAL2, DPYD, JMJD1C
## TMEM117, FHIT, PLXDC2, VPS13B, MITF, ZEB2, DENND4C, MAML2, ATG7, ZNF438
## TPRG1, COP1, FRMD4B, ELMO1, ATXN1, JARID2, FMNL2, TCF12, ERC1, ARID1B
## Negative: MT-ATP8, MTRNR2L1, FABP3, FABP4, GAPDH, HAMP, S100A9, C1orf56, TUBA1B, PVALB
## AL136097.2, CARD16, FAH, BIRC7, C1QC, UTS2, MT2A, MT1A, LILRA1, BLVRB
## NUPR1, BEX3, GCHFR, AC025154.2, TXN, MT1G, PRDX6, FXYD2, CCL23, MYL9
## PC_ 2
## Positive: CYSTM1, FAH, GDF15, PSAT1, CD164, FDX1, ATP6V1D, BCAT1, CCPG1, SAT1
## SLAMF9, HES2, HEBP2, TXN, PHGDH, CSTA, NUPR1, GARS, TCEA1, STMN1
## RETREG1, CTSL, SCD, MARCO, S100A10, CLGN, RILPL2, SNHG5, SNHG32, NMB
## Negative: SPRED1, HLA-DRA, RCBTB2, KCNMA1, HLA-DRB1, CD74, HLA-DPA1, GCLC, FGL2, MRC1
## HLA-DPB1, SLCO2B1, C1QA, AC020656.1, RTN1, LINC02345, PDGFC, VAMP5, TGFBI, MERTK
## HLA-DRB5, STAC, CCDC102B, DPYD, ATP8B4, MX1, RNF128, CRIP1, NFATC3, AIF1
## PC_ 3
## Positive: NCAPH, CRABP2, TM4SF19, GAL, DUSP2, RGCC, CHI3L1, CCL22, TRIM54, TMEM114
## ACAT2, RGS20, TCTEX1D2, LINC02244, NUMB, CCND1, IL1RN, MGST1, SERTAD2, AC092353.2
## POLE4, CLU, ZNF366, GPC4, SYNGR1, SLC20A1, GCLC, OCSTAMP, AC067751.1, CRIP1
## Negative: MARCKS, CD14, BTG1, TGFBI, MS4A6A, HLA-DQA1, FPR3, CTSC, CD163, MEF2C
## MPEG1, TIMP1, AIF1, HLA-DPB1, C1QC, HLA-DQB1, NUPR1, SELENOP, NAMPT, RNASE1
## HLA-DQA2, ALDH2, TCF4, EPB41L3, ARL4C, ZNF331, HIF1A, SAMSN1, MS4A7, CLEC4A
## PC_ 4
## Positive: MT-ATP8, MTRNR2L1, XIST, PDE4D, AL035446.2, HIV-Polp15p31, AL365295.1, AC073359.2, LY86-AS1, CCDC26
## RAD51B, EPHB1, PRKCE, CLVS1, LINC02320, C5orf17, FTX, LINC01473, AC008591.1, FHIT
## DMGDH, HIV-LTRU5, MIR646HG, BCAS3, ZBTB20, AF117829.1, AC079142.1, LINC01135, LINC02649, PKD1L1
## Negative: CTSB, HLA-DRB1, TUBA1B, HSP90B1, CD74, GAPDH, ACTG1, HLA-DPA1, IFI30, HLA-DRA
## RGCC, HSPA5, C15orf48, AIF1, TUBB, LGMN, C1QA, CYP1B1, HLA-DPB1, FBP1
## CA2, TUBA1A, CCL3, TXN, HLA-DQB1, S100A4, MMP9, LYZ, RNASE6, CTSL
## PC_ 5
## Positive: TYMS, PCLAF, TK1, MKI67, MYBL2, CENPM, RRM2, BIRC5, CLSPN, SHCBP1
## CEP55, DIAPH3, CDK1, CENPK, PRC1, CENPF, ESCO2, NUSAP1, TOP2A, HELLS
## NCAPG, ANLN, FAM111B, TPX2, CCNA2, KIF11, ASPM, CENPU, MAD2L1, HMMR
## Negative: RGCC, MMP9, IFI30, CTSB, PTGDS, LINC02244, MGST1, CHI3L1, HLA-DRB1, LYZ
## C15orf48, HSPA5, CD74, FABP3, GCLC, GCHFR, CRABP2, TXN, HLA-DPA1, S100A10
## NCAPH, HLA-DRA, ISG15, TMEM176B, CDKN2A, HSP90B1, RALA, PLAUR, CA2, SAT1DimHeatmap(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 SNNmono <- FindClusters(mono, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono <- RunUMAP(mono, dims = 1:4)## 11:55:34 UMAP embedding parameters a = 0.9922 b = 1.112## 11:55:34 Read 28280 rows and found 4 numeric columns## 11:55:34 Using Annoy for neighbor search, n_neighbors = 30## 11:55:34 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:55:37 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd59b606f
## 11:55:37 Searching Annoy index using 1 thread, search_k = 3000
## 11:55:46 Annoy recall = 100%
## 11:55:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:55:50 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:55:51 Commencing optimization for 200 epochs, with 933992 positive edges
## 11:55:58 Optimization finishedDimPlot(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)
Cell counts
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 | 612 | 0 | 0 | 0 |
| alv_active3 | 1 | 30 | 682 | 0 | 0 | 0 |
| alv_bystander1 | 0 | 31 | 1105 | 2 | 0 | 0 |
| alv_bystander2 | 0 | 29 | 1304 | 0 | 1 | 1 |
| alv_bystander3 | 0 | 18 | 692 | 1 | 0 | 0 |
| alv_bystander4 | 0 | 31 | 729 | 3 | 0 | 0 |
| alv_latent1 | 0 | 9 | 664 | 0 | 0 | 0 |
| alv_latent2 | 0 | 22 | 1700 | 1 | 0 | 0 |
| alv_latent3 | 1 | 18 | 1605 | 0 | 0 | 0 |
| alv_latent4 | 0 | 42 | 1650 | 0 | 0 | 0 |
| alv_mock1 | 0 | 17 | 957 | 0 | 0 | 0 |
| alv_mock2 | 0 | 20 | 784 | 0 | 0 | 0 |
| alv_mock3 | 0 | 47 | 1305 | 1 | 1 | 0 |
| mdm_active1 | 0 | 17 | 863 | 0 | 0 | 0 |
| mdm_active2 | 0 | 3 | 476 | 0 | 0 | 0 |
| mdm_active3 | 0 | 13 | 407 | 0 | 0 | 0 |
| mdm_active4 | 0 | 4 | 432 | 0 | 0 | 0 |
| mdm_bystander1 | 1 | 44 | 1423 | 0 | 0 | 0 |
| mdm_bystander2 | 0 | 31 | 1274 | 0 | 0 | 0 |
| mdm_bystander3 | 0 | 44 | 472 | 2 | 0 | 0 |
| mdm_latent1 | 0 | 20 | 1172 | 0 | 0 | 0 |
| mdm_latent2 | 0 | 11 | 1268 | 0 | 0 | 0 |
| mdm_latent3 | 0 | 14 | 353 | 1 | 0 | 0 |
| mdm_mock1 | 0 | 5 | 791 | 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 | 78 | 3016 | 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.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0615764 | 0.00000 | 0.00000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0681199 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
| Dendritic cells | 3.5978836 | 1.290323 | 4.2075736 | 2.7240773 | 2.1722846 | 2.531646 | 4.0629096 | 1.337296 | 1.2768427 | 1.1083744 | 2.48227 | 1.74538 | 2.487562 | 3.4711965 | 1.931818 | 0.6263048 | 3.095238 | 0.9174312 | 2.9972752 | 2.375479 | 8.4942085 | 1.677852 | 0.8600469 | 3.8043478 | 0.6281407 | 1.062215 | 2.259887 | 0.4907975 | 2.5177534 |
| Monocytes | 96.0846561 | 98.709677 | 95.6521739 | 97.1001757 | 97.6779026 | 97.327708 | 95.5439056 | 98.662704 | 98.6651190 | 98.8300493 | 97.51773 | 98.25462 | 97.512438 | 96.3810931 | 98.068182 | 99.3736952 | 96.904762 | 99.0825688 | 96.9346049 | 97.624521 | 91.1196911 | 98.322148 | 99.1399531 | 95.9239130 | 99.3718593 | 98.937785 | 97.740113 | 99.5092025 | 97.3531311 |
| Neutrophils | 0.2116402 | 0.000000 | 0.0000000 | 0.1757469 | 0.0000000 | 0.140647 | 0.3931848 | 0.000000 | 0.0580383 | 0.0000000 | 0.00000 | 0.00000 | 0.000000 | 0.0738552 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.3861004 | 0.000000 | 0.0000000 | 0.2717391 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
| NK cells | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0749064 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.00000 | 0.00000 | 0.000000 | 0.0738552 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
| Progenitors | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0749064 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.00000 | 0.00000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 | 0.000000 | 0.0000000 | 0.0322789 |
Focus on MOCK vs ACTIVE - exclude latent and bystander
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 countsmono_focus_mdm <- RunPCA(mono_focus_mdm, features = VariableFeatures(object = mono_focus_mdm))## Warning: The following 43 features requested have zero variance; running
## reduction without them: NMRK2, KALRN, LINC00278, IGFN1, PPP2R2B, HPGDS, CKAP2L,
## FRMD4B, ADAMTS6, SH3RF3, MIR181A1HG, SLC51A, CD79A, LINC02739, CPEB2,
## STON1-GTF2A1L, NCAPD2, AL592156.2, NHSL1, LINC00342, AL139220.2, TPST1,
## AP000442.1, CTHRC1, GDF3, AK8, ERVMER61-1, HLA-C, ZNF208, GTSF1, GPR143, DUSP6,
## AC117473.1, LINC00382, GAK, RAPGEF1, LINC00563, CLYBL-AS2, AL133297.1,
## AC106820.2, CARMIL3, RALA, RBM46## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: HIST1H2BJ, OASL, CENPA, ENG, HIST1H3J,
## SSBP3, TVP23A, PCBP3, LINC00885, NCALD, SCAI, IFT80, AL355388.1, C5AR2, PIWIL2,
## AC013287.1, AL009177.1, AC106793.1, AL358944.1, KLF12, CCND2, AC108134.2,
## SLC2A14, TTK, FANCM, LINC01320, STXBP1, PPP6R2, TNRC6C, ZFPM2-AS1, HCAR2, RBM6,
## MT-ND1, AC064805.2, PNPLA7, KLRG1, AL353759.1, GPR155, TENM1, PARD3, SETD2,
## RASIP1, ZBTB41, EIF2B3, AGPAT5, AC015660.2, AC025262.1, NUP93, CLPX, LINC02853,
## FAM118B, ZZEF1, ESR2, CWC22, ATP6V0A2, MARCH1, SH3GL1, GYG2, AC034213.1,
## LINC02585, P2RY12, RPS4Y1, PDIA4, LRRK2, ANKRD34B, AC078923.1, PCA3,
## AL360015.1, AL157911.1, CLIC5, SF3B3, PKIA-AS1, AC007364.1, ASTL, NR2F1-AS1,
## MT-ND5, DOK5, AC005740.5, TMEM220-AS1, CCNB2, AURKC, AC007091.1, IFT140,
## AC133550.2, DUSP16, RAP1GDS1, AL355388.3, AL357873.1, AC072022.2, CHRNB4,
## ZNF276, BTBD2, AP000462.1, AC093766.1, FAM184A, PDGFD, SRRM2-AS1, AC108879.1,
## RORA, MAP4K1, MAP3K15, CPLX1, LINC01054, CREM, LPL, SIGLEC10, WWP1, AL021917.1,
## MT-ND2, LRRC9, PSME2, AL031658.1, JAKMIP3, TTC21B-AS1, STAT5B, VGLL3,
## AC099552.1, ACMSD, AC006427.2, ZNF385D-AS1, CFAP52, SEMA6D, TENM2, ANK3, MRC2,
## LDHA, GNG4, AP000446.1, CU638689.5, AC100849.2, LINC01010, RGS18, DAPK2,
## ACTR3C, DPEP3, SETBP1, CCDC7, AL592295.3, OBI1-AS1, AC079062.1, APCDD1L, NLRP4,
## GRM7, PACS1, FHOD3, MMP28, AC092811.1, AC105001.1, AC009315.1, AC240274.1,
## P2RY13, CDK19, AC089983.1, DIAPH2, TMEM212-AS1, FDXACB1, ARC, AC092490.1,
## SAMD11, GAS1RR, CLEC4C, AC114977.1, RAB10, CHRM3, SF3B1, MANF, MIR155HG,
## RASSF4, SRGN, CORIN, AHR, KIF21A, SNAP25-AS1, AL132775.1, IFI6, SLC44A5,
## CARD11, TUBB3, AC019330.1, AL391807.1, TEX15, OR8D1, PRRT2, AC005393.1, VIPR1,
## AC007389.1, CRISPLD2, ADRA2B, UBASH3B, AC016074.2, TPTE2, AC026333.3,
## AP002793.1, AC084809.2, AC073050.1, TRIM71, LINC01600, LCTL, ZDHHC17, CD93,
## AC073569.1, LRRIQ4, AL049828.1, AC011365.1, Z99758.1, NWD1, CEP112, EPHA6,
## AC019117.1, PLEK, CEP135, STPG2, TMEM45A, TFG, AC137770.1, ENPEP, AC009264.1,
## AC096570.1, SESN3, DHCR24, U62317.4, DHX8, LINC00649, OR3A2, BCL2A1, ABCA13,
## LINC00237, TRIM37, DNAH2, RPS6KA6, LNX1, SRGAP3, GAL3ST4, GFRA2, TACC3, BCAR3,
## CPE, SNHG12, CREG2, LYPLAL1-DT, AC011509.2, SULF2, IGSF6, LINC01917, LINC01855,
## MLLT3, AL589740.1, HIST1H2AC, HIST1H2BG, RTKN2, PRTG, OSBPL6, SIGLEC7,
## AL445584.2, MPDZ, LINC02208, AL136418.1, MCF2L-AS1, DPF1, PAX8-AS1, AC005264.1,
## COX5B, Z94721.1, FKBP1B, LINC00958, MAPK8, RUFY2, DARS, AC021231.1, SOX15,
## PODN, LINC01353, NYX, MT-CO3, AC124017.1, MMD2, LINC02449, LINC01891, RRAS2,
## TASP1, AC104459.1, AC087683.3, HIST2H2BE, SCAPER, AC090630.1, CNDP1, HCAR3,
## AKR7A2, RAP2C-AS1, CKMT1A, MRVI1, AC015849.1, AC011476.3, TMCC1, AC021086.1,
## CACNA2D3, RFX7, TMEM71, MNAT1, AC073475.1, MCCC2, AL031710.2, AL031005.2,
## ZBTB46, CD72, TUBB4B, JAKMIP2-AS1, AP000829.1, MAP1A, ATF3, CFAP69, AC024901.1,
## AC129803.1, KCNJ1, AC010834.2, ZHX2, KIAA0825, KDM6A, PPP1R12B, AC008115.2,
## AL356379.2, LINC02112, LINC01191, GSN, MNDA, RALGAPA2, SPTLC3, LINC02015,
## COL8A2, C11orf45, GLIPR1, PTPRG-AS1, ITGA4, ERO1A, RHCE, RHOF, LGALSL,
## AC098588.3, P3H2, PDE11A, TUBB4A, TNS3, NPRL3, SCFD2, MT-ATP6, AF131216.1,
## AC103726.2, AC008555.2, SPEG, SOCS3, MREG, SPIB, AL110292.1, ALDH1A1,
## AC097654.1, MKX, BMPR1B, SLC26A7, AC114763.1, AC020637.1, PAN3, PROCR, AEBP2,
## OXSR1, MARCH3, PTPRB, AC020743.2, PRSS3, UBE2T, LINC00571, RFX3, MCM9,
## AC246817.1, LILRB2, MAFB, ADORA2B, SPAG7, LINC02851, NETO2, GMDS, TMEM131L,
## TRIO, KIAA1841, TMEM72-AS1, MAPK13, CLDN4, E2F1, AL365295.2, BFSP2, GRK3,
## RNF24, BUB1, ITM2A, ANGPTL4, NEGR1, CHM, PTPN13, SGO1, RNF180, XKR9,
## AC079584.1, AC006525.1, LINC00987, ST3GAL6, TRAF2, CALR, ACTB, PIP4K2C, PLA2G5,
## AL137076.1, AL158839.1, AL589765.6, C4BPA, AC096639.1, AC114810.1, AL133243.1,
## AC009229.3, AC012358.1, LINC01104, AC063944.3, AC092902.6, AC128709.2,
## AC021220.2, ADAMTS3, AC116351.1, AC011374.3, LINC02533, LINC02571, Z84484.1,
## MRAP2, FUT9, AL357992.1, AL078582.1, AC002480.2, TRIM74, DEFB136, PRDM14,
## CALB1, AF178030.1, AF235103.1, AL353764.1, TMEM246-AS1, PPP3R2, AL356481.2,
## AL731559.1, AL121748.1, LINC02625, SYCE1, OR10A4, LINC02751, SAA4, AC013714.1,
## AC024940.1, LINC01479, AC012464.3, C1QTNF9-AS1, SMAD9-IT1, AL161717.1,
## AL442125.2, KLHL33, CMA1, AL163953.1, AC048382.1, AC091167.5, BAIAP3, NPIPB8,
## AC012186.2, AC092378.1, AC129507.2, RAI1-AS1, AC007952.6, AC004448.3,
## AC243773.2, KRT19, AC090409.1, OR7E24, ANGPTL8, LINC01858, AC022150.2,
## AL021396.1, LINC01747, CU634019.5, NUDT11, LINC00266-4P, CFD, HCRTR2, DNAH3,
## LINC00407, AL136298.1, FILIP1L, LINC02398, AL096794.1, EXD2, KIF20B, ZC4H2,
## GLRX, HERPUD1, ATP13A3, AL645929.2, AC083836.1, PRKG2, ADIRF, STUM, STAG1,
## CDT1, SUCLG2-AS1, ARRDC3-AS1, SPOCK3, POF1B, KCNA2, NFKB1, BACH1, TFRC,
## LINC01762, ATP6V0D2, CDKAL1, LINC00511, NIPAL2, CPXM2, PCLO, ENTPD1, CD200R1,
## AC092134.3, ITGA9, ZFP36L1, MGAT5, FAM110B, TIMP4, ADAMTS10, UFL1-AS1, DUSP5,
## FAM135A, AC021546.1, MYADM, PRDM1, PFKFB3, LINC02798, LINC01572, COBL, LY96,
## RHPN2, TNIK, SPATA5, AC092844.1, CACNA2D2, AHCYL2, EDA, SH3D19, MCTP2, GADD45G,
## CASC2, BDNF-AS, PRH1, MDH1, GNG7, DRD5, FAM189A2, AC007391.2, TSPAN8,
## AC008438.2, AC079804.3, PRORP, ST3GAL3, AL158068.2, CASP1, TTLL7, DOCK10,
## ELANE, CRELD2, FANCC, PLBD1, KL, AC007262.2, PMAIP1, MOCOS, SLC41A2, LINC02457,
## PCNX2, UHRF1, AC097515.1, RALYL, MATN1, TRAF1, ITPR2, NCAPD3, NPTX2, MTUS1,
## FOXO6, LAMC1-AS1, AC099782.2, WWTR1-AS1, AC104785.1, SLC6A3, AL078599.2,
## AC010719.1, AC073210.3, AP005436.3, AP002884.1, LUM, AC073863.1, AC063943.1,
## LINC02280, AC011939.3, AC092127.1, AC104423.1, L1CAM, LGALS3, DYM, PAFAH1B1,
## AL157834.1.## PC_ 1
## Positive: FMN1, NEAT1, FHIT, RAD51B, ARL15, MALAT1, SNTB1, AL035446.2, MBD5, FTX
## FNDC3B, EXOC4, VTI1A, PDE4D, GMDS-DT, COP1, TBC1D8, JMJD1C, SLC22A15, SLC8A1
## DENND1A, CBLB, DOCK3, DENND4C, VPS13B, DOCK4, LRMDA, SBF2, SMYD3, COG5
## Negative: SLC35F1, TUBA1B, PRDX6, CRABP2, TUBA1A, GAL, MT-ATP8, PTGDS, DUSP2, C15orf48
## FABP3, LILRA1, MT2A, CRIP1, NCAPH, MYL9, S100A4, MMP7, LINC02244, CCNA1
## CA2, UCHL1, MT1E, TMEM114, RGCC, SEL1L2, ISG15, RNF128, ACAT2, SPON2
## PC_ 2
## Positive: HIV-Gagp17, HIV-Polp15p31, HIV-BaLEnv, HIV-Polprot, HIV-LTRU5, HIV-Vif, GDF15, PSAT1, HIV-Gagp2p7, BEX2
## PHGDH, SNCA, TCEA1, SLAMF9, OCIAD2, B4GALT5, G0S2, HIV-Gagp1Pol, HIV-TatEx2Rev, HIV-EGFP
## RAB38, SUPV3L1, S100A10, HIV-Vpu, UGCG, TXN, FDX1, TM4SF19, CLGN, HIV-TatEx1
## Negative: HLA-DRB1, HLA-DRA, TGFBI, CD74, CEBPD, HLA-DPA1, PTGDS, HLA-DPB1, SIPA1L1, LYZ
## KCNMA1, RCBTB2, MS4A6A, EPB41L3, RNASE6, HLA-DRB5, SSBP2, ATP8B4, AOAH, RGS2
## HLA-DQB1, SAMSN1, FOS, ST8SIA4, PDE4B, C1QA, DOCK4, TRPS1, GCLC, DPYD
## PC_ 3
## Positive: BTG1, G0S2, MARCKS, CD14, ARL4C, TGFBI, HIV-Gagp17, HIV-Polp15p31, HIV-Polprot, CLEC4E
## SEMA4A, CXCR4, MS4A6A, VMO1, RNASE1, P2RY8, MEF2C, HIV-BaLEnv, HIV-Vif, TIMP1
## TCF4, HIV-LTRU5, FUCA1, HIF1A, MPEG1, CEBPD, FPR3, HIV-Gagp2p7, MS4A7, PDK4
## Negative: NCAPH, CRABP2, RGCC, CHI3L1, LINC02244, TMEM114, TM4SF19, GPC4, RGS20, ACAT2
## NUMB, PSD3, LRCH1, GCLC, TCTEX1D2, AC005280.2, SLC28A3, AC092353.2, GAL, FNIP2
## DUSP2, FBP1, MICAL3, DOCK3, CCND1, S100A4, ASAP1, CSF1, RAI14, SLC25A48
## PC_ 4
## Positive: AC078850.1, ACTG1, TUBA1B, TUBB, FABP4, LGMN, CTSL, TPM4, SLC35F1, HSP90B1
## INSIG1, CD36, TUBA1C, FBP1, IL18BP, TUBA1A, PHLDA1, MGLL, CD300LB, CCL7
## CADM1, C3, PRDX6, CSF1, SDS, MMP9, FUCA1, ANXA1, TIMP1, TIMP3
## Negative: MT-ATP8, HIV-Polp15p31, HIV-Vif, HIV-Polprot, PDE4D, HIV-BaLEnv, HIV-Gagp2p7, HIV-LTRU5, HIV-Gagp17, HIV-Gagp1Pol
## BEX2, HIV-Vpu, HES2, CLGN, NIBAN1, ST6GALNAC3, MT1G, HIV-EGFP, S100P, HIV-TatEx1
## EPHB1, XIST, HES4, AL365295.1, C5orf17, CPD, LINC02320, HIV-TatEx2Rev, CLVS1, HIV-Vpr
## PC_ 5
## Positive: HIV-Gagp17, HIV-Vif, HIV-Polprot, HIV-BaLEnv, HIV-Polp15p31, HIV-Gagp2p7, HIV-Gagp1Pol, HIV-LTRU5, HIV-TatEx1, HIV-LTRR
## HIV-TatEx2Rev, HIV-Nef, HIV-Vpu, HIV-EGFP, HIV-EnvStart, HIV-Vpr, LINC01091, IL1RN, SLC35F1, TIMP3
## CCL7, TNFRSF9, INSIG1, C3, PCYT1A, MADD, AC005062.1, CCL3, TPM2, MACC1
## Negative: CCPG1, STMN1, RAB6B, NUPR1, HES2, PLEKHA5, PSAT1, CLEC4E, LDHB, PCDH19
## CP, CSRP2, CARD16, PLAUR, AC073359.1, SUPV3L1, SCG5, FOLR2, SOD2, NMB
## S100P, FABP3, GSTM3, MARCO, PHGDH, TCEA1, SLAMF9, BX664727.3, RETREG1, CXCR4DimHeatmap(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 SNNmono_focus_mdm <- FindClusters(mono_focus_mdm, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono_focus_mdm <- RunUMAP(mono_focus_mdm, dims = 1:4)## 11:56:10 UMAP embedding parameters a = 0.9922 b = 1.112## 11:56:10 Read 4605 rows and found 4 numeric columns## 11:56:10 Using Annoy for neighbor search, n_neighbors = 30## 11:56:10 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:56:11 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd3b04aa93
## 11:56:11 Searching Annoy index using 1 thread, search_k = 3000
## 11:56:12 Annoy recall = 100%
## 11:56:13 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:56:15 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:56:15 Commencing optimization for 500 epochs, with 162620 positive edges
## 11:56:18 Optimization finishedDimPlot(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 countsmono_focus_alv <- RunPCA(mono_focus_alv, features = VariableFeatures(object = mono_focus_alv))## Warning: The following 50 features requested have zero variance; running
## reduction without them: PRG2, AC008574.1, CD163, DNAJC5B, HIV-Gagp1Pol, CCL22,
## DIRC3, SLIT3, GPR183, LINC01500, SLC28A3, LGI2, HIV-Vpu, OLR1, TSPAN18, NCAPH,
## ARHGAP24, AL096794.1, STAC, LNCOG, SLC7A14-AS1, HIV-TatRev, GIGYF2, PARD3B,
## IL36B, AC084809.2, RGS8, AC093895.1, IGSF23, ME3, AC073569.2, AC092813.1,
## ULBP2, UST, STXBP5-AS1, NES, AC113137.1, DUSP2, LINC01198, EXOSC10, ACSBG2,
## TMEM233, CARD16, ZNF367, ITSN1, GABRR3, LINC02577, AEBP2, LINC01739, VTN## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: STEAP1B, LINC00964, NUP210L, CPNE8,
## AC083837.1, HCAR2, CD28, AC011396.2, AP002075.1, SFMBT1, PPP1R1C, EDIL3,
## ACVR2A, IQCA1, AC022031.2, SLC6A16, SPINK6, EFCAB7, LINC02068, TENM4, SRL, IL6,
## EXO1, C11orf49, IFI6, COL25A1, SMS, DNAH12, HIST1H2BG, DNAJC9, LINC02269,
## NPTX2, AC026391.1, GNLY, AC076968.2, KIAA1755, AL356010.2, EMP3, MS4A5, ACTN2,
## INO80, AC207130.1, BTNL8, LINC00589, AC135050.3, CNTNAP2, AC013391.2,
## AC092821.3, ARMC8, ADAM22, SMCHD1, AC084740.1, PLAAT3, TEK, LINC02073, NWD1,
## GCH1, TRIM2, MSRB3, GPRC5D, DSG2, MIF, AC005740.5, BANK1, ERAP2, CD226, SCFD1,
## AHCTF1, COL27A1, MARCKSL1, CTSW, CH25H, BRMS1L, AL353596.1, PMAIP1, KIF6,
## TUBB4A, MNDA, SHCBP1L, LINC00350, LINC02789, RAD51C, AL161646.1, RFX3-AS1,
## USP10, ZC3H8, NEURL3, PTPN2, ASAP3, MRPS6, AC104596.1, LINC01344, DYNC1H1, XDH,
## CCDC85A, AC010834.2, MIR4300HG, RFTN2, YJEFN3, LINC00299, MECP2, CR1, TNFAIP3,
## LPCAT1, LINC00511, AC005264.1, AC004231.1, SPOCK1, LINC01862, SH3TC2,
## AC079313.2, RPH3AL, NEK4, SHANK2, URB1, STK32B, SVEP1, CFAP74, SPACA3, XPO5,
## PPEF1, ZPLD1, AC007785.1, AC019131.1, HSF2BP, RLF, FBXO4, ANGPT1, NCMAP,
## MBOAT4, ABCB1, SMAD1-AS2, AC125421.1, LINC02466, PTX3, NIPAL2, LINC00842, TEPP,
## EML2-AS1, HNRNPM, SYT12, NDRG1, WDR90, SOX6, POU6F2, TIMELESS, LGALS1, ZSCAN32,
## GRID2, CCNC, NRG1, C2orf72, RXFP1, FCRLB, ADCY3, PMPCA, GRHL2, STS, FSD1,
## INKA2-AS1, LINC01414, FO393415.3, UPK1A-AS1, GTF2IRD2, ERI2, DDN-AS1, DCSTAMP,
## PLEK, CCDC57, PAX8-AS1, AC124852.1, HCAR3, RFX8, GNRHR, MID2, CLDN18, PTPRB,
## GATA2, CFD, ABRAXAS2, AKAP6, AC009435.1, AC041005.1, AC009292.2, HRH2,
## AC107081.2, WIPF3, AC011346.1, NNMT, PDE1C, MIR2052HG, MFSD11, ANOS1, IL17RB,
## LMNB1-DT, DACH1, PHACTR1, CSTB, TMEM37, RELN, HLTF-AS1, LINC00609, AC006441.4,
## AC022035.1, AL662884.3, C22orf34, AL513166.1, CNR2, LINC01299, NCR3, IGF2R,
## COBL, EAF2, AOC1, CKAP5, TNR, LINGO2, ZNF157, LINC01999, AC007001.1, CD5L,
## AC007100.1, SFTPD-AS1, NBAT1, CTSZ, PLCH1, AL592295.3, GOLGA7B, RGS7, SLC2A5,
## AC068228.3, SPRY4-AS1, SLC12A3, AL353595.1, DZIP1, AC025430.1, TPSAB1, SLC44A5,
## UNC5D, SLC9A2, SVOP, ENTPD1, RAB7B, PSME2, HIP1, ABLIM1, ADM, KRTAP10-4, NEIL3,
## CDK6, NDRG2, STPG2, HSPA1B, KIAA0825, NIPAL4, CTXND1, CDHR3, LINC00378, DHRS9,
## Z99758.1, KDR, LDHA, AC109454.3, AL354733.2, NSG2, AC004949.1, CLSTN2, PLXDC1,
## MRC2, LINC00461, IKZF3, MAS1, RHOF, SLC6A7, GALNT14, AL645933.3, VIPR1, CSMD2,
## TNC, BPIFC, AC093898.1, SHROOM4, LINC01358, FOXM1, MTFMT, OASL, SLC4A1,
## AC110491.2, KIF21A, MIR3142HG, AC006460.1, NRIP3, AC116563.1, AL049828.1,
## SH3TC2-DT, AC092957.1, AL360015.1, CHODL, AC010997.3, CRIM1, GAPLINC, GPRIN3,
## ITGA2, NRCAM, LINC01320, CHDH, PRKCG, ZNF385D-AS1, AC006525.1, SERINC2, SNTG1,
## GINS2, LINC02109, AC005753.2, LGALS3, SCIN, AL109930.1, LINC01900, APOM,
## ZSCAN5A, STXBP5L, CYP4F22, POU2F2, MSH4, C1orf143, MEP1A, AC084816.1, CHD5,
## FAXC, ATP1B2, LINC00571, ARHGAP6, TROAP, CNGB1, GLIPR1, PARD3, CLPB,
## AC108066.1, SLC23A1, SLA, AIM2, AC011893.1, TMEM45A, PCLO, ABCG2, GPX3, YLPM1,
## AC002454.1, AC117453.1, CHAC1, CLEC7A, CRIPT, RBL1, DENND2A, BRCA1, FHAD1,
## DPYS, AMPD1, AC093307.1, TSGA13, ARSF, HMOX1, MCAM, AC092718.1, AC093010.2,
## TEX49, COL4A5, LINC00973, CD1E, AL359220.1, RNF212, AC016587.1, PRRX2, ZNF365,
## MOSPD1, SIK2, TNFRSF11A, AC007529.2, SEMA6D, CFAP57, PPP1R16B, LIPG, TMEM132C,
## SLC23A2, TGFB3, AC099489.1, MREG, DHCR24, RHOD, ELOC, ARL9, FAF1, PCDH15,
## AC079742.1, MEI4, ID2, KCNA2, AC016831.7, SRGAP3, MEIS1, SAMD12, UBE2E2, TDRD9,
## LINC02698, EFNA2, TBC1D24, ING3, AL713998.1, GNGT1, LINC00519, LIN28B-AS1,
## AC007381.1, STMND1, AC015908.2, RBPJL, AC246817.1, LCP1, SLC22A2, LINC02752,
## XKR9, AL355981.1, CCDC141, MOBP, AC097487.1, TMEM213, LINC02621, AC087897.2,
## AL160035.1, TCEA3, AC079298.3, AC137810.1, AC090515.6, AC120498.4, AC024901.1,
## INPP4B, AHCYL2, SLC4A8, MX2, SPSB1, DNAH2, PARN, ANKH, SNAI3, SLC35F3, CEP126,
## AC104041.1, PLTP, AL591845.1, SULT6B1, KLB, STXBP6, LINC01924, AC008443.9,
## RASL10A, DEGS2, GALNTL6, WDFY4, GNAI1, TNIK, AL136456.1, MAP1A, AC099560.1,
## COL8A2, LINC01572, FRMD6-AS2, LINC02805, RALGAPA2, AC010343.3, ST3GAL6, LGR4,
## GLUL, LINC01800, AP001496.1, CPEB2-DT, ELOVL5, CAMK1, SOX5, MB21D2, AC073091.4,
## PTP4A2, FRRS1, KDM7A, PAX5, NANOS1, CFI, LINC01948, ROBO4, IAPP, WASF3-AS1,
## AC130456.2, SEPTIN4-AS1, ELF5, FUT2, FBXW7, LINC02742, SPOPL, TNFRSF12A,
## FBXO43, AC079163.2, ASMT, MARCH1, AP000424.1, BX004807.1, EOGT, PLAT, CNTNAP5,
## HPN, STUM, CLSTN3, SKA3, GALNT18, LINC02777, NR6A1, GLT1D1, PDLIM4, BICD1,
## CNGA4, AC087280.2, HECW1, OPRD1, AC239860.4, LHCGR, AC103563.9, NECTIN3-AS1,
## AC010307.4, ZBED9, AC068305.2, SCG2, AC145146.1, HIST2H2BE, IL3RA, GLCCI1,
## AL645465.1, ASTN2, DNAJC1, WDR54, AC006333.1, SYT10, MYO16-AS1, AC096531.2,
## LINC01276, AL354811.2, OR10G3, NUDT10, LINC00894, SSPO, CIDEC, GRIK5, TREM1,
## CDT1, TMEM131L, C11orf45, FA2H, RBM11, PPP1R14C, HIST1H2AC, GRXCR2, SH3PXD2B,
## LPAR1, ZNF431, CORO1A, SERPINA1, RASSF4, GRAMD2A, JAML, LIMCH1, CDH12, TRMT5,
## AMPD3, ERBB4, AC055855.2, KCNJ1, S100A6, CDC5L, KCP, GPAT3, MPDZ, LMCD1-AS1,
## FCMR, GM2A, RBPJ, EMP1, TYW1B, SYNE1, DAO, PROSER2-AS1, AF274853.1, TPPP3,
## AHRR, ACTB, SCFD2, EML4, AP001636.3, LMO4, AC005280.1, FABP6, AC011476.3,
## NCAM2, ITSN2, IGSF6, HIST1H1C, SLC25A23, ANKS1B, NME5, SLC30A10, ATP13A3,
## AC009299.3, AC021851.2, LINC01933, BTBD11, AP000812.1, PKIB, FRMD6-AS1, BARD1,
## MIR155HG, AC104809.2, AC092801.1, SLC30A8, FAM92B, RNF217-AS1, MARCH3, ZNF543,
## HTRA4, AC007262.2, AC079793.1, AC068587.4, AXL, ATP1B1, BTG3-AS1, LINC01010,
## MED13L, MMP1, SVOPL, AL033530.1, ADGRL4, AL592402.1, DMRT1, AL136119.1,
## AC018618.1, FBLN1, AL357146.1, NEBL, LUM, LINC01169.## PC_ 1
## Positive: NUPR1, FABP4, STMN1, PSAT1, PHGDH, S100P, CLGN, G0S2, IFITM3, MARCKS
## ALDH2, DDIT4, BTG1, BLVRB, CAMP, ADAMDEC1, GDF15, TRIB3, BEX3, AC073359.1
## BEX2, CD14, FOLR2, RAB6B, CLEC4A, PCDH19, TIMP1, FAH, FXYD2, GYPC
## Negative: RASAL2, DOCK3, TMEM117, AC092353.2, DPYD, ASAP1, LINC01374, PLXDC2, LRMDA, CPEB2
## ATG7, NEAT1, FMNL2, LRCH1, MITF, MAML2, ARL15, ELMO1, TPRG1, ATXN1
## NUMB, EXOC4, DENND1B, MALAT1, PPARG, FHIT, VWA8, MYO1E, RABGAP1L, MEF2A
## PC_ 2
## Positive: SLC8A1, MRC1, RCBTB2, LYZ, AL356124.1, FCGR3A, TRPS1, SAMSN1, NRP1, ATP8B4
## CTSC, PDGFC, ME1, FCHSD2, ARHGAP15, XYLT1, RARRES1, HLA-DPA1, DOCK4, SELENOP
## CAMK1D, CCDC102B, SPRED1, SLCO2B1, HLA-DRA, KCNMA1, ZEB2, FLRT2, HLA-DPB1, TGFBI
## Negative: HIV-Nef, HIV-TatEx1, HIV-LTRU5, HIV-BaLEnv, HIV-Gagp17, HIV-Polprot, HIV-LTRR, HIV-Polp15p31, IL6R-AS1, GAL
## AL157912.1, HIV-TatEx2Rev, HIV-EnvStart, HIV-Vif, HIV-Vpr, HIV-Gagp2p7, SLC20A1, HIV-EGFP, IL1RN, DUSP13
## TRIM54, CYTOR, MYL9, MIR4435-2HG, CIR1, GPC3, NMRK2, POLE4, CSF1, BIN2
## PC_ 3
## Positive: CTSL, TXN, FDX1, CD164, SAT1, TPM4, BCAT1, DNASE2B, CYSTM1, BLVRB
## MARCO, ACTG1, H2AFZ, FAH, S100A9, GAPDH, ATP6V1D, SPP1, TXNRD1, SCD
## CTSB, UPP1, LGMN, TUBB, HAMP, CCL3, FABP4, CSF1, STMN1, ANXA1
## Negative: HIV-TatEx1, HIV-Nef, HIV-LTRR, HIV-Polprot, HIV-Vpr, HIV-EnvStart, HIV-Vif, HIV-LTRU5, HIV-Polp15p31, HIV-BaLEnv
## HIV-Gagp17, HIV-Gagp2p7, HIV-TatEx2Rev, HIV-EGFP, AC067751.1, LINC02408, AC092944.1, RTN1, MMP7, CRIP1
## DPH6-DT, RNF128, ST5, LINC02345, CLU, DPYD, SGO1-AS1, KCNMA1, AC008591.1, MSI2
## PC_ 4
## Positive: PTGDS, LINC02244, SYNGR1, TRIM54, RGS20, KCNMA1, SERTAD2, C2orf92, CLU, CRIP1
## MGST1, TMEM114, MT2A, MT1E, GCLC, ZNF366, CCND1, ACAT2, OCSTAMP, MMP7
## AC067751.1, SPON2, TMEM176B, KCNQ5, PRDX6, RGCC, RGS16, LINC02801, TNFSF14, GPC4
## Negative: HIV-TatEx1, HIV-Nef, HIV-Polprot, HIV-Polp15p31, HIV-BaLEnv, HIV-LTRR, HIV-Vif, HIV-Vpr, HIV-EnvStart, HIV-LTRU5
## HIV-Gagp17, HIV-Gagp2p7, HIV-TatEx2Rev, HIV-EGFP, AIF1, FPR3, ID3, MRC1, MARCKS, IL7R
## COLEC12, ALOX5, CCL2, C1QA, SNCA, CTSC, SLC8A1, FMN1, LGMN, TMEM236
## PC_ 5
## Positive: XIST, CCDC26, LINC02320, OSBP2, MIR646HG, AL136317.2, ZFPM2, PKD1L1, SKAP1, AC008591.1
## LINC01340, AC012668.3, AC068413.1, SEL1L2, LINC01374, LINC01135, BX664727.3, QPCT, LINC00607, ANO5
## PSD3, LIX1-AS1, PLEKHA5, LINC01708, CLMP, AP005262.2, STOX2, AL365295.1, PLXNA2, TMTC1
## Negative: IL1RN, TUBA1A, MMP7, TIMP3, TUBA1B, CTSB, ALOX5AP, SLC35F1, AIF1, TMEM176A
## HLA-DRB1, CD74, TMEM176B, GAPDH, ACTG1, HLA-DRA, LINC02345, ALCAM, GCLC, HLA-DPA1
## RNF128, MMP9, RGCC, TUBB, CRIP1, S100A4, LINC00910, VAMP5, SPRED1, AC007952.4DimHeatmap(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 SNNmono_focus_alv <- FindClusters(mono_focus_alv, algorithm = 3, resolution = 0.3, verbose = FALSE)
mono_focus_alv <- RunUMAP(mono_focus_alv, dims = 1:4)## 11:56:32 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:56:32 Read 5248 rows and found 4 numeric columns
## 11:56:32 Using Annoy for neighbor search, n_neighbors = 30
## 11:56:32 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:56:32 Writing NN index file to temp file /tmp/RtmpU7K1Jd/file1562fd35b3d3e8
## 11:56:32 Searching Annoy index using 1 thread, search_k = 3000
## 11:56:34 Annoy recall = 100%
## 11:56:35 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:56:37 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:56:37 Commencing optimization for 500 epochs, with 185406 positive edges
## 11:56:41 Optimization finishedDimPlot(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)
Differential expression
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 72761 7103
## mdm_mock1|AAACGCTCATCAGCGC mac 49143 6282
## cell sample RNA_snn_res.0.5
## <character> <character> <factor>
## mdm_mock1|AAACGAATCACATACG mdm_mock1|AAACGAATCA.. mdm_mock1 4
## mdm_mock1|AAACGCTCATCAGCGC mdm_mock1|AAACGCTCAT.. mdm_mock1 13
## seurat_clusters cell_barcode
## <factor> <character>
## mdm_mock1|AAACGAATCACATACG 4 mdm_mock1|AAACGAATCA..
## mdm_mock1|AAACGCTCATCAGCGC 13 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 4
## mdm_mock1|AAACGCTCATCAGCGC 13colnames(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)))## HIV-LTRR HIV-LTRU5 HIV-Gagp17 HIV-Gagp24 HIV-Gagp2p7
## alv_active1 0 20 1 0 0
## alv_active2 0 0 0 0 0
## alv_active3 0 2 0 0 0
## alv_bystander1 0 0 0 0 0
## alv_bystander2 0 0 0 0 0
## alv_bystander3 0 0 0 0 0
## alv_bystander4 0 0 0 0 0
## alv_latent1 0 0 0 0 0
## alv_latent2 0 0 0 0 0
## alv_latent3 0 9 0 0 1
## alv_latent4 0 0 0 0 0
## alv_mock1 0 0 0 0 0
## alv_mock2 0 0 0 0 0
## alv_mock3 0 0 0 0 0
## mdm_active1 0 0 0 0 0
## mdm_active2 0 0 0 0 0
## mdm_active3 0 0 0 0 0
## mdm_active4 0 0 0 0 0
## mdm_bystander1 0 0 0 0 0
## mdm_bystander2 0 0 0 0 0
## mdm_bystander3 0 0 0 0 0
## mdm_latent1 0 0 0 0 0
## mdm_latent2 0 0 0 0 0
## mdm_latent3 0 0 0 0 0
## mdm_mock1 0 0 0 0 0
## mdm_mock2 0 0 0 0 0
## mdm_mock3 0 0 0 0 0
## mdm_mock4 0 0 0 0 0
## react6 0 0 0 0 0
## HIV-Gagp1Pol
## alv_active1 0
## alv_active2 0
## alv_active3 0
## alv_bystander1 0
## alv_bystander2 0
## alv_bystander3 0
## alv_bystander4 0
## alv_latent1 0
## alv_latent2 0
## alv_latent3 0
## alv_latent4 0
## alv_mock1 0
## alv_mock2 0
## alv_mock3 0
## mdm_active1 0
## mdm_active2 0
## mdm_active3 0
## mdm_active4 0
## mdm_bystander1 0
## mdm_bystander2 0
## mdm_bystander3 0
## mdm_latent1 0
## mdm_latent2 0
## mdm_latent3 0
## mdm_mock1 0
## mdm_mock2 0
## mdm_mock3 0
## mdm_mock4 0
## react6 0plotMDS(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))DE Prep
pbmono <- assays(pb)[["Monocytes"]]
head(pbmono)## alv_active1 alv_active2 alv_active3 alv_bystander1 alv_bystander2
## HIV-LTRR 4532 3355 2708 0 0
## HIV-LTRU5 168788 129485 107460 0 0
## HIV-Gagp17 33137 27369 17430 38 45
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 1239 1259 798 4 4
## HIV-Gagp1Pol 2199 2384 1723 6 9
## alv_bystander3 alv_bystander4 alv_latent1 alv_latent2 alv_latent3
## HIV-LTRR 0 0 92 390 411
## HIV-LTRU5 0 0 2546 13645 13744
## HIV-Gagp17 37 49 571 2405 1947
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 5 7 18 87 129
## HIV-Gagp1Pol 18 13 31 160 220
## alv_latent4 alv_mock1 alv_mock2 alv_mock3 mdm_active1 mdm_active2
## HIV-LTRR 599 31 60 206 2379 1520
## HIV-LTRU5 22798 764 1397 7389 105566 73186
## HIV-Gagp17 2780 107 186 1541 38988 23811
## HIV-Gagp24 0 0 0 0 0 0
## HIV-Gagp2p7 146 2 8 53 1512 1041
## HIV-Gagp1Pol 271 6 23 97 2120 1415
## mdm_active3 mdm_active4 mdm_bystander1 mdm_bystander2
## HIV-LTRR 1571 1665 0 0
## HIV-LTRU5 98480 56843 1 0
## HIV-Gagp17 41474 15181 169 153
## HIV-Gagp24 0 0 0 0
## HIV-Gagp2p7 2441 420 18 12
## HIV-Gagp1Pol 3167 794 14 23
## mdm_bystander3 mdm_latent1 mdm_latent2 mdm_latent3 mdm_mock1
## HIV-LTRR 2 293 327 97 37
## HIV-LTRU5 7 8200 9557 3846 587
## HIV-Gagp17 44 2226 2276 676 135
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 4 79 131 64 8
## HIV-Gagp1Pol 8 129 190 106 8
## mdm_mock2 mdm_mock3 mdm_mock4 react6
## HIV-LTRR 41 9 39 262
## HIV-LTRU5 752 118 1125 7662
## HIV-Gagp17 271 41 161 602
## HIV-Gagp24 0 0 0 0
## HIV-Gagp2p7 8 1 5 23
## HIV-Gagp1Pol 14 3 13 42dim(pbmono)## [1] 36622 29hiv <- pbmono[1:2,]
pbmono <- pbmono[3:nrow(pbmono),]Gene sets
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")DE1 Latently- vs productively-infected cells (groups 3 vs 4).
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
## HIV-Gagp17 38988 23811 41474 15181 2226
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 1512 1041 2441 420 79
## HIV-Gagp1Pol 2120 1415 3167 794 129
## HIV-Polprot 28623 18993 46285 9177 1793
## HIV-Polp15p31 79443 56460 109411 15132 4996
## mdm_latent2 mdm_latent3
## HIV-Gagp17 2276 676
## HIV-Gagp24 0 0
## HIV-Gagp2p7 131 64
## HIV-Gagp1Pol 190 106
## HIV-Polprot 2050 1200
## HIV-Polp15p31 6607 3415pb1mf <- pb1m[which(rowMeans(pb1m)>=10),]
head(pb1mf)## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_latent1
## HIV-Gagp17 38988 23811 41474 15181 2226
## HIV-Gagp2p7 1512 1041 2441 420 79
## HIV-Gagp1Pol 2120 1415 3167 794 129
## HIV-Polprot 28623 18993 46285 9177 1793
## HIV-Polp15p31 79443 56460 109411 15132 4996
## HIV-Vif 5563 4336 7532 1123 323
## mdm_latent2 mdm_latent3
## HIV-Gagp17 2276 676
## HIV-Gagp2p7 131 64
## HIV-Gagp1Pol 190 106
## HIV-Polprot 2050 1200
## HIV-Polp15p31 6607 3415
## HIV-Vif 400 264colSums(pb1mf)## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_latent1 mdm_latent2
## 31082255 23021475 26484164 13953137 35438672 40374328
## mdm_latent3
## 15438893des1m <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 132nrow(subset(de1mf,padj<0.05 & log2FoldChange<0))## [1] 231head(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 | |
|---|---|---|---|---|---|---|
| HIV-TatEx1 | 1675.6625 | 4.342862 | 0.2147464 | 20.22321 | 0 | 0 |
| HIV-Gagp17 | 18173.6276 | 4.529783 | 0.2240585 | 20.21697 | 0 | 0 |
| HIV-TatEx2Rev | 506.6022 | 4.176781 | 0.2201533 | 18.97215 | 0 | 0 |
| HIV-Vpr | 236.0864 | 4.070712 | 0.2367387 | 17.19496 | 0 | 0 |
| HIV-EGFP | 581.3032 | 4.002777 | 0.2502409 | 15.99569 | 0 | 0 |
| HIV-Vpu | 407.7671 | 4.005511 | 0.2614884 | 15.31813 | 0 | 0 |
| HIV-BaLEnv | 24514.5004 | 4.023186 | 0.3136090 | 12.82867 | 0 | 0 |
| HIV-EnvStart | 344.4164 | 4.077739 | 0.3214085 | 12.68709 | 0 | 0 |
| HIV-Nef | 4223.7485 | 4.174131 | 0.3555051 | 11.74141 | 0 | 0 |
| HIV-Gagp1Pol | 1145.1277 | 3.951316 | 0.3402085 | 11.61440 | 0 | 0 |
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 | 843.20052 | -1.510758 | 0.2052774 | -7.359591 | 0 | 0.0e+00 |
| AL031123.1 | 530.93043 | -1.002340 | 0.1449762 | -6.913822 | 0 | 0.0e+00 |
| TGFBI | 483.93836 | -2.512506 | 0.3655489 | -6.873243 | 0 | 0.0e+00 |
| HIST1H1C | 225.17397 | -2.149168 | 0.3260517 | -6.591494 | 0 | 0.0e+00 |
| EVI2A | 578.47758 | -1.190354 | 0.1818460 | -6.545945 | 0 | 0.0e+00 |
| PTGDS | 29235.16003 | -2.300458 | 0.3620596 | -6.353812 | 0 | 1.0e-07 |
| MAP1A | 91.37518 | -1.980122 | 0.3246652 | -6.098967 | 0 | 7.0e-07 |
| FLRT2 | 314.72618 | -1.926588 | 0.3252583 | -5.923254 | 0 | 1.9e-06 |
| IGF1R | 390.87715 | -1.019752 | 0.1746670 | -5.838261 | 0 | 3.1e-06 |
| IGFBP7 | 211.39803 | -1.160608 | 0.1997703 | -5.809712 | 0 | 3.4e-06 |
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 = 15050## Note: no. genes in output = 15050## Note: estimated proportion of input genes in output = 1mres1m <- 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
## HIV-Gagp17 33137 27369 17430 571 2405
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 1239 1259 798 18 87
## HIV-Gagp1Pol 2199 2384 1723 31 160
## HIV-Polprot 24408 31056 23030 307 1666
## HIV-Polp15p31 40029 60765 43633 571 2886
## alv_latent3 alv_latent4
## HIV-Gagp17 1947 2780
## HIV-Gagp24 0 0
## HIV-Gagp2p7 129 146
## HIV-Gagp1Pol 220 271
## HIV-Polprot 2606 4055
## HIV-Polp15p31 5359 7157pb1af <- pb1a[which(rowMeans(pb1a)>=10),]
head(pb1af)## alv_active1 alv_active2 alv_active3 alv_latent1 alv_latent2
## HIV-Gagp17 33137 27369 17430 571 2405
## HIV-Gagp2p7 1239 1259 798 18 87
## HIV-Gagp1Pol 2199 2384 1723 31 160
## HIV-Polprot 24408 31056 23030 307 1666
## HIV-Polp15p31 40029 60765 43633 571 2886
## HIV-Vif 3260 4567 3261 38 212
## alv_latent3 alv_latent4
## HIV-Gagp17 1947 2780
## HIV-Gagp2p7 129 146
## HIV-Gagp1Pol 220 271
## HIV-Polprot 2606 4055
## HIV-Polp15p31 5359 7157
## HIV-Vif 397 535colSums(pb1af)## alv_active1 alv_active2 alv_active3 alv_latent1 alv_latent2 alv_latent3
## 30260812 28741617 24197964 15348688 39767313 50853187
## alv_latent4
## 47787737des1a <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 396nrow(subset(de1af,padj<0.05 & log2FoldChange<0))## [1] 267head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 13586.0895 | 4.319391 | 0.2290598 | 18.85705 | 0 | 0 |
| HIV-TatEx2Rev | 555.5079 | 4.195088 | 0.2410154 | 17.40589 | 0 | 0 |
| HIV-Gagp2p7 | 582.3286 | 4.167297 | 0.2669497 | 15.61079 | 0 | 0 |
| HIV-Gagp1Pol | 1118.5777 | 4.291624 | 0.2799387 | 15.33059 | 0 | 0 |
| HIV-Vpu | 435.0793 | 4.295050 | 0.3125677 | 13.74119 | 0 | 0 |
| HIV-BaLEnv | 24086.4951 | 4.458958 | 0.3392044 | 13.14534 | 0 | 0 |
| HIV-TatEx1 | 3138.2463 | 4.399587 | 0.3435753 | 12.80531 | 0 | 0 |
| HIV-Vpr | 456.2951 | 4.261135 | 0.3661304 | 11.63830 | 0 | 0 |
| HIV-Polprot | 14017.9578 | 4.334978 | 0.3986929 | 10.87298 | 0 | 0 |
| HIV-Polp15p31 | 25913.3306 | 4.340592 | 0.4067124 | 10.67239 | 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 | 887.92176 | -1.5471671 | 0.1592876 | -9.713041 | 0 | 0e+00 |
| LYZ | 54764.96581 | -1.1187143 | 0.1237738 | -9.038375 | 0 | 0e+00 |
| H1FX | 1085.85841 | -1.2152364 | 0.1636653 | -7.425132 | 0 | 0e+00 |
| JAKMIP2 | 75.29146 | -4.0707741 | 0.5566792 | -7.312603 | 0 | 0e+00 |
| P2RY11 | 192.30659 | -1.2331198 | 0.1713095 | -7.198197 | 0 | 0e+00 |
| IFNGR2 | 2429.88536 | -0.6571357 | 0.0981228 | -6.697076 | 0 | 0e+00 |
| CSF1R | 613.64876 | -1.6775637 | 0.2512201 | -6.677666 | 0 | 0e+00 |
| SLCO3A1 | 236.46470 | -1.7743334 | 0.2714364 | -6.536830 | 0 | 0e+00 |
| GCA | 1404.12051 | -0.9011225 | 0.1386440 | -6.499540 | 0 | 0e+00 |
| GCNT1 | 293.69222 | -1.1392712 | 0.1816482 | -6.271855 | 0 | 1e-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 = 16363## Note: no. genes in output = 16363## Note: estimated proportion of input genes in output = 1mres1a <- 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.1766999 0.2181777
## 2 A1BG-AS1 -0.5211198 -0.8128086
## 3 A2M 0.6933579 -0.4365208
## 4 A2M-AS1 -0.2377587 -1.1450259
## 5 A2ML1-AS1 -1.3628116 0.6944468
## 6 A4GALT 3.8897487 -0.8326164rownames(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
## -7.6842 -0.8263 -0.0705 0.7441 10.5087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.030914 0.011129 -2.778 0.00548 **
## MDM 0.491609 0.007687 63.952 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.351 on 14753 degrees of freedom
## Multiple R-squared: 0.2171, Adjusted R-squared: 0.217
## F-statistic: 4090 on 1 and 14753 DF, p-value: < 2.2e-16cor.test(mm1$Alv,mm1$MDM)##
## Pearson's product-moment correlation
##
## data: mm1$Alv and mm1$MDM
## t = 63.952, df = 14753, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4531582 0.4784263
## sample estimates:
## cor
## 0.4658873mm1r <- 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
## -9833.9 -3192.7 -40.5 3173.3 9997.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.390e+03 6.413e+01 68.46 <2e-16 ***
## MDM 4.049e-01 7.528e-03 53.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3895 on 14753 degrees of freedom
## Multiple R-squared: 0.164, Adjusted R-squared: 0.1639
## F-statistic: 2894 on 1 and 14753 DF, p-value: < 2.2e-16cor.test(mm1r$Alv,mm1r$MDM)##
## Pearson's product-moment correlation
##
## data: mm1r$Alv and mm1r$MDM
## t = 53.793, df = 14753, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3913650 0.4183456
## sample estimates:
## cor
## 0.4049435DE2 Latently-infected vs bystander cells
MDM group.
pb2m <- pbmdm[,c(grep("bystander",colnames(pbmdm)),grep("latent",colnames(pbmdm)))]
head(pb2m)## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1
## HIV-Gagp17 169 153 44 2226
## HIV-Gagp24 0 0 0 0
## HIV-Gagp2p7 18 12 4 79
## HIV-Gagp1Pol 14 23 8 129
## HIV-Polprot 183 211 114 1793
## HIV-Polp15p31 621 644 267 4996
## mdm_latent2 mdm_latent3
## HIV-Gagp17 2276 676
## HIV-Gagp24 0 0
## HIV-Gagp2p7 131 64
## HIV-Gagp1Pol 190 106
## HIV-Polprot 2050 1200
## HIV-Polp15p31 6607 3415pb2mf <- pb2m[which(rowMeans(pb2m)>=10),]
head(pb2mf)## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1
## HIV-Gagp17 169 153 44 2226
## HIV-Gagp2p7 18 12 4 79
## HIV-Gagp1Pol 14 23 8 129
## HIV-Polprot 183 211 114 1793
## HIV-Polp15p31 621 644 267 4996
## HIV-Vif 56 34 14 323
## mdm_latent2 mdm_latent3
## HIV-Gagp17 2276 676
## HIV-Gagp2p7 131 64
## HIV-Gagp1Pol 190 106
## HIV-Polprot 2050 1200
## HIV-Polp15p31 6607 3415
## HIV-Vif 400 264colSums(pb2mf)## mdm_bystander1 mdm_bystander2 mdm_bystander3 mdm_latent1 mdm_latent2
## 40413184 35971143 15858323 35442346 40379883
## mdm_latent3
## 15441023des2m <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 14nrow(subset(de2mf,padj<0.05 & log2FoldChange<0))## [1] 0head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 807.32055 | 3.805170 | 0.2315096 | 16.436341 | 0 | 0 |
| HIV-Polp15p31 | 2676.45023 | 3.360809 | 0.2116666 | 15.877842 | 0 | 0 |
| HIV-BaLEnv | 1544.75337 | 3.599384 | 0.2681426 | 13.423396 | 0 | 0 |
| HIV-Polprot | 916.10994 | 3.331844 | 0.2699042 | 12.344542 | 0 | 0 |
| HIV-TatEx1 | 87.68851 | 3.220525 | 0.3204557 | 10.049831 | 0 | 0 |
| HIV-Vif | 183.36367 | 3.425117 | 0.3693735 | 9.272774 | 0 | 0 |
| HIV-Nef | 248.74897 | 3.019120 | 0.3462927 | 8.718404 | 0 | 0 |
| HIV-Gagp1Pol | 77.58706 | 3.323194 | 0.4192376 | 7.926755 | 0 | 0 |
| HIV-Gagp2p7 | 49.41721 | 3.109492 | 0.4562406 | 6.815465 | 0 | 0 |
| FN1 | 370.66862 | 1.233991 | 0.1819136 | 6.783389 | 0 | 0 |
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 | |
|---|---|---|---|---|---|---|
| HAUS7 | 17.64648 | -1.6097132 | 0.5350538 | -3.008507 | 0.0026253 | 0.999993 |
| UBE2C | 13.84004 | -1.9696489 | 0.8528655 | -2.309448 | 0.0209187 | 0.999993 |
| NEGR1 | 15.58160 | -1.6421309 | 0.7128872 | -2.303494 | 0.0212511 | 0.999993 |
| MAP3K6 | 67.46090 | -0.6569070 | 0.2893879 | -2.269987 | 0.0232084 | 0.999993 |
| ACTG1 | 17106.05531 | -0.3266782 | 0.1452213 | -2.249520 | 0.0244794 | 0.999993 |
| FAM124B | 22.52463 | -0.9942079 | 0.4599870 | -2.161382 | 0.0306658 | 0.999993 |
| L1TD1 | 12.88960 | -1.7106347 | 0.7953757 | -2.150726 | 0.0314979 | 0.999993 |
| KIAA1324 | 13.72621 | -1.6095274 | 0.7497198 | -2.146838 | 0.0318061 | 0.999993 |
| HSPA5 | 3142.40439 | -0.3903567 | 0.1838809 | -2.122878 | 0.0337641 | 0.999993 |
| NOXA1 | 100.56679 | -0.5152430 | 0.2445087 | -2.107258 | 0.0350952 | 0.999993 |
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() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 13nrow(subset(de2mf,padj<0.05 & log2FoldChange<0))## [1] 0head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 807.32055 | 3.814767 | 0.1454306 | 26.230839 | 0 | 0 |
| HIV-Polp15p31 | 2676.45023 | 3.329784 | 0.1520529 | 21.898850 | 0 | 0 |
| HIV-BaLEnv | 1544.75337 | 3.547144 | 0.1823231 | 19.455261 | 0 | 0 |
| HIV-Polprot | 916.10994 | 3.321343 | 0.1785778 | 18.598860 | 0 | 0 |
| HIV-Nef | 248.74897 | 3.057518 | 0.2524389 | 12.111916 | 0 | 0 |
| HIV-TatEx1 | 87.68851 | 3.211806 | 0.3244026 | 9.900677 | 0 | 0 |
| HIV-Vif | 183.36367 | 3.372521 | 0.3575460 | 9.432410 | 0 | 0 |
| HIV-Gagp1Pol | 77.58706 | 3.247310 | 0.3482884 | 9.323624 | 0 | 0 |
| FN1 | 370.66862 | 1.236975 | 0.1747896 | 7.076937 | 0 | 0 |
| HIV-EGFP | 39.12916 | 3.259385 | 0.4742418 | 6.872833 | 0 | 0 |
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 | |
|---|---|---|---|---|---|---|
| HAUS7 | 17.64648 | -1.6257934 | 0.5660163 | -2.872344 | 0.0040744 | 0.9998533 |
| VSIG4 | 137.94402 | -0.5710272 | 0.2171643 | -2.629471 | 0.0085518 | 0.9998533 |
| ACTG1 | 17106.05531 | -0.3295194 | 0.1261737 | -2.611634 | 0.0090111 | 0.9998533 |
| HSPA5 | 3142.40439 | -0.3746313 | 0.1522283 | -2.460983 | 0.0138557 | 0.9998533 |
| SNHG29 | 6662.25703 | -0.1935706 | 0.0811983 | -2.383924 | 0.0171291 | 0.9998533 |
| CASP4 | 1537.73005 | -0.2404318 | 0.1014503 | -2.369947 | 0.0177906 | 0.9998533 |
| NEGR1 | 15.58160 | -1.5199124 | 0.6495528 | -2.339937 | 0.0192870 | 0.9998533 |
| MAD2L2 | 952.61966 | -0.2524462 | 0.1085408 | -2.325818 | 0.0200283 | 0.9998533 |
| MAP3K6 | 67.46090 | -0.6446851 | 0.2822796 | -2.283853 | 0.0223802 | 0.9998533 |
| KIAA1324 | 13.72621 | -1.7619241 | 0.7803713 | -2.257802 | 0.0239580 | 0.9998533 |
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 = 15341## Note: no. genes in output = 15341## Note: estimated proportion of input genes in output = 1mres2m <- 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
## HIV-Gagp17 38 45 37 49
## HIV-Gagp24 0 0 0 0
## HIV-Gagp2p7 4 4 5 7
## HIV-Gagp1Pol 6 9 18 13
## HIV-Polprot 40 63 106 122
## HIV-Polp15p31 93 183 271 263
## alv_latent1 alv_latent2 alv_latent3 alv_latent4
## HIV-Gagp17 571 2405 1947 2780
## HIV-Gagp24 0 0 0 0
## HIV-Gagp2p7 18 87 129 146
## HIV-Gagp1Pol 31 160 220 271
## HIV-Polprot 307 1666 2606 4055
## HIV-Polp15p31 571 2886 5359 7157pb2af <- pb2a[which(rowMeans(pb2a)>=10),]
head(pb2af)## alv_bystander1 alv_bystander2 alv_bystander3 alv_bystander4
## HIV-Gagp17 38 45 37 49
## HIV-Gagp2p7 4 4 5 7
## HIV-Gagp1Pol 6 9 18 13
## HIV-Polprot 40 63 106 122
## HIV-Polp15p31 93 183 271 263
## HIV-Vif 9 7 24 22
## alv_latent1 alv_latent2 alv_latent3 alv_latent4
## HIV-Gagp17 571 2405 1947 2780
## HIV-Gagp2p7 18 87 129 146
## HIV-Gagp1Pol 31 160 220 271
## HIV-Polprot 307 1666 2606 4055
## HIV-Polp15p31 571 2886 5359 7157
## HIV-Vif 38 212 397 535colSums(pb2af)## alv_bystander1 alv_bystander2 alv_bystander3 alv_bystander4 alv_latent1
## 22135810 26459262 20201267 18043458 15347444
## alv_latent2 alv_latent3 alv_latent4
## 39761925 50848809 47782839des2a <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 14nrow(subset(de2af,padj<0.05 & log2FoldChange<0))## [1] 0head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 682.23838 | 4.587657 | 0.2553548 | 17.965816 | 0 | 0.00e+00 |
| HIV-BaLEnv | 1126.97572 | 3.973068 | 0.3506990 | 11.328995 | 0 | 0.00e+00 |
| HIV-TatEx1 | 158.62905 | 3.183806 | 0.4097671 | 7.769793 | 0 | 0.00e+00 |
| HIV-Polprot | 723.24693 | 3.602250 | 0.4856452 | 7.417452 | 0 | 0.00e+00 |
| HIV-Polp15p31 | 1363.63683 | 3.214320 | 0.4595271 | 6.994844 | 0 | 0.00e+00 |
| HIV-Nef | 960.79384 | 3.136794 | 0.4627469 | 6.778639 | 0 | 0.00e+00 |
| HIV-Gagp2p7 | 33.46548 | 3.243846 | 0.4907472 | 6.610013 | 0 | 1.00e-07 |
| HIV-Gagp1Pol | 61.53005 | 2.867191 | 0.4475297 | 6.406706 | 0 | 3.00e-07 |
| HIV-TatEx2Rev | 31.14635 | 3.293310 | 0.5210228 | 6.320855 | 0 | 5.00e-07 |
| HIV-Vif | 100.89948 | 3.122225 | 0.5657412 | 5.518822 | 0 | 5.44e-05 |
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 | |
|---|---|---|---|---|---|---|
| TOP2A | 93.87932 | -0.7914506 | 0.3207406 | -2.467572 | 0.0136033 | 0.9999184 |
| PPP1R9A | 63.81716 | -0.5345124 | 0.2335402 | -2.288738 | 0.0220946 | 0.9999184 |
| GTSE1 | 39.44756 | -1.0381862 | 0.4573980 | -2.269766 | 0.0232218 | 0.9999184 |
| AL139041.1 | 10.57704 | -1.1314320 | 0.5144665 | -2.199234 | 0.0278613 | 0.9999184 |
| CEP55 | 68.61181 | -0.7744077 | 0.3560814 | -2.174806 | 0.0296447 | 0.9999184 |
| GALK1 | 855.92475 | -0.2345995 | 0.1085333 | -2.161545 | 0.0306533 | 0.9999184 |
| KIF11 | 58.51508 | -0.6952549 | 0.3220384 | -2.158920 | 0.0308564 | 0.9999184 |
| METTL18 | 60.52893 | -0.5017172 | 0.2402173 | -2.088597 | 0.0367440 | 0.9999184 |
| HJURP | 11.46129 | -1.3801405 | 0.6740506 | -2.047532 | 0.0406058 | 0.9999184 |
| CTSW | 53.58878 | -0.5759731 | 0.2870709 | -2.006380 | 0.0448158 | 0.9999184 |
des2a$sample <- rep(1:4,2)
dds <- DESeqDataSetFromMatrix(countData = pb2af , colData = des2a, design = ~ case)## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 14nrow(subset(de2af,padj<0.05 & log2FoldChange<0))## [1] 0head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 682.23838 | 4.587657 | 0.2553548 | 17.965816 | 0 | 0.00e+00 |
| HIV-BaLEnv | 1126.97572 | 3.973068 | 0.3506990 | 11.328995 | 0 | 0.00e+00 |
| HIV-TatEx1 | 158.62905 | 3.183806 | 0.4097671 | 7.769793 | 0 | 0.00e+00 |
| HIV-Polprot | 723.24693 | 3.602250 | 0.4856452 | 7.417452 | 0 | 0.00e+00 |
| HIV-Polp15p31 | 1363.63683 | 3.214320 | 0.4595271 | 6.994844 | 0 | 0.00e+00 |
| HIV-Nef | 960.79384 | 3.136794 | 0.4627469 | 6.778639 | 0 | 0.00e+00 |
| HIV-Gagp2p7 | 33.46548 | 3.243846 | 0.4907472 | 6.610013 | 0 | 1.00e-07 |
| HIV-Gagp1Pol | 61.53005 | 2.867191 | 0.4475297 | 6.406706 | 0 | 3.00e-07 |
| HIV-TatEx2Rev | 31.14635 | 3.293310 | 0.5210228 | 6.320855 | 0 | 5.00e-07 |
| HIV-Vif | 100.89948 | 3.122225 | 0.5657412 | 5.518822 | 0 | 5.44e-05 |
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 | |
|---|---|---|---|---|---|---|
| TOP2A | 93.87932 | -0.7914506 | 0.3207406 | -2.467572 | 0.0136033 | 0.9999184 |
| PPP1R9A | 63.81716 | -0.5345124 | 0.2335402 | -2.288738 | 0.0220946 | 0.9999184 |
| GTSE1 | 39.44756 | -1.0381862 | 0.4573980 | -2.269766 | 0.0232218 | 0.9999184 |
| AL139041.1 | 10.57704 | -1.1314320 | 0.5144665 | -2.199234 | 0.0278613 | 0.9999184 |
| CEP55 | 68.61181 | -0.7744077 | 0.3560814 | -2.174806 | 0.0296447 | 0.9999184 |
| GALK1 | 855.92475 | -0.2345995 | 0.1085333 | -2.161545 | 0.0306533 | 0.9999184 |
| KIF11 | 58.51508 | -0.6952549 | 0.3220384 | -2.158920 | 0.0308564 | 0.9999184 |
| METTL18 | 60.52893 | -0.5017172 | 0.2402173 | -2.088597 | 0.0367440 | 0.9999184 |
| HJURP | 11.46129 | -1.3801405 | 0.6740506 | -2.047532 | 0.0406058 | 0.9999184 |
| CTSW | 53.58878 | -0.5759731 | 0.2870709 | -2.006380 | 0.0448158 | 0.9999184 |
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 = 15993## Note: no. genes in output = 15993## Note: estimated proportion of input genes in output = 1mres2a <- 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.07488022 -1.1629136
## 2 A1BG-AS1 -0.04236658 0.1695308
## 3 A2M -0.58712217 0.2715816
## 4 A2M-AS1 -0.24589352 -0.1005539
## 5 A2ML1-AS1 -0.65417214 -0.2578660
## 6 AAAS -0.22837395 -0.7255075rownames(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.4473 -0.2593 -0.0047 0.2505 11.1454
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.014684 0.003945 -3.722 0.000198 ***
## MDM 0.260576 0.005076 51.340 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4803 on 14828 degrees of freedom
## Multiple R-squared: 0.1509, Adjusted R-squared: 0.1509
## F-statistic: 2636 on 1 and 14828 DF, p-value: < 2.2e-16cor.test(mm2$Alv,mm2$MDM)##
## Pearson's product-moment correlation
##
## data: mm2$Alv and mm2$MDM
## t = 51.34, df = 14828, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3747418 0.4020741
## sample estimates:
## cor
## 0.3884934mm2r <- 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
## -9037.6 -3450.6 57.9 3452.6 9078.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.721e+03 6.846e+01 83.57 <2e-16 ***
## MDM 2.285e-01 7.995e-03 28.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4168 on 14828 degrees of freedom
## Multiple R-squared: 0.0522, Adjusted R-squared: 0.05214
## F-statistic: 816.7 on 1 and 14828 DF, p-value: < 2.2e-16cor.test(mm2r$Alv,mm2r$MDM)##
## Pearson's product-moment correlation
##
## data: mm2r$Alv and mm2r$MDM
## t = 28.578, df = 14828, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2131706 0.2436801
## sample estimates:
## cor
## 0.2284814DE3 Active vs mock
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
## HIV-Gagp17 38988 23811 41474 15181 135
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 1512 1041 2441 420 8
## HIV-Gagp1Pol 2120 1415 3167 794 8
## HIV-Polprot 28623 18993 46285 9177 159
## HIV-Polp15p31 79443 56460 109411 15132 427
## mdm_mock2 mdm_mock3 mdm_mock4
## HIV-Gagp17 271 41 161
## HIV-Gagp24 0 0 0
## HIV-Gagp2p7 8 1 5
## HIV-Gagp1Pol 14 3 13
## HIV-Polprot 216 51 147
## HIV-Polp15p31 621 109 259pb3mf <- pb3m[which(rowMeans(pb3m)>=10),]
head(pb3mf)## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_mock1
## HIV-Gagp17 38988 23811 41474 15181 135
## HIV-Gagp2p7 1512 1041 2441 420 8
## HIV-Gagp1Pol 2120 1415 3167 794 8
## HIV-Polprot 28623 18993 46285 9177 159
## HIV-Polp15p31 79443 56460 109411 15132 427
## HIV-Vif 5563 4336 7532 1123 24
## mdm_mock2 mdm_mock3 mdm_mock4
## HIV-Gagp17 271 41 161
## HIV-Gagp2p7 8 1 5
## HIV-Gagp1Pol 14 3 13
## HIV-Polprot 216 51 147
## HIV-Polp15p31 621 109 259
## HIV-Vif 51 11 17colSums(pb3mf)## mdm_active1 mdm_active2 mdm_active3 mdm_active4 mdm_mock1 mdm_mock2
## 31076648 23016763 26480052 13951080 29174508 21631818
## mdm_mock3 mdm_mock4
## 7551959 20738346des3m <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 126nrow(subset(de3mf,padj<0.05 & log2FoldChange<0))## [1] 233head(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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 12757.7640 | 7.432858 | 0.3050026 | 24.36982 | 0 | 0 |
| HIV-BaLEnv | 16996.1820 | 7.264108 | 0.3244289 | 22.39045 | 0 | 0 |
| HIV-Polprot | 10745.0888 | 7.235320 | 0.3344748 | 21.63188 | 0 | 0 |
| HIV-Gagp1Pol | 791.1831 | 7.380796 | 0.3974803 | 18.56896 | 0 | 0 |
| HIV-TatEx1 | 1172.5458 | 7.012038 | 0.3801895 | 18.44353 | 0 | 0 |
| HIV-TatEx2Rev | 352.2232 | 7.210588 | 0.3928546 | 18.35434 | 0 | 0 |
| HIV-Polp15p31 | 26459.4331 | 7.273481 | 0.4014867 | 18.11637 | 0 | 0 |
| HIV-EGFP | 401.8456 | 7.365772 | 0.4117606 | 17.88848 | 0 | 0 |
| HIV-Vpu | 286.1028 | 5.559431 | 0.3279902 | 16.94999 | 0 | 0 |
| HIV-EnvStart | 239.4098 | 6.724076 | 0.3973697 | 16.92146 | 0 | 0 |
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 | 843.20052 | -1.510758 | 0.2052774 | -7.359591 | 0 | 0.0e+00 |
| AL031123.1 | 530.93043 | -1.002340 | 0.1449762 | -6.913822 | 0 | 0.0e+00 |
| TGFBI | 483.93836 | -2.512506 | 0.3655489 | -6.873243 | 0 | 0.0e+00 |
| HIST1H1C | 225.17397 | -2.149168 | 0.3260517 | -6.591494 | 0 | 0.0e+00 |
| EVI2A | 578.47758 | -1.190354 | 0.1818460 | -6.545945 | 0 | 0.0e+00 |
| PTGDS | 29235.16003 | -2.300458 | 0.3620596 | -6.353812 | 0 | 1.0e-07 |
| MAP1A | 91.37518 | -1.980122 | 0.3246652 | -6.098967 | 0 | 7.0e-07 |
| FLRT2 | 314.72618 | -1.926588 | 0.3252583 | -5.923254 | 0 | 1.9e-06 |
| IGF1R | 390.87715 | -1.019752 | 0.1746670 | -5.838261 | 0 | 3.1e-06 |
| IGFBP7 | 211.39803 | -1.160608 | 0.1997703 | -5.809712 | 0 | 3.4e-06 |
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 = 14627## Note: no. genes in output = 14627## Note: estimated proportion of input genes in output = 1mres3m <- 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
## HIV-Gagp17 33137 27369 17430 107 186 1541
## HIV-Gagp24 0 0 0 0 0 0
## HIV-Gagp2p7 1239 1259 798 2 8 53
## HIV-Gagp1Pol 2199 2384 1723 6 23 97
## HIV-Polprot 24408 31056 23030 96 246 1630
## HIV-Polp15p31 40029 60765 43633 172 415 2890pb3af <- pb3a[which(rowMeans(pb3a)>=10),]
head(pb3af)## alv_active1 alv_active2 alv_active3 alv_mock1 alv_mock2 alv_mock3
## HIV-Gagp17 33137 27369 17430 107 186 1541
## HIV-Gagp2p7 1239 1259 798 2 8 53
## HIV-Gagp1Pol 2199 2384 1723 6 23 97
## HIV-Polprot 24408 31056 23030 96 246 1630
## HIV-Polp15p31 40029 60765 43633 172 415 2890
## HIV-Vif 3260 4567 3261 10 38 170colSums(pb3af)## alv_active1 alv_active2 alv_active3 alv_mock1 alv_mock2 alv_mock3
## 30254254 28737263 24194114 20460289 25558226 34906906des3a <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 910nrow(subset(de3af,padj<0.05 & log2FoldChange<0))## [1] 657head(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 | |
|---|---|---|---|---|---|---|
| HIV-EGFP | 526.4017 | 6.235724 | 0.4312771 | 14.45874 | 0 | 0 |
| AL157912.1 | 754.2726 | 2.732506 | 0.2101087 | 13.00520 | 0 | 0 |
| LAYN | 553.0524 | 2.260727 | 0.1867841 | 12.10342 | 0 | 0 |
| HIV-Vpu | 418.4397 | 5.500942 | 0.4733986 | 11.62011 | 0 | 0 |
| HES4 | 381.7114 | 2.356797 | 0.2048937 | 11.50253 | 0 | 0 |
| CDKN1A | 5305.2130 | 1.554187 | 0.1353865 | 11.47962 | 0 | 0 |
| OCIAD2 | 353.0509 | 2.481124 | 0.2216815 | 11.19229 | 0 | 0 |
| SDS | 4765.1967 | 2.107095 | 0.1956164 | 10.77157 | 0 | 0 |
| IL6R-AS1 | 455.0925 | 3.519844 | 0.3314263 | 10.62029 | 0 | 0 |
| TNFRSF9 | 170.2840 | 2.517982 | 0.2511971 | 10.02393 | 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.3566 | -1.805140 | 0.1647590 | -10.956244 | 0 | 0 |
| HIST1H1C | 791.2781 | -1.157674 | 0.1092179 | -10.599670 | 0 | 0 |
| ADA2 | 1790.5726 | -1.042049 | 0.1149916 | -9.061955 | 0 | 0 |
| CEBPD | 686.4738 | -1.475359 | 0.1637708 | -9.008680 | 0 | 0 |
| CRIM1 | 387.7032 | -1.675237 | 0.1941417 | -8.628942 | 0 | 0 |
| LYZ | 45988.3407 | -1.176674 | 0.1384751 | -8.497367 | 0 | 0 |
| RPL10A | 8112.3056 | -0.869581 | 0.1032653 | -8.420840 | 0 | 0 |
| CHST13 | 459.5533 | -1.623448 | 0.1970857 | -8.237269 | 0 | 0 |
| TGFBI | 395.5933 | -2.330214 | 0.2883685 | -8.080684 | 0 | 0 |
| PIK3CD | 584.8424 | -1.034669 | 0.1375107 | -7.524284 | 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 = 15785## Note: no. genes in output = 15785## Note: estimated proportion of input genes in output = 1mres3a <- 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.2596782 0.37934188
## 2 A1BG-AS1 -0.6990158 -1.19192851
## 3 A2M -1.1374479 -0.46329878
## 4 A2ML1-AS1 -0.3602939 0.37546545
## 5 A4GALT 3.0255735 0.06952938
## 6 AAAS 0.7912558 -0.37305290rownames(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
## -13.1001 -0.8694 -0.0950 0.7791 11.5619
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.085433 0.012104 7.058 1.76e-12 ***
## MDM 0.768818 0.008419 91.324 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.451 on 14365 degrees of freedom
## Multiple R-squared: 0.3673, Adjusted R-squared: 0.3673
## F-statistic: 8340 on 1 and 14365 DF, p-value: < 2.2e-16cor.test(mm3$Alv,mm3$MDM)##
## Pearson's product-moment correlation
##
## data: mm3$Alv and mm3$MDM
## t = 91.324, df = 14365, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5956226 0.6163157
## sample estimates:
## cor
## 0.6060717mm3r <- 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
## -10079 -2490 -79 2419 11034
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.778e+03 5.467e+01 50.82 <2e-16 ***
## MDM 6.133e-01 6.590e-03 93.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3276 on 14365 degrees of freedom
## Multiple R-squared: 0.3762, Adjusted R-squared: 0.3761
## F-statistic: 8662 on 1 and 14365 DF, p-value: < 2.2e-16cor.test(mm3r$Alv,mm3r$MDM)##
## Pearson's product-moment correlation
##
## data: mm3r$Alv and mm3r$MDM
## t = 93.069, df = 14365, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6030165 0.6234206
## sample estimates:
## cor
## 0.6133209DE4 Mock vs bystander
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
## HIV-Gagp17 135 271 41 161 169
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 8 8 1 5 18
## HIV-Gagp1Pol 8 14 3 13 14
## HIV-Polprot 159 216 51 147 183
## HIV-Polp15p31 427 621 109 259 621
## mdm_bystander2 mdm_bystander3
## HIV-Gagp17 153 44
## HIV-Gagp24 0 0
## HIV-Gagp2p7 12 4
## HIV-Gagp1Pol 23 8
## HIV-Polprot 211 114
## HIV-Polp15p31 644 267pb4mf <- pb4m[which(rowMeans(pb4m)>=10),]
head(pb4mf)## mdm_mock1 mdm_mock2 mdm_mock3 mdm_mock4 mdm_bystander1
## HIV-Gagp17 135 271 41 161 169
## HIV-Gagp1Pol 8 14 3 13 14
## HIV-Polprot 159 216 51 147 183
## HIV-Polp15p31 427 621 109 259 621
## HIV-Vif 24 51 11 17 56
## HIV-TatEx1 21 14 2 39 19
## mdm_bystander2 mdm_bystander3
## HIV-Gagp17 153 44
## HIV-Gagp1Pol 23 8
## HIV-Polprot 211 114
## HIV-Polp15p31 644 267
## HIV-Vif 34 14
## HIV-TatEx1 27 9colSums(pb4mf)## mdm_mock1 mdm_mock2 mdm_mock3 mdm_mock4 mdm_bystander1
## 29177712 21635063 7553184 20740135 40406000
## mdm_bystander2 mdm_bystander3
## 35964277 15855580des4m <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 1nrow(subset(de4mf,padj<0.05 & log2FoldChange<0))## [1] 8head(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.434153e+02 | 0.9688453 | 0.2283164 | 4.243433 | 0.0000220 | 0.0419342 |
| IFI6 | 1.226532e+04 | 1.0668393 | 0.2729669 | 3.908310 | 0.0000929 | 0.1059490 |
| CCL8 | 9.050177e+00 | 3.2167923 | 0.8683945 | 3.704298 | 0.0002120 | 0.1847813 |
| ISG15 | 9.087671e+02 | 0.7233303 | 0.2053478 | 3.522464 | 0.0004276 | 0.2754754 |
| PSME2 | 3.426673e+03 | 0.3935419 | 0.1135973 | 3.464359 | 0.0005315 | 0.3029325 |
| PLAAT3 | 2.118101e+01 | 1.5674303 | 0.4541629 | 3.451251 | 0.0005580 | 0.3062559 |
| TNFAIP6 | 1.216433e+02 | 1.4444557 | 0.4238186 | 3.408193 | 0.0006539 | 0.3341660 |
| LINC01705 | 1.345073e+02 | 1.2371128 | 0.4028518 | 3.070888 | 0.0021342 | 0.7355154 |
| APOE | 2.869035e+05 | 0.3093318 | 0.1019077 | 3.035410 | 0.0024021 | 0.7910351 |
| PSMC4 | 1.060944e+03 | 0.4357276 | 0.1443698 | 3.018135 | 0.0025434 | 0.8193475 |
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 | |
|---|---|---|---|---|---|---|
| CEP55 | 47.15237 | -2.007242 | 0.3898421 | -5.148859 | 3.00e-07 | 0.0038837 |
| CCL4L2 | 38.05543 | -4.060426 | 0.8222305 | -4.938306 | 8.00e-07 | 0.0058390 |
| CENPF | 129.49077 | -2.491816 | 0.5340136 | -4.666204 | 3.10e-06 | 0.0151557 |
| PBK | 16.92059 | -5.615874 | 1.2515739 | -4.487050 | 7.20e-06 | 0.0267543 |
| GINS2 | 31.46686 | -2.171725 | 0.5018815 | -4.327167 | 1.51e-05 | 0.0419342 |
| HELLS | 66.97033 | -1.526082 | 0.3560429 | -4.286231 | 1.82e-05 | 0.0419342 |
| TK1 | 91.40061 | -1.648684 | 0.3891027 | -4.237144 | 2.26e-05 | 0.0419342 |
| CIT | 21.61889 | -2.359432 | 0.5613257 | -4.203322 | 2.63e-05 | 0.0433088 |
| RRM2 | 51.73552 | -2.021116 | 0.4952034 | -4.081386 | 4.48e-05 | 0.0663417 |
| DIAPH3 | 27.63781 | -2.255145 | 0.5680445 | -3.970014 | 7.19e-05 | 0.0917829 |
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 = 14841## Note: no. genes in output = 14841## Note: estimated proportion of input genes in output = 1mres4m <- 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
## HIV-Gagp17 107 186 1541 38 45
## HIV-Gagp24 0 0 0 0 0
## HIV-Gagp2p7 2 8 53 4 4
## HIV-Gagp1Pol 6 23 97 6 9
## HIV-Polprot 96 246 1630 40 63
## HIV-Polp15p31 172 415 2890 93 183
## alv_bystander3 alv_bystander4
## HIV-Gagp17 37 49
## HIV-Gagp24 0 0
## HIV-Gagp2p7 5 7
## HIV-Gagp1Pol 18 13
## HIV-Polprot 106 122
## HIV-Polp15p31 271 263pb4af <- pb4a[which(rowMeans(pb4a)>=10),]
head(pb4af)## alv_mock1 alv_mock2 alv_mock3 alv_bystander1 alv_bystander2
## HIV-Gagp17 107 186 1541 38 45
## HIV-Gagp2p7 2 8 53 4 4
## HIV-Gagp1Pol 6 23 97 6 9
## HIV-Polprot 96 246 1630 40 63
## HIV-Polp15p31 172 415 2890 93 183
## HIV-Vif 10 38 170 9 7
## alv_bystander3 alv_bystander4
## HIV-Gagp17 37 49
## HIV-Gagp2p7 5 7
## HIV-Gagp1Pol 18 13
## HIV-Polprot 106 122
## HIV-Polp15p31 271 263
## HIV-Vif 24 22colSums(pb4af)## alv_mock1 alv_mock2 alv_mock3 alv_bystander1 alv_bystander2
## 20457270 25554735 34901831 22131109 26451508
## alv_bystander3 alv_bystander4
## 20197319 18039683des4a <- 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 moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
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] 5nrow(subset(de4af,padj<0.05 & log2FoldChange<0))## [1] 2head(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 | 459.3954 | 1.8606919 | 0.3242310 | 5.738785 | 0.0e+00 | 0.0001320 |
| EPSTI1 | 472.1009 | 2.9143141 | 0.5316064 | 5.482090 | 0.0e+00 | 0.0002910 |
| IFI6 | 14927.6263 | 2.3027078 | 0.4436679 | 5.190161 | 2.0e-07 | 0.0009696 |
| STAT1 | 2607.8926 | 0.8949205 | 0.1996408 | 4.482653 | 7.4e-06 | 0.0204118 |
| TRIM22 | 545.5101 | 0.9162842 | 0.2101830 | 4.359459 | 1.3e-05 | 0.0300840 |
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 | |
|---|---|---|---|---|---|---|
| HIV-Gagp17 | 216.30501 | -3.2208126 | 0.6560128 | -4.909680 | 0.0000009 | 0.0031573 |
| F13A1 | 34.61936 | -5.0096404 | 1.1666584 | -4.294008 | 0.0000175 | 0.0347042 |
| HIV-BaLEnv | 368.61668 | -2.5961220 | 0.6408586 | -4.051006 | 0.0000510 | 0.0737365 |
| SPP1 | 45874.10919 | -0.3908419 | 0.0967232 | -4.040830 | 0.0000533 | 0.0737365 |
| FCGBP | 23.86648 | -4.1277493 | 1.1053562 | -3.734316 | 0.0001882 | 0.1491341 |
| MTSS1 | 220.89652 | -0.7087421 | 0.2090415 | -3.390437 | 0.0006978 | 0.3715582 |
| HIV-Nef | 311.43755 | -1.5661889 | 0.5565793 | -2.813954 | 0.0048936 | 0.9998366 |
| FAM135B | 12.70228 | -1.8689513 | 0.6773992 | -2.759010 | 0.0057977 | NA |
| SCN9A | 21.67957 | -1.0850854 | 0.4048652 | -2.680115 | 0.0073597 | 0.9998366 |
| HIV-Polprot | 258.24846 | -2.3234129 | 0.9116187 | -2.548667 | 0.0108135 | 0.9998366 |
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 = 15359## Note: no. genes in output = 15359## Note: estimated proportion of input genes in output = 1mres4a <- 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.04744125 0.40612769
## 2 A1BG-AS1 -0.15136471 -0.46385822
## 3 A2M -1.39090885 -0.20667761
## 4 A2ML1-AS1 1.77490432 0.01077993
## 5 AAAS 0.74618958 -0.01659174
## 6 AACS 0.19712619 -0.69335224rownames(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.8148 -0.4849 -0.0352 0.4657 5.5277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.105910 0.005792 18.29 <2e-16 ***
## MDM 0.081008 0.006345 12.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6918 on 14370 degrees of freedom
## Multiple R-squared: 0.01122, Adjusted R-squared: 0.01115
## F-statistic: 163 on 1 and 14370 DF, p-value: < 2.2e-16cor.test(mm4$Alv,mm4$MDM)##
## Pearson's product-moment correlation
##
## data: mm4$Alv and mm4$MDM
## t = 12.767, df = 14370, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08970817 0.12203992
## sample estimates:
## cor
## 0.105902mm4r <- 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
## -7663.1 -3592.7 5.6 3600.6 7695.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.630e+03 6.902e+01 96.065 <2e-16 ***
## MDM 7.744e-02 8.317e-03 9.312 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4137 on 14370 degrees of freedom
## Multiple R-squared: 0.005998, Adjusted R-squared: 0.005928
## F-statistic: 86.7 on 1 and 14370 DF, p-value: < 2.2e-16cor.test(mm4r$Alv,mm4r$MDM)##
## Pearson's product-moment correlation
##
## data: mm4r$Alv and mm4r$MDM
## t = 9.3115, df = 14370, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06117181 0.09367413
## sample estimates:
## cor
## 0.07744355Cross comparison pathway enrichment heatmap
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 = 15419.875## Note: no. genes in output = 14060## Note: estimated proportion of input genes in output = 0.912lmres <- 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2043 | GO:0019773 proteasome core complex, alpha-subunit complex | 7 | 0.0000004 | -0.0236553 | 0.3842291 | -0.6484330 | -0.9203220 | 0.2721026 | 0.7910563 | 0.8458489 | 0.8398715 | 0.9137055 | 0.0783636 | 0.0029687 | 0.0000247 | 0.2125649 | 0.0002893 | 0.0001062 | 0.0001188 | 1.880506 | 0.6802925 | 0.0000172 |
| 3692 | GO:0045656 negative regulation of monocyte differentiation | 5 | 0.0041548 | -0.7651512 | -0.4216720 | 0.4028602 | 0.7613945 | -0.8781644 | -0.6022768 | -0.4549413 | -0.5978655 | 0.0030457 | 0.1025173 | 0.1187783 | 0.0031928 | 0.0006716 | 0.0196896 | 0.0781357 | 0.0206060 | 1.789862 | 0.5839620 | 0.0408788 |
| 2055 | GO:0019864 IgG binding | 6 | 0.0001792 | -0.8054646 | -0.5805227 | -0.6754423 | -0.4893032 | -0.7538779 | -0.7463593 | 0.5432854 | 0.1523173 | 0.0006332 | 0.0137984 | 0.0041673 | 0.0379431 | 0.0013836 | 0.0015448 | 0.0211954 | 0.5182555 | 1.767813 | 0.4953862 | 0.0034676 |
| 3120 | GO:0036402 proteasome-activating activity | 6 | 0.0002016 | 0.1303306 | 0.5216783 | -0.8234429 | -0.6095773 | 0.1662872 | 0.8219487 | 0.6671173 | 0.5398226 | 0.5804128 | 0.0269103 | 0.0004771 | 0.0097169 | 0.4806300 | 0.0004886 | 0.0046565 | 0.0220332 | 1.666878 | 0.6010110 | 0.0037661 |
| 2593 | GO:0032395 MHC class II receptor activity | 6 | 0.0000001 | -0.8578578 | -0.4018310 | 0.1870879 | -0.1942507 | -0.3388834 | -0.7122290 | 0.9674588 | -0.4596082 | 0.0002733 | 0.0883069 | 0.4274774 | 0.4099981 | 0.1506105 | 0.0025168 | 0.0000405 | 0.0512364 | 1.655129 | 0.5769239 | 0.0000027 |
| 811 | GO:0005839 proteasome core complex | 18 | 0.0000000 | -0.1036256 | 0.1667695 | -0.6676004 | -0.8221763 | 0.2077972 | 0.4841112 | 0.8100935 | 0.7624666 | 0.4466915 | 0.2207109 | 0.0000009 | 0.0000000 | 0.1270245 | 0.0003770 | 0.0000000 | 0.0000000 | 1.635654 | 0.6079965 | 0.0000000 |
| 4621 | GO:0070508 cholesterol import | 5 | 0.0022323 | 0.7318534 | 0.7695055 | -0.1760085 | 0.1607542 | 0.7475916 | 0.6210317 | 0.4609747 | -0.5562291 | 0.0045954 | 0.0028829 | 0.4955539 | 0.5336549 | 0.0037908 | 0.0161791 | 0.0742658 | 0.0312504 | 1.628199 | 0.4926886 | 0.0260713 |
| 1458 | GO:0008541 proteasome regulatory particle, lid subcomplex | 8 | 0.0000686 | 0.1359949 | 0.5980288 | -0.7048641 | -0.6677875 | 0.2514233 | 0.7700861 | 0.6241638 | 0.5337852 | 0.5054158 | 0.0033996 | 0.0005553 | 0.0010721 | 0.2182141 | 0.0001618 | 0.0022347 | 0.0089399 | 1.627778 | 0.5797647 | 0.0015285 |
| 3279 | GO:0042613 MHC class II protein complex | 12 | 0.0000000 | -0.7915718 | -0.5520952 | -0.0209876 | -0.4147565 | -0.3786660 | -0.7293090 | 0.8711916 | -0.1914389 | 0.0000020 | 0.0009279 | 0.8998424 | 0.0128634 | 0.0231450 | 0.0000121 | 0.0000002 | 0.2509406 | 1.604606 | 0.5298997 | 0.0000000 |
| 374 | GO:0002503 peptide antigen assembly with MHC class II protein complex | 11 | 0.0000000 | -0.7761083 | -0.5601499 | -0.0385728 | -0.4351846 | -0.3381347 | -0.7180776 | 0.8623260 | -0.1403464 | 0.0000083 | 0.0012960 | 0.8247206 | 0.0124536 | 0.0521815 | 0.0000372 | 0.0000007 | 0.4203287 | 1.581221 | 0.5244772 | 0.0000000 |
| 3941 | GO:0048245 eosinophil chemotaxis | 7 | 0.0000806 | 0.6525805 | 0.7160342 | 0.4465849 | 0.4398146 | 0.6447124 | 0.7443149 | 0.0059062 | -0.4226754 | 0.0027901 | 0.0010351 | 0.0407593 | 0.0439101 | 0.0031378 | 0.0006486 | 0.9784147 | 0.0528167 | 1.574735 | 0.4102053 | 0.0017405 |
| 2249 | GO:0030292 protein tyrosine kinase inhibitor activity | 5 | 0.0081547 | -0.8881537 | -0.6694130 | -0.1873924 | -0.2771540 | -0.8133618 | -0.5747279 | 0.1918321 | 0.0891782 | 0.0005826 | 0.0095352 | 0.4680958 | 0.2832044 | 0.0016335 | 0.0260468 | 0.4576200 | 0.7298731 | 1.544501 | 0.4073437 | 0.0663701 |
| 2044 | GO:0019774 proteasome core complex, beta-subunit complex | 10 | 0.0000000 | -0.2329680 | 0.0506192 | -0.6570249 | -0.8243274 | 0.1144626 | 0.2774662 | 0.7787046 | 0.7141068 | 0.2021366 | 0.7816803 | 0.0003208 | 0.0000063 | 0.5308799 | 0.1287294 | 0.0000200 | 0.0000920 | 1.540926 | 0.5816657 | 0.0000013 |
| 509 | GO:0004185 serine-type carboxypeptidase activity | 5 | 0.0006972 | -0.4305514 | -0.7292067 | -0.4098897 | -0.5764354 | -0.2667094 | -0.2756172 | 0.5881323 | 0.7727215 | 0.0954834 | 0.0047450 | 0.1124801 | 0.0256066 | 0.3017420 | 0.2858825 | 0.0227606 | 0.0027678 | 1.519046 | 0.5460508 | 0.0103500 |
| 3146 | GO:0038094 Fc-gamma receptor signaling pathway | 8 | 0.0007325 | -0.7920225 | -0.5865891 | -0.4400263 | -0.1786578 | -0.6797965 | -0.6783554 | 0.3632045 | 0.1428266 | 0.0001046 | 0.0040655 | 0.0311577 | 0.3816141 | 0.0008693 | 0.0008916 | 0.0752779 | 0.4842711 | 1.507161 | 0.4236976 | 0.0106825 |
| 494 | GO:0004045 aminoacyl-tRNA hydrolase activity | 5 | 0.0335599 | 0.4520384 | 0.5388403 | -0.4544575 | -0.5354536 | 0.6277481 | 0.6993810 | 0.5698328 | 0.2406119 | 0.0800552 | 0.0369312 | 0.0784530 | 0.0381343 | 0.0150628 | 0.0067624 | 0.0273451 | 0.3515166 | 1.501244 | 0.4901989 | 0.1679748 |
| 3289 | GO:0042719 mitochondrial intermembrane space protein transporter complex | 6 | 0.0092078 | 0.6126844 | 0.6680660 | -0.2104027 | 0.0061667 | 0.7944120 | 0.4642332 | 0.6373512 | -0.3402353 | 0.0093516 | 0.0045982 | 0.3721754 | 0.9791336 | 0.0007515 | 0.0489399 | 0.0068596 | 0.1489886 | 1.494852 | 0.4421493 | 0.0719076 |
| 4698 | GO:0071139 resolution of DNA recombination intermediates | 5 | 0.0693363 | -0.7335041 | -0.3362647 | -0.3366346 | 0.2305656 | -0.7298613 | -0.7233440 | 0.0275347 | -0.5904945 | 0.0045043 | 0.1929029 | 0.1924137 | 0.3719903 | 0.0047076 | 0.0050922 | 0.9150960 | 0.0222205 | 1.490957 | 0.3682664 | 0.2620653 |
| 2625 | GO:0032489 regulation of Cdc42 protein signal transduction | 5 | 0.0795511 | -0.7271291 | -0.5942796 | 0.1538385 | 0.1616364 | -0.7319957 | -0.7438919 | -0.3817432 | -0.0789043 | 0.0048655 | 0.0213782 | 0.5514032 | 0.5314116 | 0.0045875 | 0.0039675 | 0.1393682 | 0.7599747 | 1.474050 | 0.3947076 | 0.2809024 |
| 3664 | GO:0045569 TRAIL binding | 5 | 0.0000573 | 0.5554322 | 0.2128068 | 0.4731555 | 0.3337033 | 0.8467734 | 0.5629456 | 0.1925720 | 0.6379367 | 0.0314933 | 0.4099467 | 0.0669300 | 0.1963152 | 0.0010407 | 0.0292662 | 0.4558869 | 0.0134990 | 1.472075 | 0.2227802 | 0.0013057 |
| 848 | GO:0005942 phosphatidylinositol 3-kinase complex | 6 | 0.0145353 | -0.6493762 | -0.5779612 | 0.2303733 | 0.5976235 | -0.5741426 | -0.7075091 | 0.0205398 | -0.3969926 | 0.0058762 | 0.0142222 | 0.3285135 | 0.0112439 | 0.0148753 | 0.0026885 | 0.9305787 | 0.0922063 | 1.467715 | 0.4818199 | 0.0967782 |
| 82 | GO:0000413 protein peptidyl-prolyl isomerization | 5 | 0.0260507 | 0.1174671 | 0.5441622 | -0.6581430 | -0.8464034 | 0.2368837 | 0.2304233 | 0.7293205 | 0.0117111 | 0.6492327 | 0.0351058 | 0.0108159 | 0.0010460 | 0.3590297 | 0.3722856 | 0.0047385 | 0.9638324 | 1.449385 | 0.5456352 | 0.1437921 |
| 1728 | GO:0014909 smooth muscle cell migration | 5 | 0.0531647 | 0.6546709 | 0.4877268 | 0.0836855 | -0.5075916 | 0.7214657 | 0.5501672 | 0.3595162 | 0.4464603 | 0.0112401 | 0.0589490 | 0.7459178 | 0.0493559 | 0.0052082 | 0.0331393 | 0.1638991 | 0.0838513 | 1.443225 | 0.3974251 | 0.2251176 |
| 3343 | GO:0043009 chordate embryonic development | 5 | 0.0185867 | -0.4094913 | -0.4713340 | -0.3285806 | 0.5720811 | -0.6394166 | -0.6628958 | -0.4883529 | -0.3954891 | 0.1128299 | 0.0679874 | 0.2032721 | 0.0267420 | 0.0132842 | 0.0102583 | 0.0586246 | 0.1256759 | 1.438146 | 0.3912911 | 0.1149614 |
| 4557 | GO:0070106 interleukin-27-mediated signaling pathway | 8 | 0.0000000 | -0.4279640 | -0.4482102 | 0.1068887 | -0.3374075 | 0.6368488 | 0.0750783 | 0.8907095 | 0.5940436 | 0.0360851 | 0.0281535 | 0.6006594 | 0.0984499 | 0.0018129 | 0.7131189 | 0.0000128 | 0.0036194 | 1.437630 | 0.5234869 | 0.0000000 |
| 2338 | GO:0030836 positive regulation of actin filament depolymerization | 5 | 0.0365654 | 0.2951121 | 0.8367556 | -0.4430167 | -0.2144290 | 0.4142725 | 0.6591960 | 0.5519886 | -0.3447456 | 0.2531691 | 0.0011933 | 0.0862667 | 0.4063869 | 0.1086885 | 0.0106901 | 0.0325617 | 0.1819155 | 1.434990 | 0.4890421 | 0.1762624 |
| 2626 | GO:0032494 response to peptidoglycan | 5 | 0.0343920 | -0.6559516 | -0.1751832 | -0.6289719 | -0.0261402 | -0.6085663 | -0.6242191 | 0.2679616 | -0.5922874 | 0.0110820 | 0.4975775 | 0.0148669 | 0.9193810 | 0.0184445 | 0.0156406 | 0.2994778 | 0.0218180 | 1.428237 | 0.3549864 | 0.1706680 |
| 5439 | GO:1903238 positive regulation of leukocyte tethering or rolling | 6 | 0.0030707 | -0.7879370 | -0.6156254 | 0.0955125 | -0.0205398 | -0.6793084 | -0.6729045 | 0.1048100 | -0.3038281 | 0.0008300 | 0.0090171 | 0.6854010 | 0.9305787 | 0.0039564 | 0.0043112 | 0.6566534 | 0.1975087 | 1.423707 | 0.3761212 | 0.0327955 |
| 2514 | GO:0031726 CCR1 chemokine receptor binding | 6 | 0.0000051 | 0.5605284 | 0.5354585 | 0.3863194 | 0.4504056 | 0.5897491 | 0.4824961 | 0.3481571 | -0.6057113 | 0.0174238 | 0.0231301 | 0.1012971 | 0.0560746 | 0.0123628 | 0.0406982 | 0.1397517 | 0.0101892 | 1.421873 | 0.3924636 | 0.0001545 |
| 1244 | GO:0007221 positive regulation of transcription of Notch receptor target | 6 | 0.0102317 | -0.6458185 | -0.2970922 | 0.3827143 | 0.4489825 | -0.7309663 | -0.3191025 | -0.6441582 | -0.3332147 | 0.0061530 | 0.2076297 | 0.1045243 | 0.0568557 | 0.0019299 | 0.1759048 | 0.0062862 | 0.1575577 | 1.419687 | 0.4541522 | 0.0762537 |
| 3261 | GO:0042555 MCM complex | 11 | 0.0000149 | -0.5308886 | -0.1236775 | -0.0956328 | 0.1744090 | -0.6823585 | -0.6680450 | -0.4793288 | -0.7115162 | 0.0022979 | 0.4776183 | 0.5829314 | 0.3166139 | 0.0000889 | 0.0001247 | 0.0059127 | 0.0000438 | 1.408762 | 0.3316884 | 0.0004071 |
| 2127 | GO:0022624 proteasome accessory complex | 17 | 0.0000000 | 0.0684746 | 0.4489446 | -0.6710859 | -0.5239495 | 0.1609971 | 0.6610076 | 0.6076588 | 0.4492965 | 0.6250760 | 0.0013531 | 0.0000017 | 0.0001840 | 0.2505647 | 0.0000024 | 0.0000144 | 0.0013414 | 1.401816 | 0.5049305 | 0.0000002 |
| 4968 | GO:0097100 supercoiled DNA binding | 6 | 0.0103035 | -0.8066031 | -0.4078791 | -0.3344955 | -0.2658081 | -0.7853280 | -0.4626678 | 0.3191262 | 0.1759167 | 0.0006221 | 0.0836207 | 0.1559672 | 0.2595680 | 0.0008637 | 0.0497073 | 0.1758727 | 0.4555870 | 1.401109 | 0.4033688 | 0.0765901 |
| 3278 | GO:0042612 MHC class I protein complex | 9 | 0.0000000 | -0.1644406 | -0.5298555 | -0.5121897 | -0.3332226 | 0.3338236 | -0.1484671 | 0.8429689 | 0.6615188 | 0.3930383 | 0.0059150 | 0.0077987 | 0.0834768 | 0.0829213 | 0.4406179 | 0.0000119 | 0.0005889 | 1.401017 | 0.5291543 | 0.0000000 |
| 3102 | GO:0036150 phosphatidylserine acyl-chain remodeling | 7 | 0.0313505 | -0.6238322 | -0.6243405 | 0.0961970 | 0.1457848 | -0.6495715 | -0.7088776 | -0.3066859 | -0.2866800 | 0.0042604 | 0.0042292 | 0.6594426 | 0.5042296 | 0.0029187 | 0.0011622 | 0.1600255 | 0.1890720 | 1.382088 | 0.3415181 | 0.1617057 |
| 3109 | GO:0036261 7-methylguanosine cap hypermethylation | 8 | 0.0000693 | -0.2603188 | 0.5653644 | -0.5517364 | -0.5937055 | 0.0952000 | 0.4824936 | 0.7793375 | 0.0550811 | 0.2023598 | 0.0056220 | 0.0068864 | 0.0036386 | 0.6410649 | 0.0181238 | 0.0001348 | 0.7873596 | 1.377159 | 0.5148798 | 0.0015384 |
| 810 | GO:0005838 proteasome regulatory particle | 8 | 0.0005805 | 0.0096250 | 0.3962603 | -0.6439297 | -0.5458831 | 0.1316361 | 0.6481818 | 0.6052341 | 0.4508077 | 0.9624057 | 0.0522976 | 0.0016106 | 0.0075039 | 0.5191564 | 0.0014993 | 0.0030329 | 0.0272537 | 1.369946 | 0.4983447 | 0.0088246 |
| 2418 | GO:0031123 RNA 3’-end processing | 9 | 0.0016171 | -0.4063847 | -0.6972457 | 0.1600282 | 0.5485651 | -0.4131853 | -0.7029709 | -0.2334591 | -0.4014661 | 0.0347774 | 0.0002919 | 0.4058602 | 0.0043764 | 0.0318513 | 0.0002601 | 0.2252734 | 0.0370342 | 1.363232 | 0.4280643 | 0.0201996 |
| 657 | GO:0005216 monoatomic ion channel activity | 5 | 0.0019779 | 0.4727855 | 0.6879118 | -0.2717182 | -0.4412522 | 0.7074351 | 0.5489150 | 0.2332408 | -0.1273995 | 0.0671437 | 0.0077241 | 0.2927535 | 0.0875259 | 0.0061529 | 0.0335414 | 0.3664675 | 0.6218100 | 1.355613 | 0.4516970 | 0.0237763 |
| 2392 | GO:0031048 regulatory ncRNA-mediated heterochromatin formation | 7 | 0.0104014 | -0.6060424 | -0.2626892 | -0.2835084 | 0.3786380 | -0.6498155 | -0.5596670 | 0.0788342 | -0.6593305 | 0.0054918 | 0.2288148 | 0.1940131 | 0.0828039 | 0.0029081 | 0.0103429 | 0.7179929 | 0.0025202 | 1.355159 | 0.3807816 | 0.0772177 |
| 5425 | GO:1902975 mitotic DNA replication initiation | 5 | 0.0002895 | -0.5688367 | 0.2902170 | 0.2234223 | 0.2217147 | -0.6093348 | -0.4640768 | -0.4148702 | -0.7463963 | 0.0276161 | 0.2611271 | 0.3869886 | 0.3906277 | 0.0182972 | 0.0723379 | 0.1081794 | 0.0038470 | 1.350095 | 0.4289682 | 0.0049836 |
| 3157 | GO:0038156 interleukin-3-mediated signaling pathway | 6 | 0.0082110 | -0.8769034 | -0.3148333 | 0.2087662 | -0.0829420 | -0.8494616 | -0.4223234 | 0.0233860 | -0.0583701 | 0.0001991 | 0.1817576 | 0.3759061 | 0.7249977 | 0.0003137 | 0.0732408 | 0.9209887 | 0.8044681 | 1.349977 | 0.3997742 | 0.0664011 |
| 3896 | GO:0047631 ADP-ribose diphosphatase activity | 6 | 0.0150921 | 0.4170580 | 0.2280727 | -0.3243679 | -0.3430103 | 0.5778900 | 0.6872302 | 0.5193776 | 0.5427399 | 0.0768954 | 0.3333671 | 0.1688814 | 0.1457013 | 0.0142341 | 0.0035543 | 0.0275907 | 0.0213255 | 1.348847 | 0.4062484 | 0.0993346 |
| 2128 | GO:0022625 cytosolic large ribosomal subunit | 52 | 0.0000000 | -0.6126719 | -0.1827747 | -0.2278780 | -0.7080212 | -0.5047198 | -0.2091634 | 0.6332892 | 0.3806781 | 0.0000000 | 0.0226886 | 0.0044961 | 0.0000000 | 0.0000000 | 0.0091153 | 0.0000000 | 0.0000021 | 1.344046 | 0.4706219 | 0.0000000 |
| 909 | GO:0006198 cAMP catabolic process | 7 | 0.0028394 | -0.2878186 | -0.7189212 | 0.2813939 | 0.7486251 | -0.2706590 | -0.6401582 | -0.0327841 | -0.2674874 | 0.1873211 | 0.0009876 | 0.1973594 | 0.0006032 | 0.2150011 | 0.0033567 | 0.8806189 | 0.2204247 | 1.339773 | 0.4808669 | 0.0308419 |
| 5255 | GO:0140948 histone H3K9 monomethyltransferase activity | 5 | 0.0567051 | -0.6920669 | -0.7094273 | -0.1018143 | 0.0932480 | -0.6809107 | -0.5524155 | -0.1510210 | -0.0240911 | 0.0073624 | 0.0060100 | 0.6934179 | 0.7180604 | 0.0083697 | 0.0324276 | 0.5587163 | 0.9256816 | 1.339214 | 0.3381680 | 0.2343407 |
| 4219 | GO:0051383 kinetochore organization | 5 | 0.0235042 | -0.3923301 | -0.2951121 | -0.6318463 | 0.2255852 | -0.4446674 | -0.7004625 | 0.2920669 | -0.5721380 | 0.1287260 | 0.2531691 | 0.0144156 | 0.3824088 | 0.0851020 | 0.0066776 | 0.2580996 | 0.0267269 | 1.338714 | 0.3777838 | 0.1346514 |
| 5534 | GO:1904948 midbrain dopaminergic neuron differentiation | 5 | 0.0336407 | -0.6859481 | -0.3155745 | -0.4158093 | 0.3191604 | -0.5602704 | -0.6169050 | 0.3152046 | -0.3827677 | 0.0079006 | 0.2217374 | 0.1073832 | 0.2165303 | 0.0300432 | 0.0169004 | 0.2222797 | 0.1383096 | 1.336118 | 0.3962429 | 0.1680865 |
| 2926 | GO:0035005 1-phosphatidylinositol-4-phosphate 3-kinase activity | 6 | 0.0318359 | -0.0637304 | -0.7175181 | 0.2803235 | 0.6746122 | -0.2375599 | -0.6113325 | -0.4973910 | -0.2285470 | 0.7869266 | 0.0023365 | 0.2344466 | 0.0042139 | 0.3136490 | 0.0095090 | 0.0348773 | 0.3323626 | 1.335055 | 0.4685799 | 0.1621712 |
| 4989 | GO:0097250 mitochondrial respirasome assembly | 6 | 0.0078791 | 0.1798539 | 0.5307386 | -0.1162658 | -0.7729709 | 0.5131161 | 0.4716807 | 0.5783881 | 0.1405768 | 0.4455632 | 0.0243694 | 0.6219227 | 0.0010415 | 0.0295185 | 0.0454238 | 0.0141507 | 0.5510146 | 1.328570 | 0.4589351 | 0.0650476 |
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)
Session information
For reproducibility.
save.image("scanalysis.Rdata")
sessionInfo()## R version 4.4.2 (2024-10-31)
## 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.2 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.2 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.2 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.2
## [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.8
## [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