PART 1: Prepare data from PREDICT cohort

Let’s load the PREDICT Seurat object and the markers of all end clusters from tiered clustering (Supp. Table XY). Keep in mind that the .rds object is really huge and might consume as much as ~40 GB of your RAM. These are the 25 end clusters driving the anti-TNF response explanatory PC2. These clusters will be used for marker selection.

# seu <- readRDS(path_predict_seu)

markers_state <- as.data.frame(read_delim(paste0(path_data,"Supp_tables/CD14_03_cellstate_markerstop50sig.tsv"), "\t", escape_double = FALSE, trim_ws = TRUE))


top_pos_pc2 <-  paste0("CD.",c("T.MT-CO2.MT-CO1","EC.GNAT3.TRPM5","EC.RBP2.CYP3A4","cDC2.CLEC10A.FCGR2B","Myeloid.MKI67.IGKC","Goblet.FCGBP.HES6","EC.ADH1C.EDN1","Mcell.CSRP2.SPIB","EC.FABP6.PLCG2","EC.MTRNR2L1.MT-ND3","T.MT-CO2.CCR7","EC.ADH1C.RPS4Y1","T.IFI44L.PTGER4","cDC2.CD1C.NDRG2","Goblet.FCGBP.ITLN1","EC.FABP1.ADIRF","EC.GSTA2.TMPRSS15","cDC2.CLEC10A.CD1E","Mac.CXCL8.HES1","cDC2.CD1C.CD207","Mono/cDC2.CLEC10A.SDS","T.LAG3.BATF","Secretory.GSTA1.REG1B","EC.ADH1C.GSTA1","EpithStem.LINC00176.RPS4Y1"))

top_neg_pc2 <-  paste0("CD.",c(
  "Mono.FCN1.LYST","NK.MKI67.GZMA","Mono.CXCL10.TNF","T.CARD16.GB2","Mono.FCN1.S100A4","Mono/Mac.CXCL10.CXCL11","NK.GNLY.GZMB","EC.OLFM4.MT-ND2","NK.GNLY.IFNG","T.MKI67.IL22","NK.GNLY.FCER1G","Mac/DC.CXCL10.CLEC4E","T.MKI67.IFNG","T.MKI67.FOXP3","Mono/Mac.CXCL10.FCN1","T.GNLY.CSF2","EC.NUPR1.LCN2","Mac.CXCL3.APOC1","pDC.IRF7.IL3RA","Mono.CXCL3.FCN1","Mac.C1QB.CD14","Goblet.RETNLB.ITLN1","T.CCL20.RORA","T.MAF.CTLA4","cDC2.CD1C.AREG"))

kableExtra::kable(data.frame(PC2_top25_neg = top_neg_pc2[order(top_neg_pc2)],
                             PC2_top25_pos = top_pos_pc2[order(top_pos_pc2)])) %>%
  kable_styling(full_width = F, font_size = 10, bootstrap_options = c("hover", "condensed", "responsive")) %>%  column_spec(1:2, background = "white", color = "black")

We will select top 10 markers of each of these clusters. The selection is based on rank score function (see Methods). Duplicate hits are only counted once and genes that appear as hits for both positive and negative PC2 clusters are filtered out.

The final gene list for PC2_top25_neg contains 92 genes:

markers_top_neg_pcs <- markers_state %>% dplyr::filter(cluster %in% c(top_neg_pc2))
markers_top_pos_pcs <- markers_state %>% dplyr::filter(cluster %in% c(top_pos_pc2)) 

markers_top_neg_pcs <- rank_score_func(markers_top_neg_pcs)
markers_top_pos_pcs <- rank_score_func(markers_top_pos_pcs)

markers_top_neg_pcs <-  markers_top_neg_pcs %>% group_by(cluster) %>% slice_max(order_by = score, n = 10)
markers_top_pos_pcs <- markers_top_pos_pcs %>% group_by(cluster) %>% slice_max(order_by = score, n = 10)

markers_top_neg_pcs_genes <- unique(markers_top_neg_pcs$gene)
markers_top_pos_pcs_genes <- unique(markers_top_pos_pcs$gene)

markers_top_neg_pcs_genes_filt <- markers_top_neg_pcs_genes[which(markers_top_neg_pcs_genes %in% markers_top_pos_pcs_genes == F)]
markers_top_pos_pcs_genes_filt <- markers_top_pos_pcs_genes[which(markers_top_pos_pcs_genes %in% markers_top_neg_pcs_genes == F)]

knitr::kable(matrix(c(markers_top_neg_pcs_genes_filt[order(markers_top_neg_pcs_genes_filt)], rep("",4)), ncol = 8)) %>%
  kable_styling(full_width = F, font_size = 10,bootstrap_options = c("hover", "condensed", "responsive")) %>%  column_spec(1:8, background = "white", color = "black")

The final gene list for PC2_top25_pos contains 103 genes:

knitr::kable(matrix(c(markers_top_pos_pcs_genes_filt[order(markers_top_pos_pcs_genes_filt)], ""), ncol = 8)) %>%
  kable_styling(full_width = F, font_size = 10,bootstrap_options = c("hover", "condensed", "responsive")) %>%  column_spec(1:8, background = "white", color = "black")

PART 2: RISK cohort

To load data from RISK cohort, please download the bulk_data.rd file from this github repository. The path to the file in the repository is martin_et_al_cell_2019/input/external_cohorts_data/bulk_data.rd. For convenience, we also provide this file together with our code. We assume that this file is stored in the /data/RISK/ folder.

martin_bulk_data <- load(paste0(path_data,"RISK/bulk_data.rd"))

# Metadata are stored in the slot design
bulk_metadata <- bulk$design

# Count matrix is stored in the slot raw
bulk_matrix <- bulk$raw

# Filtering - only genes with reported expression and only RISK patients with anti-TNF response info
bulk_risk_response_names <- rownames(bulk_metadata %>% dplyr::filter(is.na(response)==F))
bulk_matrix_filt <- bulk_matrix[which(is.na(rowSums(bulk_matrix))==F),which(colnames(bulk_matrix) %in% bulk_risk_response_names)]


response <- bulk_metadata$response
names(response) <- rownames(bulk_metadata)

response_filt <- response[which(names(response) %in% colnames(bulk_matrix_filt))]
response_filt <- response_filt[match(names(response_filt),colnames(bulk_matrix_filt))]

# Create Seurat object and add response metadata
RISK_seurat <- CreateSeuratObject(bulk_matrix_filt)
RISK_seurat <- AddMetaData(RISK_seurat, response_filt, col.name = "response")

# Calculate fold changes for all genes in responders vs. non-responders
Idents(RISK_seurat) <- RISK_seurat$response
fc_risk_seurat <- FoldChange(RISK_seurat, ident.1 = "R", ident.2 = "NR")

# The calculated fold changes will be used as input for GSEA
fc_risk_seurat.input <- as.numeric(fc_risk_seurat$avg_log2FC)
names(fc_risk_seurat.input) <- rownames(fc_risk_seurat)
fc_risk_seurat.input <- fc_risk_seurat.input[which(!is.na(fc_risk_seurat.input))]

# Prerank genes with similar fold change
set.seed(123)
fc_risk_seurat.input <- fc_risk_seurat.input[rank(fc_risk_seurat.input,  ties.method = "random")]

fgseaRes_risk_neg <- fgsea(pathways = list(our_list = markers_top_neg_pcs_genes_filt), 
                           stats = fc_risk_seurat.input)

p1 <- plotEnrichment(markers_top_neg_pcs_genes_filt, fc_risk_seurat.input) + labs(title="\n\nTOP neg PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_risk_neg$pval, digits = 2)))


fgseaRes_risk_pos <- fgsea(pathways = list(our_list = markers_top_pos_pcs_genes_filt), 
                           stats    = fc_risk_seurat.input)

p2 <- plotEnrichment(markers_top_pos_pcs_genes_filt, fc_risk_seurat.input) + labs(title="\n\nTOP pos PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_risk_pos$pval, digits = 2)))

cowplot::plot_grid(p1,p2)

PART 3: DATASET E_MTAB_7604

The second dataset can be obtained from Array-Express under the accession number E-MTAB-7604. It should be unpacked and stored in the folder data/E_MTAB_7604/.

emtab_7604 <- c('GC058859', 'GC058865', 'GC058866', 'GC058872', 'GC058874', 'GC058911', 'GC050149', 'GC050152',
                'GC050158', 'GC050162', 'GC050173', 'GC050174', 'GC050260', 'GC050265', 'GC050272', 'GC050281', 'GC050282', 'GC050312', 'GC050318', 'GC050320', 'GC050321', 'GC050322', 'GC058827', 'GC058833', 'GC058835', 'GC058840', 'GC058849', 'GC058851', 'GC058889', 'GC058897', 'GC058904', 'GC058910', 'GC049971', 'GC049992', 'GC049997', 'GC049998', 'GC050011', 'GC050016', 'GC050024', 'GC050040', 'GC050130', 'GC050131', 'GC050193')

# Load raw counts from all patients and construct a dataframe
df <- data.frame(gene = "")
for(j in emtab_7604){
  k <-  read_tsv(paste0(path_data,"E_MTAB_7604/",j,".count"), col_names = FALSE, 
                 trim_ws = TRUE)
  df_temp <- data.frame(gene = k$X1, value = k$X2)
  colnames(df_temp)[2] <- j
  df <- right_join(df, df_temp, by = "gene")
}

# Convert dataframe to matrix
mtx <- as.matrix(df %>% dplyr::select(-gene))
rownames(mtx) <- df$gene

# Load metadata
E_MTAB_7604_sdrf <- read_delim(paste0(path_data,"E_MTAB_7604/E-MTAB-7604.sdrf.txt"), 
                               "\t", escape_double = FALSE, trim_ws = TRUE)

md_emtab_7604 <- as.data.frame(E_MTAB_7604_sdrf)
md_emtab_7604$patient <- substr(md_emtab_7604$`Source Name`,1,8)
md_emtab_7604 <- data.frame(patient = md_emtab_7604$patient,
                            treatment = md_emtab_7604$`Factor Value[treatment]`,
                            response = md_emtab_7604$`Characteristics[clinical history]`,
                            site = md_emtab_7604$`Characteristics[organism part]`)

# Counts of patients in the dataset
knitr::kable(md_emtab_7604 %>% group_by(treatment, response, site) %>% tally()) %>%
  kable_styling(full_width = F, font_size = 10,bootstrap_options = c("hover", "condensed", "responsive")) %>%  column_spec(1:4, background = "white", color = "black")

# Prepare DeSeq input
patients.ileum <- (md_emtab_7604 %>% dplyr::filter(site == "ileum"))$patient
mtx.ileum <- mtx[,which(colnames(mtx) %in% patients.ileum)]
featureData <- data.frame(gene=rownames(mtx.ileum))

# Response metadata
ileum.response <- c(md_emtab_7604 %>% dplyr::filter(patient %in% colnames(mtx.ileum)))$response
names(ileum.response) <- c(md_emtab_7604 %>% dplyr::filter(patient %in% colnames(mtx.ileum)))$patient
ileum.response <- ileum.response[match(colnames(mtx.ileum),names(ileum.response))]

coldata <- data.frame(sample_type = factor(if_else(ileum.response == "responder","R","NR"),
                                           levels = c("R","NR")))

# Create a DESeq dataset

dds.ileum <- DESeq2::DESeqDataSetFromMatrix(countData = mtx.ileum,
                                            colData = coldata,
                                            design= ~ sample_type) 

rownames(dds.ileum) <- featureData$gene

# Run DESeq2
dds.ileum <- DESeq(dds.ileum)
dds.ileum

## class: DESeqDataSet 
## dim: 56643 17 
## metadata(1): version
## assays(6): counts mu ... replaceCounts replaceCooks
## rownames(56643): 5S_rRNA 7SK ... __not_aligned __alignment_not_unique
## rowData names(23): baseMean baseVar ... maxCooks replace
## colnames(17): GC058865 GC058866 ... GC050024 GC050040
## colData names(3): sample_type sizeFactor replaceable

# Retrieve the results
res.ileum <- results(dds.ileum)
res.ileum_df <- data.frame(gene = rownames(res.ileum),
                           pval = res.ileum@listData$pvalue,
                           logfc = res.ileum@listData$log2FoldChange,
                           padj = res.ileum@listData$padj)

# Prepare GSEA input
fc_emtab7604_ileum.input <- as.numeric(res.ileum_df$logfc)
names(fc_emtab7604_ileum.input) <- res.ileum_df$gene
fc_emtab7604_ileum.input <- fc_emtab7604_ileum.input[which(!is.na(fc_emtab7604_ileum.input))]

# Prerank genes with similar FC
set.seed(123)
fc_emtab7604_ileum.input <- fc_emtab7604_ileum.input[rank(fc_emtab7604_ileum.input,  ties.method = "random")]

# Run GSEA Top PC2 neg markers
fgseaRes_emtab7604_ileum_neg <- fgsea(pathways = list(our_list = markers_top_neg_pcs_genes_filt), 
                                      stats    = fc_emtab7604_ileum.input)


p3 <- plotEnrichment(markers_top_neg_pcs_genes_filt, fc_emtab7604_ileum.input) + labs(title="\n\nTOP neg PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_emtab7604_ileum_neg$pval, digits = 2)))

# Run GSEA Top PC2 pos markers

fgseaRes_emtab7604_ileum_pos <- fgsea(pathways = list(our_list = markers_top_pos_pcs_genes_filt), 
                                      stats    = fc_emtab7604_ileum.input)

p4 <- plotEnrichment(markers_top_pos_pcs_genes_filt, fc_emtab7604_ileum.input) + labs(title="\n\nTOP pos PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_emtab7604_ileum_pos$pval, digits = 2)))

cowplot::plot_grid(p3,p4)

# Prepare DeSeq input
patients.colon <- (md_emtab_7604 %>% dplyr::filter(site == "colon"))$patient
mtx.colon <- mtx[,which(colnames(mtx) %in% patients.colon)]
featureData <- data.frame(gene=rownames(mtx.colon))

# Response metadata
colon.response <- c(md_emtab_7604 %>% dplyr::filter(patient %in% colnames(mtx.colon)))$response
names(colon.response) <- c(md_emtab_7604 %>% dplyr::filter(patient %in% colnames(mtx.colon)))$patient
colon.response <- colon.response[match(colnames(mtx.colon),names(colon.response))]

coldata <- data.frame(sample_type = factor(if_else(colon.response == "responder","R","NR"),
                                           levels = c("R","NR")))

# Create a DESeq dataset

dds.colon <- DESeq2::DESeqDataSetFromMatrix(countData = mtx.colon,
                                            colData = coldata,
                                            design= ~ sample_type) 

rownames(dds.colon) <- featureData$gene

# Run DESeq2
dds.colon <- DESeq(dds.colon)
dds.colon

## class: DESeqDataSet 
## dim: 56643 26 
## metadata(1): version
## assays(6): counts mu ... replaceCounts replaceCooks
## rownames(56643): 5S_rRNA 7SK ... __not_aligned __alignment_not_unique
## rowData names(23): baseMean baseVar ... maxCooks replace
## colnames(26): GC058859 GC058874 ... GC050131 GC050193
## colData names(3): sample_type sizeFactor replaceable

# Retrieve the results
res.colon <- results(dds.colon)
res.colon_df <- data.frame(gene = rownames(res.colon),
                           pval = res.colon@listData$pvalue,
                           logfc = res.colon@listData$log2FoldChange,
                           padj = res.colon@listData$padj)

# Prepare GSEA input
fc_emtab7604_colon.input <- as.numeric(res.colon_df$logfc)
names(fc_emtab7604_colon.input) <- res.colon_df$gene
fc_emtab7604_colon.input <- fc_emtab7604_colon.input[which(!is.na(fc_emtab7604_colon.input))]

# Prerank genes with similar FC
set.seed(123)
fc_emtab7604_colon.input <- fc_emtab7604_colon.input[rank(fc_emtab7604_colon.input,  ties.method = "random")]


# Run GSEA Top PC2 neg markers
fgseaRes_emtab7604_colon_neg <- fgsea(pathways = list(our_list = markers_top_neg_pcs_genes_filt), 
                                      stats    = fc_emtab7604_colon.input)


p5 <- plotEnrichment(markers_top_neg_pcs_genes_filt, fc_emtab7604_colon.input) + labs(title="\n\nTOP neg PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_emtab7604_colon_neg$pval, digits = 2)))

# Run GSEA Top PC2 pos markers

fgseaRes_emtab7604_colon_pos <- fgsea(pathways = list(our_list = markers_top_pos_pcs_genes_filt), 
                                      stats    = fc_emtab7604_colon.input)

p6 <- plotEnrichment(markers_top_pos_pcs_genes_filt, fc_emtab7604_colon.input) + labs(title="\n\nTOP pos PC2 clusters", subtitle = paste0("p-value = ",format.pval(fgseaRes_emtab7604_colon_pos$pval, digits = 2)))

cowplot::plot_grid(p5,p6)

Save the plots

p_final <- cowplot::plot_grid(
  cowplot::plot_grid(p1, p2, ncol = 2),
  cowplot::plot_grid(p3, p4, ncol = 2),
  cowplot::plot_grid(p5, p6, ncol = 2),
  ncol = 1,
  labels = c("RISK","E-MTAB-7604 ileum","E-MTAB-7604 colon"))

dir.create(paste0("results/"))

ggsave(p_final,
       filename = paste0("results/GSEA_cohorts.png"),
       device = "png",
       width = 9,
       height = 8)

ggsave(p_final,
       filename = paste0("results/GSEA_cohorts.eps"),
       device = "eps",
       width = 9,
       height = 8)

SessionInfo

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Fedora 32 (Server Edition)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.12.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] DESeq2_1.30.1               SummarizedExperiment_1.20.0
##  [3] Biobase_2.50.0              MatrixGenerics_1.2.1       
##  [5] matrixStats_0.58.0          GenomicRanges_1.42.0       
##  [7] GenomeInfoDb_1.26.2         IRanges_2.24.1             
##  [9] S4Vectors_0.28.1            BiocGenerics_0.36.0        
## [11] kableExtra_1.3.4            fgsea_1.16.0               
## [13] ggVennDiagram_1.1.1         forcats_0.5.1              
## [15] stringr_1.4.0               dplyr_1.0.6                
## [17] purrr_0.3.4                 readr_1.4.0                
## [19] tidyr_1.1.3                 tibble_3.1.2               
## [21] ggplot2_3.3.3               tidyverse_1.3.0            
## [23] cowplot_1.1.1               DT_0.17                    
## [25] SeuratObject_4.0.2          Seurat_4.0.3               
## [27] Matrix_1.3-4               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.1             reticulate_1.18        tidyselect_1.1.0      
##   [4] RSQLite_2.2.3          AnnotationDbi_1.52.0   htmlwidgets_1.5.3     
##   [7] grid_4.0.3             BiocParallel_1.24.1    Rtsne_0.15            
##  [10] munsell_0.5.0          codetools_0.2-18       ica_1.0-2             
##  [13] future_1.21.0          miniUI_0.1.1.1         withr_2.4.1           
##  [16] colorspace_2.0-0       highr_0.8              knitr_1.31            
##  [19] rstudioapi_0.13        ROCR_1.0-11            tensor_1.5            
##  [22] listenv_0.8.0          labeling_0.4.2         GenomeInfoDbData_1.2.4
##  [25] polyclip_1.10-0        farver_2.1.0           bit64_4.0.5           
##  [28] parallelly_1.26.0      vctrs_0.3.8            generics_0.1.0        
##  [31] xfun_0.23              R6_2.5.0               locfit_1.5-9.4        
##  [34] RVenn_1.1.0            bitops_1.0-6           spatstat.utils_2.2-0  
##  [37] cachem_1.0.5           DelayedArray_0.16.2    assertthat_0.2.1      
##  [40] promises_1.2.0.1       scales_1.1.1           gtable_0.3.0          
##  [43] globals_0.14.0         goftest_1.2-2          rlang_0.4.10          
##  [46] genefilter_1.72.1      systemfonts_1.0.1      splines_4.0.3         
##  [49] lazyeval_0.2.2         spatstat.geom_2.2-0    broom_0.7.5           
##  [52] yaml_2.2.1             reshape2_1.4.4         abind_1.4-5           
##  [55] modelr_0.1.8           backports_1.2.1        httpuv_1.5.5          
##  [58] tools_4.0.3            ellipsis_0.3.2         spatstat.core_2.1-2   
##  [61] jquerylib_0.1.4        RColorBrewer_1.1-2     ggridges_0.5.3        
##  [64] Rcpp_1.0.6             plyr_1.8.6             zlibbioc_1.36.0       
##  [67] RCurl_1.98-1.2         rpart_4.1-15           deldir_0.2-10         
##  [70] pbapply_1.4-3          zoo_1.8-9              haven_2.3.1           
##  [73] ggrepel_0.9.1          cluster_2.1.1          fs_1.5.0              
##  [76] magrittr_2.0.1         data.table_1.14.0      scattermore_0.7       
##  [79] lmtest_0.9-38          reprex_1.0.0           RANN_2.6.1            
##  [82] fitdistrplus_1.1-3     hms_1.0.0              patchwork_1.1.1       
##  [85] mime_0.10              evaluate_0.14          xtable_1.8-4          
##  [88] XML_3.99-0.5           readxl_1.3.1           gridExtra_2.3         
##  [91] compiler_4.0.3         KernSmooth_2.23-18     crayon_1.4.1          
##  [94] htmltools_0.5.1.1      mgcv_1.8-34            later_1.1.0.1         
##  [97] geneplotter_1.68.0     lubridate_1.7.10       DBI_1.1.1             
## [100] dbplyr_2.1.0           MASS_7.3-53.1          cli_2.5.0             
## [103] igraph_1.2.6           pkgconfig_2.0.3        plotly_4.9.3          
## [106] spatstat.sparse_2.0-0  xml2_1.3.2             svglite_2.0.0         
## [109] annotate_1.68.0        bslib_0.2.5.1          webshot_0.5.2         
## [112] XVector_0.30.0         rvest_0.3.6            digest_0.6.27         
## [115] sctransform_0.3.2      RcppAnnoy_0.0.18       spatstat.data_2.1-0   
## [118] rmarkdown_2.7          cellranger_1.1.0       leiden_0.3.7          
## [121] fastmatch_1.1-0        uwot_0.1.10            shiny_1.6.0           
## [124] lifecycle_1.0.0        nlme_3.1-152           jsonlite_1.7.2        
## [127] viridisLite_0.3.0      fansi_0.4.2            pillar_1.6.1          
## [130] lattice_0.20-41        fastmap_1.1.0          httr_1.4.2            
## [133] survival_3.2-7         glue_1.4.2             png_0.1-7             
## [136] bit_4.0.4              stringi_1.5.3          sass_0.4.0            
## [139] blob_1.2.1             memoise_2.0.0          irlba_2.3.3           
## [142] future.apply_1.7.0