Introduction

The processed scRNA data dowloaded from GEO were alwalys like this.

Files of multiple samples were deposited together. But actually what we want is like follows:

Code exmaple

cd data/GSE184362  # Switch to data directory
ls *matrix.mtx.gz | awk '{split($1,a, "_matrix");printf "mkdir "a[1]"; ""mv "a[1]"* ./"a[1]"/\n"}' > script.sh
ls | grep GSM | awk '{printf "mv "$1"/"$1"_barcodes.tsv.gz "$1"/barcodes.tsv.gz; mv "$1"/"$1"_features.tsv.gz "$1"/features.tsv.gz; mv "$1"/"$1"_matrix.mtx.gz "$1"/matrix.mtx.gz\n"}' > script2.sh
bash script.sh
bash script2.sh

References

None

Original

None

Introduction

I found process was shut down due to low memory when I merged multiple seurat objects from different patients. Here is my solution. Firstly, I saved all expression matrix (in dgCMatrix format) from 10X directory into one list file. You can filter cells this time if you had correspondingly metadata. Secondly we got union genes from all matrix, then we added 0 to each matrix so that each matrix with expression of union genes. Thirly, we merge expression of all samples by using cbind instead of merge function. Lastly, we used this matrix to create seurat object.

Code exmaple

library(Seurat)
library(data.table)
library(ggplot2)
library(tidyverse)
library(magrittr)
library(Matrix)

setwd("/home/hxzk/project/sc_gbm/")

obj.ls <- list()

metadata <- fread("data/Metadata/Meta_Data_GBMatlas.txt")
setnames(metadata, "V1", "NAME")

# Read 10X matrix of all samples
for(f in dir("data/Matrix/", pattern = "^GSM")){
   Sample <- f
   P <- strsplit(Sample, split = "_")[[1]][1]
   tmp_mt <- Read10X(paste0("data/Matrix/", f))
   colnames(tmp_mt) <- paste0(P, "_", colnames(tmp_mt))
   tmp_mt <- tmp_mt[,colnames(tmp_mt) %in% metadata$NAME]
   obj.ls[[f]] <- tmp_mt
}

# Get union genes
get_all_genes <- function(obj.ls){
   all_genes <- rownames(obj.ls[[names(obj.ls)[1]]])
   for(f in names(obj.ls)[2:length(names(obj.ls))]){
      all_genes <- union(all_genes, rownames(obj.ls[[f]]))
   }
   return(all_genes)
}

# add 0
add_zero <- function(mt, all_genes){
   gene_left <- setdiff(all_genes, rownames(mt))
   left_mt <- as(matrix(0, ncol = ncol(mt), nrow = length(gene_left)), "dgCMatrix")
   colnames(left_mt)  <- colnames(mt)
   rownames(left_mt) <- gene_left
   mt <- rbind(mt, left_mt)
   return(mt[all_genes,])
}

# Cbind matrix 
cbind_dgC_lst <- function(dgc_lst){
    merge_mt <- dgc_lst[[names(dgc_lst)[1]]]
    for(i in names(dgc_lst)[2:length(names(dgc_lst))]){
       merge_mt <- cbind(merge_mt, dgc_lst[[i]])
   }
   return(merge_mt)
}

all_genes <- get_all_genes(obj.ls)
new_obj.ls <- lapply(obj.ls, add_zero, all_genes)
rm(obj.ls)
gc()
all_mt <- cbind_dgC_lst(new_obj.ls)
rm(new_obj.ls)
gc()

# Create seurat object
combined <- CreateSeuratObject(counts = all_mt, min.cells = 3)
# You can also use metadata you had
metadata <- metadata[rownames(combined@meta.data),, on = "NAME"]

References

None

Original

None

Introduction

Several data are shared in NCBI website, so knowning the methods for downloading data is necessary, especialy for people working on data mining.

Code exmaple

1. Perfetch data
   After we downloaded SRR_Acc_List.txt from ncbi data, we can use prefetch to download sra file which were listed in SRR_Acc_List.txt.
   1.1 Dowload accession list file, install sratoolkit from ncbi website
   1.2 Use prefetch in sratoolkit to download all sra file in

   ./prefetch --option-file  /yjdata/data/Candidemia/SRR_Acc_List.txt -O /yjdata/data/Candidemia/sra/
   # Continue download
   ./prefetch --option-file  /yjdata/data/Candidemia/SRR_Acc_List.txt -O /yjdata/data/Candidemia/sra/ -r yes

2. fasterq dump
   In the past， we uncompressed sra to fastq by fastq-dump. But now fastq-dump is replaced by fasterq-dump, which are multi-threaded operation and more faster.

   # Prepare codes
   less Healthy_RunTable.txt | awk 'BEGIN{FS = ",";OFS = "\t"}{print "echo "$1";""/yjdata/hxzk/utilities/sratoolkit.3.0.0-centos_linux64/bin/fasterq-dump --split-3 -O /yjdata/data/Candidemia/fq/ /yjdata/data/Candidemia/sra/"$1"/"$1".sra"}'

3. STAR
   STAR is used to map fastq data to reference genome and GTF

   # Build STAR reference index
   STAR --runThreadN 20 --runMode genomeGenerate \
   --genomeDir /home/hxzk/download/star/GRCh38_Ensembl99/ \
   --genomeFastaFiles /home/hxzk/download/star/GRCh38_Ensembl99/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
   --sjdbGTFfile /home/hxzk/download/star/GRCh38_Ensembl99/Homo_sapiens.GRCh38.99.gtf \
   --sjdbOverhang 100

   # star.sh
   hg=$1
   fq1=$2
   fq2=$3
   outprefix=$4
   mkdir $outFileNamePrefix
   STAR --runThreadN 20 --genomeDir $hg \
   --runMode alignReads \
   --quantMode TranscriptomeSAM GeneCounts \
   --readFilesIn  $fq1 $fq2 \
   --outSAMtype BAM SortedByCoordinate \
   --outFileNamePrefix $outprefix"/"

   # Prepared code
   less sample_table/Candidemia_Healthy_RunTable.txt | awk 'BEGIN{FS = ","}{print "bash star.sh /home/hxzk/download/star/GRCh38_Ensembl99/ /yjdata/data/Candidemia/fq/"$1"_1.fastq /yjdata/data/Candidemia/fq/"$1"_2.fastq ~/project/Candidemia/result/bam/"$1"/"}'

References

The manual of STAR
NCBIwebsite: https://www.ncbi.nlm.nih.gov/

Introduction

The following are the codes for classifying the samples according to the distance from the centroids of reference class.

Code exmaple

library(data.table)
centroids <- read.table("data/other/GSE10886/pam50_centroids.txt")
gene_name <- fread("data/other/GSE10886/pam_gene_name.txt", header = T)
rownames(centroids) <- gene_name[data.table(pcrID = rownames(centroids)), GeneName, on = "pcrID"]  # Get gene name

# Get intersect genes between my data and centoid data
inter_gene <- intersect(rownames(centroids), rownames(mt))
dat <- mt[inter_gene,]
centroids <- centroids[inter_gene,]

# Calculaed dist
d1=apply(centroids,2,function(i){
  dist(t(cbind(i,dat)))[1:ncol(dat)]
})
# Select the classe with min distance
p1=apply(d1,1, function(x){
  colnames(d1)[which.min(x)]
})

 p1  # Predicted class

References

None

Original

None

Donation

Introduction

The following are the codes for GO and KEGG enrichment.

Code exmaple

library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(magrittr)
library(data.table)
library(org.Hs.eg.db)

setwd("/home/hxzk/project/XXX/KEGG/")

up_dt <- fread("data/up_table.txt")
down_dt <- fread("data/down_table.txt")

# Transform symbol to entrez id, here "gene name" is the column of gene symbol
xx <- unlist(as.list(org.Hs.egALIAS2EG))
map_dt <- data.table(Symbol = names(xx), ID = xx)
map_dt <- map_dt[!duplicated(Symbol),]
up_dt$id <- map_dt[up_dt$gene name,ID, on = "Symbol"]
down_dt$id <- map_dt[down_dt$gene name,ID, on = "Symbol"]

# Remove data without gene id
gene.ls <- list(Up = up_dt[!is.na(id),]$id, Down = down_dt[!is.na(id)]$id)

# Run KEGG and GO enrichment
compKEGG <- compareCluster(geneCluster   = gene.ls,
                             fun           = "enrichKEGG",
                             pvalueCutoff  = 0.05,
                              pAdjustMethod = "BH",
                            organism = "hsa") # or organism = "mouse" for mouse 

compGO <- compareCluster(geneCluster   = gene.ls,
                           fun           = "enrichGO",
                           pvalueCutoff  = 0.05,
                           pAdjustMethod = "BH",
                   OrgDb = org.Hs.eg.db,
                   ont = 'BP') # or OrgDb = org.Mm.eg.db for mouse

up_eKEGG <- enrichKEGG(
    gene = gene.ls[["Up"]],
    pvalueCutoff = 0.05,
    organism = 'hsa'
  )
down_eKEGG <- enrichKEGG(
    gene = gene.ls[["Down"]],
    pvalueCutoff = 0.05,
    organism = 'hsa'
  )
up_eGO <- enrichGO(
    gene          = gene.ls[["Up"]],
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    readable      = TRUE
  )
down_eGO <- enrichGO(
    gene          = gene.ls[["Down"]],
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.05,
    readable      = TRUE
  )

get_GeneSymbol <- function(x){
   x <- strsplit(x, split = "/", fixed = T)[[1]]
   symbol <- map_dt[ID %in% x, Symbol]
   return(paste(symbol, collapse = "/"))
}

# add Gene symbo
up_eKEGG@result$Gene_symbol <- unlist(lapply(up_eKEGG@result$geneID, get_GeneSymbol))
down_eKEGG@result$Gene_symbol <- unlist(lapply(down_eKEGG@result$geneID, get_GeneSymbol))

fwrite(up_eKEGG@result, file = "result/up_KEGG.txt", sep = "\t", col.names = T)
fwrite(down_eKEGG@result, file = "result/down_KEGG.txt", sep = "\t", col.names = T)
fwrite(up_eGO@result, file = "result/up_GO.txt", sep = "\t", col.names = T)
fwrite(down_eGO@result, file = "result/down_GO.txt", sep = "\t", col.names = T)

Barplot codes for GO enrichment

library(data.table)
library(ggplot2)

setwd("/home/hxzk/project/Yuanyuan/KEGG/")

up_eGO_dt <- fread(file = "result/up_GO.txt")
down_eGO_dt <- fread(file = "result/down_GO.txt")

up_eGO_dt$Type <- "Up"
down_eGO_dt$Type <- "Down"

# Show top 10 up-regulated and down-regulated genes
dp <- rbind(up_eGO_dt[1:10,.(Term= Description, Pvalue=pvalue, Type)],
           down_eGO_dt[1:10,.(Term=Description, Pvalue=pvalue, Type)])
dp$Type <- factor(dp$Type, levels = c("Up", "Down"))
dp <- dp[nrow(dp):1,]  # Up first
dp$Term <- factor(dp$Term, levels = dp$Term)

p <- ggplot(dp, aes(x = Term, y = -log10(Pvalue))) +
 geom_bar(aes(fill = Type), stat = "identity") + theme_bw()+
   labs(x="Term", y=expression(-Log[10]*" (p value)"))+ coord_flip()

ggsave(p, filename = "result/GO_enrichment.pdf", width = 7, height = 3.3)

References

None

Original

None

Donation

Introduction

The following are the codes for drawing a volcano plot with ggplot2

Code exmaple

#Load package and read file
library(data.table)
library(ggplot2)
library(ggrepel)

setwd("/home/hxzk/project/XXX/KEGG/")

dt_pvalue <- fread("data/P_value_table.txt")

#Filter data
dt_pvalue <- dt_pvalue[!is.na(Pvalue),][, 1:12][!is.na(log2FC)][!is.infinite(log2FC)]

# Mark sinificant, positive, negagive or gene name to shown
dt_pvalue$changed <- factor(ifelse(dt_pvalue$Pvalue < 0.05 & abs(dt_pvalue$log2FC) > 0.6,
                              ifelse(dt_pvalue$log2FC > 0.6,'positive','negative'),'no'))
#dt_pvalue$selectedgene <- ifelse(dt_pvalue$gene name %in% c("BCL2L1", "TAB2", "WDR6"),dt_pvalue$gene name,NA)
ggplot(dt_pvalue,aes(log2FC,-log10(Pvalue),
                color=factor(changed),
                size=factor(changed)))+
  geom_point()+
  labs(x=expression(Log[2]*" Fold Change"),
       y=expression(-Log[10]*" (p value)"))+
  theme_classic(base_size = 15)+
  scale_color_manual(values = c('deepskyblue','grey','tomato'))+
  scale_size_manual(values = c(2,1,2))+
  theme(legend.title = element_blank(),
        legend.position = c(0.2,0.9),
        legend.background = element_rect(fill='transparent'))+
  geom_text_repel(aes(label=selectedgene), color="black",size=4,
                  box.padding=unit(0.5, "lines"),
                  point.padding=NA,
                  segment.colour = "black")+
  geom_hline(yintercept = -log10(0.05),linetype=2,cex=1)+  #add lines
  geom_vline(xintercept = c(-0.6,0.6),linetype=2,cex=1)+
  theme(axis.text= element_text(colour = "black"),
        panel.border = element_rect(size=1,fill='transparent'))
ggsave("result/volcano.pdf", width = 5, height = 4)

References

None

Original

None

Donation

Reference 
1.  Zhu, W.,  et  al.,  Pattern  recognition  receptor-initiated  innate  antiviral  responses  in  mouse 
epididymal epithelial cells. J Immunol, 2015. 194(10): p. 4825-35. 
2.  Carvajal, G., et al., Impaired male fertility and abnormal epididymal epithelium differentiation 
in mice lacking CRISP1 and CRISP4. Sci Rep, 2018. 8(1): p. 17531. 
3.  Krapf, D., et al., CSrc is necessary for epididymal development and is incorporated into sperm 
during epididymal transit. Developmental Biology, 2012. 369(1): p. 43-53. 
4.  Shum, W.W., et al., Regulation of luminal acidification in the male reproductive tract via cell-
cell crosstalk. J Exp Biol, 2009. 212(Pt 11): p. 1753-61. 
5.  Pietrement, C., et al., Role of NHERF1, Cystic Fibrosis Transmembrane Conductance Regulator, and cAMP in the Regulation of Aquaporin 9. Journal of Biological Chemistry, 2008. 283(5): p. 
2986-2996. 
6.  Shum, W.W., et al., Epithelial basal cells are distinct from dendritic cells and macrophages in 
the mouse epididymis. Biol Reprod, 2014. 90(5): p. 90. 
7.  Shum,  W.W.,  et  al.,  Transepithelial  projections  from  basal  cells  are  luminal  sensors  in 
pseudostratified epithelia. Cell, 2008. 135(6): p. 1108-17. 
8.  Shum, W.W.C., et al., Regulation of luminal acidification in the male reproductive tract via cell-
cell crosstalk. Journal of Experimental Biology. 212(11): p. 1753-1761. 
9.  Blomqvist, S.R.,  et  al., Epididymal  expression  of  the  forkhead  transcription  factor Foxi1  is 
required for male fertility. EMBO J, 2006. 25(17): p. 4131-41. 
10.  Distribution  of  immune  cells  in  the  epididymis  of  the  aging Brown Norway  rat  is  segment-
specific and related to the luminal content. 1999. 61(3): p. 705-14. 
11.  Hermann,  B.P.,  et  al.,  The  Mammalian  Spermatogenesis  Single-Cell  Transcriptome,  from 
Spermatogonial Stem Cells to Spermatids. Cell Rep, 2018. 25(6): p. 1650-1667 e8. 
12.  Wang,  M. ,  et  al . ,  Single-Cell RNA Sequencing Analysis Reveals Sequential Cell Fate Transition 
during Human Spermatogenesis. Cell Stem Cell, 2018. 23(4): p. 599-614 e4. 
13.  Kalluri, A.S.,  et  al.,  Single-Cell Analysis  of  the Normal Mouse Aorta  Reveals  Functionally 
Distinct Endothelial Cell Populations. Circulation, 2019. 140(2): p. 147-163. 
14.  Xie, T.,  et  al.,  Single-Cell Deconvolution  of Fibroblast Heterogeneity  in Mouse Pulmonary 
Fibrosis. Cell Rep, 2018. 22(13): p. 3625-3640. 
15.  Guo, J., et al., The adult human testis transcriptional cell atlas. Cell Res, 2018. 28(12): p. 1141-
1157. 
16.  Rebourcet, D., et al., Sertoli cells control peritubular myoid cell fate and support adult Leydig 
cell development in the prepubertal testis. Development, 2014. 141(10): p. 2139-49. 
17.  Mould,  K.J.,  et  al.,  Single  cell  RNA  sequencing  identifies  unique  inflammatory  airspace 
macrophage subsets. JCI Insight, 2019. 4(5). 
18.  Zimmerman, K.A.,  et  al.,  Single-Cell RNA  Sequencing  Identifies Candidate Renal Resident 
Macrophage Gene Expression Signatures across Species. J Am Soc Nephrol, 2019. 30(5): p. 
767-781. 
19.  DeFalco, T., et al., Macrophages Contribute  to  the Spermatogonial Niche  in  the Adult Testis. 
Cell Rep, 2015. 12(7): p. 1107-19. 
20.  Masaki  and  T.,  Heterogeneity  of  antigen  expression  explains  controversy  over  glomerular 
macrophage accumulation in mouse glomerulonephritis. Nephrology Dialysis Transplantation, 
2003. 18(1): p. 178-181. 
21.  Kalucka,  J., et al., Single-Cell Transcriptome Atlas of Murine Endothelial Cells. Cell, 2020. 
180(4): p. 764-779 e20. 
22.  Lukowski, S.W., et al., Single-Cell Transcriptional Profiling of Aortic Endothelium Identifies a 
Hierarchy  from Endovascular Progenitors  to Differentiated Cells. Cell Rep, 2019. 27(9): p. 
2748-2758 e3. 

References

Shi J, Fok KL, Dai P, Qiao F, Zhang M, Liu H, Sang M, Ye M, Liu Y, Zhou Y, Wang C, Sun F, Xie G, Chen H. Spatio-temporal landscape of mouse epididymal cells and specific mitochondria-rich segments defined by large-scale single-cell RNA-seq. Cell Discov. 2021 May 18;7(1):34. doi: 10.1038/s41421-021-00260-7. PMID: 34001862; PMCID: PMC8129088.

Original code

https://github.com/broadinstitute/single_cell_classification

Introduction

H5ad files were uploaded as available data in some paper of scRNA. But I always analysis scRNA data by seurat package in R. As we known, there is a function named read10X_h5 for reading filtered matrix in seurat package. But seurat package does not include any function to read h5ad file, which were from pipeline of scanpy.So some command is needed to solve this problem.

Code exmaple

Installation of reticulate and renv，intall needed package of R and python

# Install renv
install.packages("renv")

#We can then ask {renv} to create an environment for us:
renv::init()

#Install {reticulate} package installed in our new renv environment:
renv::install("reticulate")

#We also want to use a Python environment:
renv::use_python()

#install the R packages used during this tutorial:
pkgs <- c(
    "renv",
    "reticulate",
    "png",
    "ggplot2",
    "BiocManager",
    "Seurat"
)

bioc_pkgs <- c(
    "SingleCellExperiment",
    "scater",
    "multtest"
)

# If you are using an {renv} environment
renv::install(pkgs)
# Install Bioconductor packages
BiocManager::install(bioc_pkgs, update = FALSE)
py_pkgs <- c(
    "scanpy",
    "python-igraph",
    "louvain"
)
reticulate::py_install(py_pkgs)

# This records the changes you have made so they can be restored later if necessary
# When you load R in same work path, this enviroment will be restored automatically .
renv::snapshot()

Load needed packge and data

library(Seurat)
library(data.table)
library(reticulate)
library(renv)
library(ggplot2)

Read scanpy model

sc <- import("scanpy")
setwd("/home/hxzk/project/Qu/Heart/")

Read read h5ad file by scanpy model

annData = sc$read_h5ad("data/figshare_8273102/Heart_droplet.h5ad")

Create seurat object from h5ad data

# Create seurat object
exprs <- t(data.matrix(annData$X))
colnames(exprs) <- annData$obs_names$to_list()
rownames(exprs) <- annData$var_names$to_list()
# Create the Seurat object
seurat <- CreateSeuratObject(exprs)
# Set the expression assay
seurat <- SetAssayData(seurat, "data", exprs)
# Add observation metadata
seurat <- AddMetaData(seurat, annData$obs)
# Zscore normalization if we need
norm_expr <- apply(exprs, 2, function(x) (x - mean(x))/sd(x))
# Add fetaure metadata
seurat[["RNA"]][["n_cells"]] <- annData$var["n_cells"]
# Add embedding
embedding <- annData$obsm["X_umap"]
rownames(embedding) <- annData$obs_names$to_list()
colnames(embedding) <- c("umap_1", "umap_2")
seurat[["umap"]] <- CreateDimReducObject(embedding, key = "umap_")

We can plot some figures as we have transformed h5ad file to seurat object successfully

# Check batch effect
pdf("result/figure/batch_age.pdf", width = 4, height = 3.5)
DimPlot(seurat, reduction = "umap",group.by = "batch")
dev.off()

# Show cluster and annotation
pdf("result/figure/cluster.pdf", width = 12, height = 5)
p1 <- DimPlot(seurat, reduction = "umap",group.by = c("louvain"), label = T)
p2 <- DimPlot(seurat, reduction = "umap",group.by = "free_annotation", label = T)
p1 + p2
dev.off()

References

https://theislab.github.io/scanpy-in-R/#introduction

Donation

Introduction

  Here codes from P. Bischoff was showed. It concluded all step in analysis of scRNA, which could be refered and learned by all of us.

Explanation of codes

Preprocessing, cell type annotation and subsetting

preprocessing

• load cellranger count matrices
• filtering, normalization, dimensionality reduction and clustering
• manual annotation of main cell types (epithelial, immune stromal)
• add clinical metadata
• calculate cell cycle phase
• subset epithelial, immune and stromal cells
• plots for figure 1, supp. fig. 1

cell type annotation

• load epithelial, immune, and stromal subsets
• rerun PCA, dimensionality reduction and clustering
• assign tumor and normal cell clusters in epithelial subset
• calculate expression of marker genes and gene signatures
• manual annotation of epithelial, immune and stromal cell types
• plots for supp. fig. 4

Separate analysis of epithelial, immune and stromal subsets

epithelial analyses

• (use seurat object epi_anno to infer copy number aberrations, code can be found here: https://github.com/bischofp/single_cell_lung_adenocarcinoma)
• epithelial cell type quantification
• selected epithelial marker genes
• subset tumor epithelial cells
• differential gene expression analysis
• infer pathway activities
• rerun PCA and low-dimensional UMAPs
• plot features of tumor cells sorted along PC1
• plots for figures 2, supp. fig. 3

stromal analyses

• stromal cell type quantification
• selected stromal marker genes
• differential gene expression analysis
• infer pathway activities
• plots for figure 3, supp. fig. 5

immune analyses

• immune cell type quantification
• selected immune marker genes
• differential gene expression analysis
• calculate gene signature module scores
• plots for figure 4, supp. fig. 6

Integrative analysis of all subsets

correlation cell types

• load epithelial, immune and stromal subsets
• count cell types
• PCA and correlation analysis of cell counts
• plots for figures 3, 4, 5, supp. fig. 7

TCGA dataset

• load TCGA data
• load epithelial, immune and stromal subsets
• perform differential gene expression analysis
• perform ssGSEA on TCGA data using marker genes from scRNAseq
• correlation analysis of ssGSEA enrichment scores
• stratify patients based on ssGSEA enrichment scores, perform survival analysis
• plots for figure 5, supp. fig. 7

#####################preprocessing

### load libraries
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(ggplot2)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(progeny)
library(readr)
library(stringr)

nFeature_lower <- 500
nFeature_upper <- 10000
nCount_lower <- 1000
nCount_upper <- 100000
pMT_lower <- 0
pMT_upper <- 30
pHB_lower <- 0
pHB_upper <- 5

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467")

# Data loading and QC

### sample list
samples <- read_excel("../data/metadata/patients_metadata.xlsx", range = cell_cols("A:A")) %>% .$sample_id

### import cellranger files from different data sets
for (i in seq_along(samples)){
  assign(paste0("scs_data", i), Read10X(data.dir = paste0("../data/cellranger/", samples[i], "/filtered_feature_bc_matrix")))
}

### create seurat objects from cellranger files
for (i in seq_along(samples)){
  assign(paste0("seu_obj", i), CreateSeuratObject(counts = eval(parse(text = paste0("scs_data", i))), project = samples[i], min.cells = 3))
}

### merge data sets
seu_obj <- merge(seu_obj1, y = c(seu_obj2, seu_obj3, seu_obj4, seu_obj5, seu_obj6, seu_obj7, seu_obj8, seu_obj9, seu_obj10, seu_obj11, seu_obj12, seu_obj13, seu_obj14, seu_obj15, seu_obj16, seu_obj17, seu_obj18, seu_obj19, seu_obj20), add.cell.ids = samples, project = "lung")

### calculate mitochondrial, hemoglobin and ribosomal gene counts
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^MT-", col.name = "pMT")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^HBA|^HBB", col.name = "pHB")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^RPS|^RPL", col.name = "pRP")

qcparams <- c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")
for (i in seq_along(qcparams)){
  print(VlnPlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident", pt.size = 0))
}
for (i in seq_along(qcparams)){
  print(RidgePlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident"))
}

VlnPlot(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"), pt.size = 0, group.by = "orig.ident", ncol = 1, log = T)
ggsave2("SuppFig1B.pdf", path = "../results", width = 30, height = 20, units = "cm")

### clear environment
remove(seu_obj1)
remove(seu_obj2)
remove(seu_obj3)
remove(seu_obj4)
remove(seu_obj5)
remove(seu_obj6)
remove(seu_obj7)
remove(seu_obj8)
remove(seu_obj9)
remove(seu_obj10)
remove(seu_obj11)
remove(seu_obj12)
remove(seu_obj13)
remove(seu_obj14)
remove(seu_obj15)
remove(seu_obj16)
remove(seu_obj17)
remove(seu_obj18)
remove(seu_obj19)
remove(seu_obj20)

remove(scs_data1)
remove(scs_data2)
remove(scs_data3)
remove(scs_data4)
remove(scs_data5)
remove(scs_data6)
remove(scs_data7)
remove(scs_data8)
remove(scs_data9)
remove(scs_data10)
remove(scs_data11)
remove(scs_data12)
remove(scs_data13)
remove(scs_data14)
remove(scs_data15)
remove(scs_data16)
remove(scs_data17)
remove(scs_data18)
remove(scs_data19)
remove(scs_data20)

# Data Filtering 

qc_std_plot_helper <- function(x) x + 
  scale_color_viridis() +
  geom_point(size = 0.01, alpha = 0.3)

qc_std_plot <- function(seu_obj) {
  qc_data <- as_tibble(FetchData(seu_obj, c("nCount_RNA", "nFeature_RNA", "pMT", "pHB", "pRP")))
  plot_grid(

    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pMT))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pHB))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pRP))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),

    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pMT, color = nFeature_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pHB, color = nFeature_RNA))) + 
      geom_hline(yintercept = pHB_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pHB_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pRP, color = nFeature_RNA))) + 
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),

    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pMT, color = nCount_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pHB, color = nCount_RNA))) + 
      geom_hline(yintercept = pHB_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pHB_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pRP, color = nCount_RNA))) + 
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),

    qc_std_plot_helper(ggplot(qc_data, aes(pRP, pMT, color = nCount_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(pRP, pMT, color = nFeature_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2),

    ggplot(gather(qc_data, key, value), aes(key, value)) +
      geom_violin() +
      facet_wrap(~key, scales = "free", ncol = 5),

    ncol = 3, align = "hv"
  )
}

## Before filtering

seu_obj_unfiltered <- seu_obj

#saveRDS(seu_obj_unfiltered, file = "seurat_objects/all_unfiltered.RDS")

qc_std_plot(seu_obj_unfiltered)

ggsave2("SuppFig1A.png", path = "../results", width = 30, height = 30, units = "cm")

## After filtering

#seu_obj_unfiltered <- readRDS("all_unfiltered.RDS")

seu_obj <- subset(seu_obj_unfiltered, subset = nFeature_RNA > nFeature_lower & nFeature_RNA < nFeature_upper & nCount_RNA > nCount_lower & nCount_RNA < nCount_upper & pMT < pMT_upper & pHB < pHB_upper)

qc_std_plot(seu_obj)

seu_obj_unfiltered
seu_obj

# Data normalization

seu_obj <- SCTransform(seu_obj, verbose = T, vars.to.regress = c("nCount_RNA", "pMT"), conserve.memory = T)

#saveRDS(seu_obj, file = "seurat_objects/all_SCTransform.RDS")

# Dimensionality reduction

seu_obj <- RunPCA(seu_obj)
ElbowPlot(seu_obj, ndims = 50)

#for (i in c(15, 20)) {
#  umaptest <- RunUMAP(seu_obj, dims = 1:i, verbose = F)
#  print(DimPlot(umaptest, reduction = "umap", group.by = "orig.ident") + labs(title = paste0(i, " dimensions")))
#  print(FeaturePlot(umaptest, features = c("EPCAM", "PTPRC"), sort.cell = T))
#  print(FeaturePlot(umaptest, features = c("MARCO", "KIT"), sort.cell = T))
#  print(FeaturePlot(umaptest, features = c("FOXJ1", "AGER"), sort.cell = T))
#  print(FeaturePlot(umaptest, features = c("JCHAIN", "VWF"), sort.cell = T))
#  remove(umaptest)
#}

seu_obj <- RunUMAP(seu_obj, dims = 1:15, verbose = T)
seu_obj <- FindNeighbors(seu_obj, dims = 1:15)

for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  seu_obj <- FindClusters(seu_obj, resolution = i)
  print(DimPlot(seu_obj, reduction = "umap") + labs(title = paste0("resolution: ", i)))
}

for (i in c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")) {
  print(FeaturePlot(seu_obj, features = i, coord.fixed = T, sort.cell = T))
}

#DimPlot(seu_obj, group.by = "orig.ident")
#DimPlot(seu_obj, group.by = "SCT_snn_res.0.5", label = T)

# Main cell type annotation

mainmarkers <- c("PECAM1", "VWF", "ACTA2", "JCHAIN", "MS4A1", "PTPRC", "CD68", "KIT", "EPCAM", "CDH1", "KRT7", "KRT19")

for (i in seq_along(mainmarkers)) {
  FeaturePlot(seu_obj, features = mainmarkers[i], coord.fixed = T, order = T, cols = viridis(10))
  #ggsave2(paste0("FeaturePlot_mainmarkers_", mainmarkers[i], ".png"), path = "output/annotation", width = 10, height = 10, units = "cm")
}

DotPlot(seu_obj, features = mainmarkers, group.by = "SCT_snn_res.0.2") + 
  coord_flip() + 
  scale_color_viridis()
#ggsave2("DotPlot_mainmarkers.png", path = "output/annotation", width = 30, height = 8, units = "cm")

DimPlot(seu_obj, group.by = "SCT_snn_res.0.2", label = T, label.size = 5)
#ggsave2("DimPlot_all_clusters.png", path = "output/annotation", width = 20, height = 20, units = "cm")

Idents(seu_obj) <- seu_obj$SCT_snn_res.0.2
annotation_curated_main <- read_excel("../data/curated_annotation/curated_annotation_main.xlsx")
new_ids_main <- annotation_curated_main$main_cell_type
names(new_ids_main) <- levels(seu_obj)
seu_obj <- RenameIdents(seu_obj, new_ids_main)
seu_obj@meta.data$main_cell_type <- Idents(seu_obj)

# Add metadata

metatable <- read_excel("../data/metadata/patients_metadata.xlsx")

metadata <- FetchData(seu_obj, "orig.ident")
metadata$cell_id <- rownames(metadata)
metadata$sample_id <- metadata$orig.ident
metadata <- left_join(x = metadata, y = metatable, by = "sample_id")
rownames(metadata) <- metadata$cell_id

seu_obj <- AddMetaData(seu_obj, metadata = metadata)

# Cell cycle scoring

### add cell cycle, cc.genes loaded with Seurat

s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

score_cc <- function(seu_obj) {
  seu_obj <- CellCycleScoring(seu_obj, s.genes, g2m.genes)
  seu_obj@meta.data$CC.Diff <- seu_obj@meta.data$S.Score - seu_obj@meta.data$G2M.Score
  return(seu_obj)
}

seu_obj <- score_cc(seu_obj)

FeatureScatter(seu_obj, "G2M.Score", "S.Score", group.by = "Phase", pt.size = .1) +
  coord_fixed(ratio = 1)

# Subset, rescale and save RDS files

###subset and rescale

#saveRDS(seu_obj, file = "seurat_objects/all.RDS")

Idents(seu_obj) <- seu_obj@meta.data$main_cell_type

epi <- subset(seu_obj, idents = "Epithelial")
imm <- subset(seu_obj, idents = "Immune")
str <- subset(seu_obj, idents = "Stromal")

epi <- ScaleData(epi)
imm <- ScaleData(imm)
str <- ScaleData(str)

epi
imm
str

#saveRDS(epi, file = "seurat_objects/epi.RDS")
#saveRDS(imm, file = "seurat_objects/imm.RDS")
#saveRDS(str, file = "seurat_objects/str.RDS")

# Plots for figure 1

###r plots for figure 1

DimPlot(seu_obj, group.by = "tissue_type", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1B.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(seu_obj, group.by = "patient_id", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1C.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(seu_obj, group.by = "main_cell_type", cols = use_colors, pt.size = 0.1)
ggsave2("Fig1D_umap.png", path = "../results", width = 15, height = 15, units = "cm")

cell_types <- FetchData(seu_obj, vars = c("sample_id", "main_cell_type", "tissue_type")) %>% 
  mutate(main_cell_type = factor(main_cell_type, levels = c("Stromal", "Immune", "Epithelial"))) %>% 
  mutate(sample_id = factor(sample_id, levels = rev(c("p018t", "p019t", "p023t", "p024t", "p027t", "p028t", "p030t", "p031t", "p032t", "p033t", "p034t", "p018n", "p019n", "p027n", "p028n", "p029n", "p030n", "p031n", "p032n", "p033n", "p034n"))))

ggplot(data = cell_types) + 
  geom_bar(mapping = aes(x = sample_id, fill = main_cell_type, ), position = "fill", width = 0.75) +
  scale_fill_manual(values = use_colors) +
  coord_flip()
ggsave2("Fig1D_barplot.pdf", path = "../results", width = 15, height = 30, units = "cm")

remove(seu_obj)

##########################################cell type annotation

### load libraries
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(ggplot2)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(readr)
library(stringr)
library(progeny)
library(scales)

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467")

# Load data

#epi <- readRDS("seurat_objects/epi.RDS")
#imm <- readRDS("seurat_objects/imm.RDS")
#str <- readRDS("seurat_objects/str.RDS")

# Rerun PCA, reclustering

### epithelial subclustering

epi <- RunPCA(epi)
ElbowPlot(epi,  ndims = 50)

#for (i in c(10, 15, 20, 25)){
#  umaptest <- RunUMAP(epi, dims = 1:i, verbose = F)
#  print(DimPlot(umaptest, reduction = "umap", group.by = "patient_id", split.by = "tissue_type") + labs(title = paste0(i, "dimensions")))
#  remove(umaptest)
#}

epi <- RunUMAP(epi, dims = 1:20)
epi <- FindNeighbors(epi, dims = 1:20)
for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  epi <- FindClusters(epi, resolution = i)
  print(DimPlot(epi, reduction = "umap", label = T) + labs(title = paste0("resolution: ", i)))
}

Idents(epi) <- epi@meta.data$SCT_snn_res.1

### immune subclustering

imm <- RunPCA(imm)
ElbowPlot(imm,  ndims = 50)

#for (i in c(10, 15, 20, 25)){
#  umaptest <- RunUMAP(imm, dims = 1:i, verbose = F)
#  print(DimPlot(umaptest, reduction = "umap", group.by = "patient_id", split.by = "tissue_type") + labs(title = paste0(i, " dimensions")))
#  remove(umaptest)
#}

imm <- RunUMAP(imm, dims = 1:20)
imm <- FindNeighbors(imm, dims = 1:20)
for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  imm <- FindClusters(imm, resolution = i)
  print(DimPlot(imm, reduction = "umap") + labs(title = paste0("resolution: ", i)))
}

Idents(imm) <- imm@meta.data$SCT_snn_res.0.5

### stromal sublustering

str <- RunPCA(str)
ElbowPlot(str, ndims = 50)

#for (i in c(5, 10, 15, 20, 25, 30)){
#  umaptest <- RunUMAP(str, dims = 1:i, verbose = F)
#  print(DimPlot(umaptest, reduction = "umap", group.by = "patient_id", split.by = "tissue_type") + labs(title = paste0(i, " dimensions")))
#  print(DimPlot(umaptest, reduction = "umap", group.by = "tissue_type") + labs(title = paste0(i, " dimensions")))
#  remove(umaptest)
#}

str <- RunUMAP(str, dims = 1:20)
str <- FindNeighbors(str, dims = 1:20)
for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  str <- FindClusters(str, resolution = i)
  print(DimPlot(str, reduction = "umap") + labs(title = paste0("resolution: ", i)))
}

Idents(str) <- str@meta.data$SCT_snn_res.1

# Define normal and tumor cell clusters

DimPlot(epi, group.by = "SCT_snn_res.1", label = T, repel = T, split.by = "tissue_type")
ggsave2("SuppFig3A.png", path = "../results", width = 30, height = 15, units = "cm")

### compare proportion of cells in a cluster to all epithelial cells for tumor and normal separately, clusters overrepresented in normal samples are supposed to be cell of normal lung parenchyma, all other clusters are supposed to be tumor cells

epi_clusters <- FetchData(epi, vars = c("SCT_snn_res.1", "tissue_type"))

count_tumor <- epi_clusters %>% filter(tissue_type == "Tumor") %>% count() %>% as.numeric()
count_normal <- epi_clusters %>% filter(tissue_type == "Normal") %>% count() %>% as.numeric()

epi_counts <- epi_clusters %>% group_by(tissue_type) %>% count(SCT_snn_res.1)
proportion_tumor <- epi_counts %>% filter(tissue_type == "Tumor") %>% mutate(proportion = n/count_tumor)
proportion_normal <- epi_counts %>% filter(tissue_type == "Normal") %>% mutate(proportion = n/count_normal)

proportion_epi <- full_join(proportion_normal, proportion_tumor, by = "SCT_snn_res.1") %>% 
  mutate(proportion.x = ifelse(is.na(proportion.x), 0,  proportion.x)) %>%  
  mutate(proportion.y = ifelse(is.na(proportion.y), 0,  proportion.y)) %>%
  mutate(tissue_type.x = "Normal") %>%
  mutate(tissue_type.y = "Tumor") %>%
  mutate(cluster_type = ifelse(proportion.x > proportion.y, "Normal", "Tumor"))

cluster_type_data <- left_join(x = epi_clusters, y = proportion_epi, by = "SCT_snn_res.1")
rownames(cluster_type_data) <- rownames(epi_clusters)

epi <- AddMetaData(epi, select(cluster_type_data, cluster_type))

### Bar plot for figure 2

n1 <- select(proportion_epi, c(tissue_type.x, SCT_snn_res.1, proportion.x)) %>%
  mutate(tissue_type = tissue_type.x) %>% 
  mutate(proportion = proportion.x) %>%
  mutate(tissue_type.x = NULL) %>%
  mutate(proportion.x = NULL)
t1 <- select(proportion_epi, c(tissue_type.y, SCT_snn_res.1, proportion.y)) %>%
  mutate(tissue_type = tissue_type.y) %>% 
  mutate(proportion = proportion.y) %>%
  mutate(tissue_type.y = NULL) %>%
  mutate(proportion.y = NULL)

proportion_epi2 <- rbind(n1, t1)

ggplot(proportion_epi2, aes(fill = tissue_type, y = proportion, x = SCT_snn_res.1)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  scale_fill_manual(values = use_colors)

ggsave2("SuppFig3B.pdf", path = "../results", width = 40, height = 20, units = "cm")

# Cell type scoring

myeloid_markers <- c("S100A12", "FCN1", "S100A8", "S100A9", "CD14", "CTSS","VCAN", "LYZ", 
                     "MARCO", "FCGR1A", "C1QA", "APOC1", "LGMN", "CTSB", "FCGR3A", 
                     "MAFB", "MAF", "CX3CR1", "ITGAM", "CSF1R",
                     "FABP4", "MCEMP1", 
                     "IL1B", "CXCL8", 
                     "APOE", "CD163", "C1QB", "C1QC", 
                     "FCER1A", "CD1C", "CLEC9A", 
                     "LILRA4", "CLEC4C", "JCHAIN", "IL3RA", "NRP1", 
                     "CLEC10A", "PTCRA", "CCR7", "LAMP3", 
                     "ITGAX", "CD68", "MKI67", "CDK1", "EPCAM")

tcell_nk_markers <- c("CD3E", "CD4", "FOXP3", "IL7R", "IL2RA", "CD40LG", "CD8A", "CCL5", "NCR1", "NKG7", "GNLY", "NCAM1", "KLRD1", "KLRB1", "CD69", "KLRG1", "MKI67", "CDK1", "EPCAM")

bcell_plasma_mast_markers <- c("MS4A1", "CD19", "CD79A", "JCHAIN", "IGHA1", "IGHA2", "IGHG1", "IGHG2", "IGHG3", "IGHG4", "IGKC", "IGLC2", "IGLC3", "CPA3", "KIT", "MS4A2", "GATA2",  "MKI67", "CDK1", "EPCAM")

DotPlot(imm, features = myeloid_markers, group.by = "SCT_snn_res.0.5") + coord_flip()
#ggsave2("DotPlot_myeloid_markers.png", path = "output", width = 20, height = 20, units = "cm")

DotPlot(imm, features = tcell_nk_markers, group.by = "SCT_snn_res.0.5") + coord_flip()
#ggsave2("DotPlot_T_NK_markers.png", path = "output", width = 20, height = 20, units = "cm")

DotPlot(imm, features = bcell_plasma_mast_markers, group.by = "SCT_snn_res.0.5") + coord_flip()
#ggsave2("DotPlot_B_Plasma_markers.png", path = "output", width = 20, height = 20, units = "cm")

DimPlot(imm, group.by = "SCT_snn_res.0.5", label = T, split.by = "tissue_type")
#ggsave2("DimPlot_immune_clusters.png", path = "output", width = 30, height = 15, units = "cm")

DimPlot(imm, group.by = "patient_id", split.by = "tissue_type", cols = use_colors)
#ggsave2("DimPlot_immune_patients.png", path = "output", width = 30, height = 15, units = "cm")

## Habermann et al.
#https://www.biorxiv.org/content/10.1101/753806v1

habermann_epi <- c("ABCA3", "SFTPB", "SFTPC", "AGER", "PDPN",  "KRT5", "TRP63", "NGFR", "SCGB1A1", "MUC5B", "KRT17", "FOXJ1", "TMEM190", "CAPS", "CHGA", "CALCA", "ASCL1", "PTPRC", "EPCAM")

habermann_imm <- c("CD3E", "CD4", "FOXP3", "IL7R", "IL2RA", "CD40LG", "CD8A", "CCL5", "NCR1", "KLRB1", "NKG7", "LYZ", "CD68", "ITGAX", "MARCO", "FCGR1A", "FCGR3A", "C1QA", "APOC1", "S100A12", "FCN1", "S100A9", "CD14", "FCER1A", "CD1C", "CD16", "CLEC9A", "LILRA4", "CLEC4C", "JCHAIN", "IGHG1", "IGLL5", "MS4A1", "CD19", "CD79A", "CPA3", "KIT", "MKI67", "CDK1", "EPCAM")

habermann_oth <- c("VWF", "PECAM1", "CCL21", "PROX1", "ACTA2", "MYH11", "PDGFRB", "WT1", "UPK3B", "LUM", "PDGFRA", "MYLK", "HAS1", "PLIN2", "FAP", "PTPRC", "EPCAM")

### Epithelial genes according to Habermann et al.

### manual annotation of normal cell types & detection of immune cell contaminated clusters

DotPlot(SubsetData(epi, subset.name = "cluster_type", accept.value = "Normal"), features = habermann_epi, group.by = "SCT_snn_res.1") + theme(axis.text.x = element_text(angle = 90, hjust = 1))

### manual detection of immune cell contaminated clusters

DotPlot(SubsetData(epi, subset.name = "cluster_type", accept.value = "Tumor"), features = habermann_epi, group.by = "SCT_snn_res.1") + theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Immune genes according to Habermann et al.

DotPlot(imm, features = habermann_imm, group.by = "SCT_snn_res.0.5") + coord_flip()

#for (i in seq_along(habermann_imm)) {
#  plotlist <- list()
#  plotlist[1] <- FeaturePlot(imm, features = habermann_imm[i], sort.cell = T, combine = F)
#  plotlist[2] <- VlnPlot(imm, features = habermann_imm[i], pt.size = 0, combine = F)
#  print(CombinePlots(plots = plotlist))
#}

### Stromal genes according to Habermann et al.

DotPlot(str, features = habermann_oth, group.by = "SCT_snn_res.1") + 
  coord_flip()

#for (i in seq_along(habermann_oth)) {
#  plotlist <- list()
#  plotlist[1] <- FeaturePlot(str, features = habermann_oth[i], sort.cell = T, combine = F)
#  plotlist[2] <- VlnPlot(str, features = habermann_oth[i], pt.size = 0, combine = F)
#  print(CombinePlots(plots = plotlist, ncol = 3))
#}

## Travaglini et al.
#https://www.biorxiv.org/content/10.1101/742320v1

#load marker gene lists

sheets <- paste0("Cluster ", c(1:58))
sheets <- sheets[-43]

signaturelist <- list()

for (i in seq_along(sheets)) {
  a <- read_excel("../data/data/media-3.xlsx", sheet = sheets[[i]])
  a <- filter(a, a$...2 > 0.7 & a$...4 < 0.3)
  signaturelist <- c(signaturelist,a[1])
  remove(a)
}

#generate list with names of module scores in seurat object

names_of_modulescores <- c()
for (i in seq_along(signaturelist)){
  names_of_modulescores <- c(names_of_modulescores, paste0("T_", names(signaturelist[i]), i))
}

names_of_modulescores <- gsub(names_of_modulescores, pattern = " ", replacement = ".", fixed = TRUE)
names_of_modulescores <- gsub(names_of_modulescores, pattern = "+", replacement = ".", fixed = TRUE)
names_of_modulescores <- gsub(names_of_modulescores, pattern = "/", replacement = ".", fixed = TRUE)

#names_of_modulescores_unfiltered <- c()
#for (i in seq_along(signaturelist)){
#  names_of_modulescores_unfiltered <- c(names_of_modulescores_unfiltered, paste0(names(signaturelist[i]), "_unfiltered", i))
#}
#names_of_modulescores_unfiltered <- gsub(names_of_modulescores_unfiltered, pattern = " ", replacement = ".", fixed = TRUE)
#names_of_modulescores_unfiltered <- gsub(names_of_modulescores_unfiltered, pattern = "+", replacement = ".", fixed = TRUE)
#names_of_modulescores_unfiltered <- gsub(names_of_modulescores_unfiltered, pattern = "/", replacement = ".", fixed = TRUE)
#signature_list_updated <- list()
#for (i in seq_along(sheets)) {
#  signature_list_updated[[i]] <- checkGeneSymbols(signature_list[[i]])
#}

#calculate module scores for different subsets

epi <- AddModuleScore(object = epi, features = signaturelist, name = paste0("T_", names(signaturelist)))
imm <- AddModuleScore(object = imm, features = signaturelist, name = paste0("T_", names(signaturelist)))
str <- AddModuleScore(object = str, features = signaturelist, nbin = 12 , name = paste0("T_", names(signaturelist)))

## Vieira Braga et al.
#https://www.nature.com/articles/s41591-019-0468-5

#load marker gene lists

teichmann_signatures_epi <- read.csv("../data/data/Fig1_DE_Lung_atlas_epithelial.csv")
teichmann_signatures_epi$gene <- as.character(teichmann_signatures_epi$gene)
teichmann_epi <- levels(teichmann_signatures_epi$cluster)

teichmann_signatures_imm <- read.csv("../data/data/Fig2_DE_Lung_atlas_immune.csv")
teichmann_signatures_imm$gene <- as.character(teichmann_signatures_imm$gene)
teichmann_imm <- levels(teichmann_signatures_imm$cluster)

signaturelist2 <- list()

for (i in seq_along(teichmann_epi)) {
  signaturelist2 <- c(signaturelist2, teichmann_signatures_epi %>% filter(cluster == teichmann_epi[i], pct.2 < 0.3, avg_logFC > 0.7) %>% select(gene))
}

for (i in seq_along(teichmann_imm)) {
  signaturelist2 <- c(signaturelist2, teichmann_signatures_imm %>% filter(cluster == teichmann_imm[i], pct.2 < 0.3, avg_logFC > 0.7) %>% select(gene))
}

names(signaturelist2) <- gsub(c(teichmann_epi, teichmann_imm), pattern = "_", replacement = " ")

#generate list with names of module scores in seurat object

names_of_modulescores2 <- c()
for (i in seq_along(signaturelist2)){
  names_of_modulescores2 <- c(names_of_modulescores2, paste0("VB_", names(signaturelist2[i]), i))
}
names_of_modulescores2 <- gsub(names_of_modulescores2, pattern = " ", replacement = ".", fixed = TRUE)

#calculate module scores for different subsets

epi <- AddModuleScore(epi, features = signaturelist2, name = paste0("VB_", names(signaturelist2)))
imm <- AddModuleScore(imm, features = signaturelist2, name = paste0("VB_", names(signaturelist2)))
str <- AddModuleScore(str, features = signaturelist2, nbin = 12, name = paste0("VB_", names(signaturelist2)))

# Curated cell type annotation and subsetting

###epithelial

annotation_curated_epi <- read_excel("../data/curated_annotation/curated_annotation_epi.xlsx")
epi_anno <- epi
new_ids_epi <- annotation_curated_epi$cell_type_epi
names(new_ids_epi) <- levels(epi_anno)
epi_anno <- RenameIdents(epi_anno, new_ids_epi)
epi_anno@meta.data$cell_type_epi <- Idents(epi_anno)

epi_anno <- subset(epi_anno, subset = cell_type_epi != "Immune_contamination")
epi_anno <- ScaleData(epi_anno)

###immune

annotation_curated_imm <- read_excel("../data/curated_annotation/curated_annotation_imm.xlsx")
imm_anno <- imm
new_ids_imm <- annotation_curated_imm$cell_type_imm
names(new_ids_imm) <- levels(imm_anno)
imm_anno <- RenameIdents(imm_anno, new_ids_imm)
imm_anno@meta.data$cell_type_imm <- Idents(imm_anno)

imm_anno <- subset(imm_anno, subset = cell_type_imm != "Epithelial_contamination")
imm_anno <- ScaleData(imm_anno)

###stromal

annotation_curated_str <- read_excel("../data/curated_annotation/curated_annotation_str.xlsx")
str_anno <- str
new_ids_str <- annotation_curated_str$cell_type_str
names(new_ids_str) <- levels(str_anno)
str_anno <- RenameIdents(str_anno, new_ids_str)
str_anno@meta.data$cell_type_str <- Idents(str_anno)

str_anno <- subset(str_anno, subset = cell_type_str != "Immune/Epithelial contamination")
str_anno <- ScaleData(str_anno)

###Travaglini et al.

names_of_modulescores_original <- names(signaturelist)

#Epithelial

epi_type <- FetchData(epi_anno, vars = c(names_of_modulescores))

for(i in seq_along(names_of_modulescores_original)) {
  colnames(epi_type)[i] <- names_of_modulescores_original[i]
}

epi_type %>%
  merge(FetchData(epi_anno, vars = c("cluster_type", "cell_type_epi")), by = 0) %>%
  filter(cluster_type == "Normal") %>%
  group_by(cell_type_epi) %>%
  pivot_longer(cols = names_of_modulescores_original, names_to = "T_cell_type") %>%
  group_by(cell_type_epi, T_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = T_cell_type, y = cell_type_epi, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5B_epithelial.pdf", path = "../results", width = 30, height = 15, units = "cm")

#Immune

imm_type <- FetchData(imm_anno, vars = c(names_of_modulescores))

for(i in seq_along(names_of_modulescores_original)) {
  colnames(imm_type)[i] <- names_of_modulescores_original[i]
}

imm_type %>%
  merge(FetchData(imm_anno, vars = "cell_type_imm"), by = 0) %>%
  group_by(cell_type_imm) %>%
  pivot_longer(cols = names_of_modulescores_original, names_to = "T_cell_type") %>%
  group_by(cell_type_imm, T_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = T_cell_type, y = cell_type_imm, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5B_immune.pdf", path = "../results", width = 30, height = 30, units = "cm")

#Stromal

str_type <- FetchData(str_anno, vars = c(names_of_modulescores))

for(i in seq_along(names_of_modulescores_original)) {
  colnames(str_type)[i] <- names_of_modulescores_original[i]
}

str_type %>%
  merge(FetchData(str_anno, vars = "cell_type_str"), by = 0) %>%
  group_by(cell_type_str) %>%
  pivot_longer(cols = names_of_modulescores_original, names_to = "T_cell_type") %>%
  group_by(cell_type_str, T_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = T_cell_type, y = cell_type_str, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5B_stromal.pdf", path = "../results", width = 30, height = 20, units = "cm")

###Vieira Braga et al.

names_of_modulescores_original2 <- names(signaturelist2)

#Epithelial

epi_type <- FetchData(epi_anno, vars = c(names_of_modulescores2))

for(i in seq_along(names_of_modulescores_original2)) {
  colnames(epi_type)[i] <- names_of_modulescores_original2[i]
}

epi_type %>%
  merge(FetchData(epi_anno, vars = c("cluster_type", "cell_type_epi")), by = 0) %>%
  filter(cluster_type == "Normal") %>%
  group_by(cell_type_epi) %>%
  pivot_longer(cols = names_of_modulescores_original2, names_to = "VB_cell_type") %>%
  group_by(cell_type_epi, VB_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = VB_cell_type, y = cell_type_epi, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5A_epithelial.pdf", path = "../results", width = 20, height = 10, units = "cm")

#Immune

imm_type <- FetchData(imm_anno, vars = c(names_of_modulescores2))

for(i in seq_along(names_of_modulescores_original2)) {
  colnames(imm_type)[i] <- names_of_modulescores_original2[i]
}

imm_type %>%
  merge(FetchData(imm_anno, vars = "cell_type_imm"), by = 0) %>%
  group_by(cell_type_imm) %>%
  pivot_longer(cols = names_of_modulescores_original2, names_to = "VB_cell_type") %>%
  group_by(cell_type_imm, VB_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = VB_cell_type, y = cell_type_imm, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5A_immune.pdf", path = "../results", width = 20, height = 30, units = "cm")

#Stromal

str_type <- FetchData(str_anno, vars = c(names_of_modulescores2))

for(i in seq_along(names_of_modulescores_original2)) {
  colnames(str_type)[i] <- names_of_modulescores_original2[i]
}

str_type %>%
  merge(FetchData(str_anno, vars = "cell_type_str"), by = 0) %>%
  group_by(cell_type_str) %>%
  pivot_longer(cols = names_of_modulescores_original2, names_to = "VB_cell_type") %>%
  group_by(cell_type_str, VB_cell_type) %>%
  summarise(mean = mean(value)) %>%
  mutate(mean = rescale(mean)) %>%
  ggplot() +
  geom_tile(aes(x = VB_cell_type, y = cell_type_str, fill = mean))+
  scale_fill_gradientn(colors = c("blue", "white", "red"),
                       breaks = c(0, 1),
                       labels = c("0", "1")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), axis.text.y = element_text(hjust = 1))
ggsave2("SuppFig5A_stromal.pdf", path = "../results", width = 20, height = 20, units = "cm")

#hallmark signatures

broad_pws <- read_lines("../data/data/h.all.v6.2.symbols.gmt") %>%
  lapply(str_split, "\\t") %>% 
  unlist(recursive = F) %>% 
  lapply(function(x) setNames(list(x[-c(1:2)]), x[1])) %>% 
  unlist(recursive = F)

epi_anno <- AddModuleScore(object = epi_anno, features = broad_pws, name = names(broad_pws))
imm_anno <- AddModuleScore(object = imm_anno, features = broad_pws, name = names(broad_pws))
str_anno <- AddModuleScore(object = str_anno, features = broad_pws, name = names(broad_pws), nbin = 12)

#progeny signatures

epi_anno <- progeny(epi_anno, scale = F, organism="Human", top=500, perm=1, return_assay=T)
epi_anno <- ScaleData(epi_anno, assay = "progeny")

### obviously too many cells for progeny function, progeny scores added after immune cells split into myeloid and lymphoid subset
#imm_anno <- progeny(imm_anno, scale = F, organism="Human", top=500, perm=1, return_assay=T)
#imm_anno <- ScaleData(imm_anno, assay = "progeny")

str_anno <- progeny(str_anno, scale = F, organism="Human", top=500, perm=1, return_assay=T)
str_anno <- ScaleData(str_anno, assay = "progeny")

#saveRDS(epi_anno, file = "seurat_objects/epi_anno.RDS")
#saveRDS(imm_anno, file = "seurat_objects/imm_anno.RDS")
#saveRDS(str_anno, file = "seurat_objects/str_anno.RDS")
remove(epi)
remove(imm)
remove(str)

#########################################epithelial analyses

### load libraries
library(ggplot2)
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(stringr)
library(cowplot)
library(scales)
library(tibble)
library(gplots)
library(RColorBrewer)

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467",
  AT1 = "#2B8CBE",
  AT2 = "#045A8D",
  Club = "#006D2C",
  DifferentiatingCiliated = "#31A354",
  Ciliated = "#74C476",
  Neuroendocrine = "#8856A7",
  lepidic = "#0B775E",
  acinar = "#74A089",
  mucinuous (papillary) = "#E2D200",
  (micro)papillary = "#CEAB07",
  solid = "#B40F20",
  sarcomatoid = "#5B1A18",
  CNN = "chartreuse4",
  CNA = "orange")

#load data

#epi_anno <- readRDS("seurat_objects/epi_anno.RDS")

epi_anno@meta.data$cell_type_epi <- factor(epi_anno@meta.data$cell_type_epi, levels = c("AT2",
                                                                                        "AT1",
                                                                                        "Club",
                                                                                        "Ciliated",
                                                                                        "Neuroendocrine",
                                                                                        "Tumor"))

#add inferCNV clone scores (for inferCNV, use seurat object epi_anno and the following code: https://github.com/bischofp/single_cell_lung_adenocarcinoma, computation time ~12h)

scna_scores <- read.delim("../data/inferCNV_output/infercnv_clone_scores_nsclc.tsv") %>% filter(tissue_type == "Tumor") %>% filter(!is.na(cna_clone))
rownames(scna_scores) <- scna_scores$cell_id
scna_scores <- scna_scores %>% select(cna_clone) %>% mutate(cna_clone = as.character(cna_clone))
epi_anno <- AddMetaData(epi_anno, scna_scores)

scna_data <- FetchData(epi_anno, c("tissue_type", "cna_clone"))
scna_data <- scna_data %>% mutate(cna_clone = ifelse(is.na(cna_clone), "CNN", cna_clone))
epi_anno <- AddMetaData(epi_anno, scna_data)

#some plots

DotPlot(SubsetData(epi_anno, subset.name = "cluster_type", accept.value = "Normal"), features = c("ABCA3", "SFTPC", "AGER", "PDPN",  "KRT5", "TRP63", "NGFR", "SCGB1A1", "MUC5B", "FOXJ1", "TMEM190", "CHGA", "CALCA"), group.by = "cell_type_epi") + 
  coord_flip() + 
  scale_color_viridis()
#ggsave2("DotPlot_markergenes_epi_cell_type.emf", path = "../results", width = 11, height = 8, units = "cm")
ggsave2("Fig2B.png", path = "../results", width = 11, height = 8, units = "cm")

DimPlot(epi_anno, group.by = "cell_type_epi", cols = use_colors, pt.size = 0.5)
ggsave2("Fig2A_celltype.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(epi_anno, group.by = "patient_id", cols = use_colors, pt.size = 0.5)
ggsave2("Fig2A_patients.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(epi_anno, group.by = "tissue_type", cols = use_colors, pt.size = 0.5)
ggsave2("Fig2A_tissuetype.png", path = "../results", width = 15, height = 15, units = "cm")

epi_cell_counts <- FetchData(epi_anno, vars = c("tissue_type", "cell_type_epi", "cna_clone", "cluster_type")) %>%
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal")))

ggplot(data = epi_cell_counts, aes(x = tissue_type, fill = cell_type_epi)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  scale_y_reverse() +
  coord_flip()
ggsave2("Fig2A_barplot.pdf", path = "../results", width = 20, height = 5, units = "cm")

ggplot(data = epi_cell_counts, aes(x = tissue_type, fill = cna_clone)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  scale_y_reverse() +
  coord_flip()
ggsave2("SuppFig4_barplot_cna_clone.pdf", path = "../results", width = 20, height = 5, units = "cm")

ggplot(data = epi_cell_counts, aes(x = tissue_type, fill = cluster_type)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("cyan3", "darkorange2")) +
  coord_flip()
ggsave2("SuppFig4_barplot_cluster_type.pdf", path = "../results", width = 20, height = 5, units = "cm")

DimPlot(epi_anno, group.by = "cna_clone", cols = use_colors, pt.size = 0.5)
ggsave2("SuppFig4_umap_cna_clone.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(epi_anno, group.by = "cluster_type", cols = c("cyan3", "darkorange2"), pt.size = 0.5)
ggsave2("SuppFig4_umap_cluster_type.png", path = "../results", width = 15, height = 15, units = "cm")

###subset tumor epithelial cells

epi_tumor <- SubsetData(SubsetData(epi_anno, subset.name = "cluster_type", accept.value = "Tumor"), subset.name = "tissue_type", accept.value = "Tumor")

epi_tumor <- ScaleData(epi_tumor)

###tumor cell marker genes

Idents(epi_tumor) <- epi_tumor@meta.data$patient_id

markers <- FindAllMarkers(epi_tumor, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)

top_TC_markers <- markers %>% group_by(cluster) %>% top_n(10, wt = avg_logFC)

DoHeatmap(epi_tumor, features = top_TC_markers$gene, group.by = "patient_id", draw.lines = F, group.colors = use_colors) +
  scale_fill_viridis()
#ggsave2("HeatMap_Tumor.pdf", path = "output/fig2", width = 30, height = 30, units = "cm")
ggsave2("Fig2C.png", path = "../results", width = 30, height = 30, units = "cm")

###mitotic activity

mitotic_activity <- FetchData(epi_tumor, c("tissue_type", "cell_type_epi", "Phase", "sample_id")) %>%
  mutate(sample_id = factor(sample_id, levels = c("p034t", "p033t", "p032t", "p031t", "p030t", "p027t", "p024t", "p023t", "p019t", "p018t")))

ggplot(mitotic_activity, aes(x = sample_id, fill = Phase)) +
  geom_bar(position = "fill", width = 0.75) +
  scale_fill_manual(values = use_colors) +
  coord_flip()
ggsave2("SuppFig3D.pdf", path = "../results", width = 12, height = 10, units = "cm")

###progeny pathway scores

#clustered heatmap progeny scores

progeny_scores <- as.data.frame(t(GetAssayData(epi_tumor, assay = "progeny", slot = "scale.data")))
progeny_scores$cell_id <- rownames(progeny_scores)
progeny_scores <- gather(progeny_scores, Pathway, Activity, -cell_id)

cells_clusters <- FetchData(epi_tumor, c("sample_id", "cluster_type")) %>% filter(str_detect(sample_id, "t"))
cells_clusters$cell_id <- rownames(cells_clusters)

progeny_scores <- inner_join(progeny_scores, cells_clusters)

summarized_progeny_scores <- progeny_scores %>% 
  group_by(Pathway, sample_id) %>% 
  summarise(avg = mean(Activity), std = sd(Activity)) %>%
  pivot_wider(id_cols = Pathway, names_from = sample_id, values_from = avg) %>%
  column_to_rownames("Pathway") %>%
  as.matrix()

pdf("../results/Fig2E.pdf", width = 6, height = 8)
heatmap.2(summarized_progeny_scores, trace = "none", density.info = "none", col = bluered(100))
dev.off()

###correlation mutational status ~ progeny scores

summarized_progeny_scores_mutations <- summarized_progeny_scores %>%
  t() %>%
  as.data.frame()

summarized_progeny_scores_mutations$Sample <- rownames(summarized_progeny_scores_mutations)
summarized_progeny_scores_mutations <- summarized_progeny_scores_mutations %>%
  mutate(KRAS_status = ifelse(Sample %in% c("p018t", "p023t", "p030t", "p031t", "p032t", "p033t"), "mutated", "wildtype")) %>%
  mutate(TP53_status = ifelse(Sample %in% c("p023t", "p027t"), "mutated", "wildtype")) %>%
  mutate(PIK3CA_status = ifelse(Sample %in% c("p031t"), "mutated", "wildtype")) %>%
  mutate(KRAS_status = factor(TP53_status, levels = c("wildtype", "mutated"))) %>%
  mutate(TP53_status = factor(TP53_status, levels = c("wildtype", "mutated")))

ggplot(summarized_progeny_scores_mutations, aes(x = KRAS_status, y = MAPK)) +
  geom_boxplot() +
  ggtitle(paste0(t.test(formula = MAPK~KRAS_status, data = summarized_progeny_scores_mutations, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig3C_KRAS~MAPK.pdf", path = "../results", width = 8, height = 8, units = "cm")

ggplot(summarized_progeny_scores_mutations, aes(x = KRAS_status, y = EGFR)) +
  geom_boxplot() +
  ggtitle(paste0(t.test(formula = EGFR~KRAS_status, data = summarized_progeny_scores_mutations, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig3C_KRAS~EGFR.pdf", path = "../results", width = 8, height = 8, units = "cm")

ggplot(summarized_progeny_scores_mutations, aes(x = TP53_status, y = p53)) +
  geom_boxplot() +
  ggtitle(paste0(t.test(formula = p53~TP53_status, data = summarized_progeny_scores_mutations, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig3C_p53~TP53.pdf", path = "../results", width = 8, height = 8, units = "cm")

### rerun PCA

epi_pca <- epi_tumor

epi_pca <- RunPCA(epi_pca)

DimPlot(epi_pca, reduction = "pca", group.by = "patient_id", dims = c(1,2))
DimPlot(epi_pca, reduction = "pca", group.by = "patient_id", dims = c(3,4))

DimPlot(epi_pca, reduction = "pca", group.by = "histo_subtype")

DimHeatmap(epi_pca, dims = 1, cells = 1000, balanced = T, fast = F, nfeatures = 60) +
  scale_fill_viridis()
#ggsave2("DimHeatmap_epitumor_PC1.pdf", path = "output/fig2", width = 10, height = 20, units = "cm")
ggsave2("SuppFig3F.png", path = "../results", width = 10, height = 20, units = "cm")

### low-dimensional UMAPs

DimPlot(epi_pca, group.by = "patient_id")
DimPlot(epi_pca, group.by = "histo_subtype")

for (i in c(4, 6, 8)){
  umaptest <- RunUMAP(epi_pca, dims = 1:i, verbose = F)
  print(DimPlot(umaptest, reduction = "umap", group.by = "patient_id", cols = use_colors) + labs(title = paste0(i, " dimensions")) + coord_fixed())
  ggsave2(paste0("SuppFig3E_", i, "dimensions_patient.png"), path = "../results", width = 15, height = 10, units = "cm")
  print(DimPlot(umaptest, reduction = "umap", group.by = "histo_subtype", cols = use_colors) + labs(title = paste0(i, " dimensions")) + coord_fixed())
  ggsave2(paste0("SuppFig3E_", i, "dimensions_histo.png"), path = "../results", width = 15, height = 10, units = "cm")
  print(FeaturePlot(umaptest, features = "PC_1", order = T, cols = viridis(10)) + labs(title = paste0(i, " dimensions")) + coord_fixed())
  ggsave2(paste0("SuppFig3E_", i, "dimensions_PC1.png"), path = "../results", width = 15, height = 10, units = "cm")
  remove(umaptest)
}

### bin cells along differentation gradients

# tumor cell signatures "alveolar/club-like" = tumor_signature, "undifferentiated" = anti_tumor_signature

epi_pca <- AddModuleScore(epi_pca, features = list(c("SCGB3A1", "PIGR", "NAPSA", "C4BPA", "SCGB3A2", "HLA-DRA", "CD74", "ADGRF5", "C16orf89", "FOLR1", "SELENBP1", "HLA-DRB1", "ID4", "MGP", "AQP3", "CA2", "LHX9", "HLA-DPB1", "FMO5", "GKN2", "C5", "MUC1", "NPC2", "RNASE1", "PIK3C2G", "SFTA2", "SLC34A2", "HLA-DPA1", "FGFR3", "PGC")), name = "tumor_signature")
epi_pca <- AddModuleScore(epi_pca, features = list(c("DSG2", "CAMK2N1", "FAM3C", "KRT7", "IFI27", "SLC2A1", "MARCKS", "PLAU", "AHNAK2", "PERP", "S100A4", "KRT19", "COL6A1", "UACA", "COL17A1", "CDA", "TPM2", "S100A16", "KRT8", "PRSS23", "DST", "LAMC2", "S100P", "PRSS3", "LAMA3", "DSP", "ITGA3", "MDK", "FAM83A", "ITGB4")), name = "anti_tumor_signature")

epi_tumor_data <- FetchData(epi_pca, c(colnames(epi_pca@meta.data), paste0("PC_", 1:10)))

progeny_scores_data <- as.data.frame(t(GetAssayData(epi_tumor, assay = "progeny", slot = "scale.data")))
progeny_scores_data$cell_id <- rownames(progeny_scores_data)

epi_tumor_data <- full_join(epi_tumor_data, progeny_scores_data, by = "cell_id")

progeny_names <- rownames(epi_pca@assays$progeny)

epi_tumor_data$bin <- cut_number(epi_tumor_data$PC_1, n = 10, labels = c(1:10))

ggplot(epi_tumor_data, mapping = aes(x = bin)) +
  geom_bar()

# patients along PC1

ggplot(epi_tumor_data, mapping = aes(x = bin, fill = factor(patient_id, levels = c("p032", "p018", "p019", "p024", "p031", "p030", "p033", "p023", "p027", "p034")))) +
  geom_bar() +
  theme(legend.title = element_blank()) +
  scale_fill_manual(values = use_colors)
ggsave2("SuppFig3G.pdf", path = "../results", width = 30, height = 30, units = "cm")

# histological subtypes along PC1

ggplot(epi_tumor_data, mapping = aes(x = bin, fill = factor(histo_subtype, levels = c("lepidic", "acinar", "mucinuous (papillary)", "(micro)papillary", "solid", "sarcomatoid")))) +
  geom_bar() +
  theme(legend.title = element_blank()) +
  scale_fill_manual(values = use_colors)
ggsave2("Fig2F.pdf", path = "../results", width = 30, height = 30, units = "cm")

# cell cycle phase along PC1

ggplot(epi_tumor_data, mapping = aes(x = bin, fill = Phase)) +
  geom_bar() +
  scale_fill_manual(values = use_colors)
ggsave2("SuppFig3H.pdf", path = "../results", width = 30, height = 30, units = "cm")

# normal cell type signatures along PC1

epi_tumor_data %>%
  select(c(bin, VB_Basal.11, VB_Basal.22, VB_Ciliated.13, VB_Ciliated.24, VB_Club5, VB_Goblet.26, VB_Goblet.17, VB_Ionocytes8, VB_Type.2.alveolar9, VB_Type.1.alveolar10)) %>%
  pivot_longer(cols = c(VB_Basal.11, VB_Basal.22, VB_Ciliated.13, VB_Ciliated.24, VB_Club5, VB_Goblet.26, VB_Goblet.17, VB_Ionocytes8, VB_Type.2.alveolar9, VB_Type.1.alveolar10), names_to = "cell_type") %>%
  group_by(bin, cell_type) %>%
  summarise(mean = mean(value), sd = sd(value)) %>%
  ggplot(aes(x = as.numeric(bin), y = mean, color = cell_type)) +
  geom_line(size = 1) +
  scale_color_brewer(type = "qual", palette = "Paired") +
  scale_x_continuous(breaks = c(1:20))
ggsave2("Fig2G.pdf", path = "../results", width = 30, height = 30, units = "cm")

# progeny pathway scores along PC1

getcolors <- colorRampPalette(brewer.pal(14, "Dark2"))

epi_tumor_data %>%
  select(bin, progeny_names) %>%
  pivot_longer(cols = progeny_names, names_to = "progeny") %>%
  group_by(bin, progeny) %>%
  summarise(mean = mean(value), sd = sd(value)) %>%
  ggplot(aes(x = as.numeric(bin), y = mean, color = progeny)) +
  geom_line(size = 1) +
  scale_color_manual(values = getcolors(14)) +
  scale_x_continuous(breaks = c(1:20))
ggsave2("Fig2H.pdf", path = "../results", width = 30, height = 30, units = "cm")

# tumor cell signatures along PC1

epi_tumor_data %>%
  select(bin, tumor_signature1, anti_tumor_signature1) %>%
  pivot_longer(cols = c(tumor_signature1, anti_tumor_signature1), names_to = "genes") %>%
  group_by(bin, genes) %>%
  summarise(mean = mean(value), sd = sd(value)) %>%
  ggplot(aes(x = as.numeric(bin), y = mean, color = genes)) +
  geom_line(size = 1) +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width = 0.3) +
  scale_x_continuous(breaks = c(1:20))
ggsave2("SuppFig3I.pdf", path = "../results", width = 30, height = 30, units = "cm")

#######################################stromal analyses

### load libraries
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(ggplot2)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(readr)
library(stringr)
library(gplots)
library(grid)
library(rlang)
library(tibble)

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467",
  Endothelial1 = "#FED976",
  Endothelial2 = "#FEB24C",
  Endothelial3 = "#fd8d3C",
  Endothelial4 = "#FC4E2A",
  Endothelial5 = "#E31A1C",
  Endothelial6 = "#BD0026",
  Endothelial7 = "#800026",
  Lymphaticendothelial = "salmon",
  Fibroblast1 = "#2166AC",
  Fibroblast2 = "#4393C3",
  Myofibroblast1 = "#5AAE61",
  Myofibroblast2 = "#1B7837",
  Smoothmuscle1 = "#9970AB",
  Smoothmuscle2 = "#762A83",
  Mesothelial = "#40004B")

#str_anno <- readRDS("seurat_objects/str_anno.RDS")

str_anno@meta.data$cell_type_str <- factor(str_anno@meta.data$cell_type_str, levels = c("Endothelial1",
                                                                                        "Endothelial2",
                                                                                        "Endothelial3",
                                                                                        "Endothelial4",
                                                                                        "Endothelial5",
                                                                                        "Endothelial6",
                                                                                        "Endothelial7",
                                                                                        "Lymphaticendothelial",
                                                                                        "Fibroblast1",
                                                                                        "Fibroblast2",
                                                                                        "Myofibroblast1",
                                                                                        "Myofibroblast2",
                                                                                        "Smoothmuscle1",
                                                                                        "Smoothmuscle2",
                                                                                        "Mesothelial"))

DimPlot(str_anno, group.by = "tissue_type", cols = use_colors)
#ggsave2("DimPlot_str_Normal_Tumor.pdf", path = "output/fig3", width = 15, height = 15, units = "cm")

DimPlot(str_anno, group.by = "patient_id", cols = use_colors, pt.size = 0.5)
ggsave2("SuppFig1C_str_patients.pdf", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(str_anno, group.by = "cell_type_str", label = F, split.by = "tissue_type", cols = use_colors, pt.size = 0.5)
ggsave2("Fig3A_umap.pdf", path = "../results", width = 30, height = 15, units = "cm")

DotPlot(str_anno, features = c("WT1", "UPK3B", "MYH11", "PDGFRB", "ACTA2", "MYLK", "LUM", "PDGFRA", "CCL21", "PROX1", "PECAM1", "VWF"), group.by = "cell_type_str") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  coord_flip() + 
  scale_color_viridis()
ggsave2("Fig3B.pdf", path = "../results", width = 16, height = 12, units = "cm")

###subsetting

str_endo <- subset(str_anno, subset = cell_type_str %in% c("Endothelial1",
                                                           "Endothelial2",
                                                           "Endothelial3",
                                                           "Endothelial4",
                                                           "Endothelial5",
                                                           "Endothelial6",
                                                           "Endothelial7",
                                                           "Lymphaticendothelial"))
str_endo <- ScaleData(str_endo)

str_fibro <- subset(str_anno, subset = cell_type_str %in% c("Fibroblast1",
                                                            "Fibroblast2",
                                                            "Myofibroblast1",
                                                            "Myofibroblast2",
                                                            "Smoothmuscle1",
                                                            "Smoothmuscle2",
                                                            "Mesothelial"))

str_fibro <- ScaleData(str_fibro)

endo_counts <- FetchData(str_endo, vars = c("tissue_type", "cell_type_str", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal")))

endo_counts_tbl <- endo_counts %>%
  dplyr::count(cell_type_str, patient_id, tissue_type)
write_csv(endo_counts_tbl, path = "../results/SuppTable1.csv")

fibro_counts <- FetchData(str_fibro, vars = c("tissue_type", "cell_type_str", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal"))) 

fibro_counts_tbl <- fibro_counts %>%
  dplyr::count(cell_type_str, patient_id, tissue_type)
write_csv(fibro_counts_tbl, path = "../results/SuppTable2.csv")

ggplot(data = endo_counts, aes(x = tissue_type, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig3A_barplot_endothelial.pdf", path = "../results", width = 20, height = 5, units = "cm")

ggplot(data = fibro_counts, aes(x = tissue_type, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig3A_barplot_fibroblastic.pdf", path = "../results", width = 20, height = 5, units = "cm")

endo_counts %>%
  filter(tissue_type == "Tumor") %>%
  ggplot(aes(x = sample_id, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig3A_barplot_endothelial_per_patient.pdf", path = "../results", width = 30, height = 30, units = "cm")

fibro_counts %>%
  filter(tissue_type == "Tumor") %>%
  ggplot(aes(x = sample_id, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig3A_barplot_fibroblastic_per_patient.pdf", path = "../results", width = 30, height = 30, units = "cm")

endo_counts %>%
  filter(tissue_type == "Normal") %>%
  ggplot(aes(x = sample_id, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("SuppFig6A_endothelial.pdf", path = "../results", width = 30, height = 30, units = "cm")

fibro_counts %>%
  filter(tissue_type == "Normal") %>%
  ggplot(aes(x = sample_id, fill = cell_type_str)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("SuppFig6A_fibroblastic.pdf", path = "../results", width = 30, height = 30, units = "cm")

#heatmap wrapper function
###code from https://github.com/satijalab/seurat/issues/2201

DoMultiBarHeatmap <- function (object, 
                               features = NULL, 
                               cells = NULL, 
                               group.by = "ident", 
                               additional.group.by = NULL, 
                               additional.group.sort.by = NULL, 
                               cols.use = NULL,
                               group.bar = TRUE, 
                               disp.min = -2.5, 
                               disp.max = NULL, 
                               slot = "scale.data", 
                               assay = NULL, 
                               label = TRUE, 
                               size = 5.5, 
                               hjust = 0, 
                               angle = 45, 
                               raster = TRUE, 
                               draw.lines = TRUE, 
                               lines.width = NULL, 
                               group.bar.height = 0.02, 
                               combine = TRUE) 
{
  cells <- cells %||% colnames(x = object)
  if (is.numeric(x = cells)) {
    cells <- colnames(x = object)[cells]
  }
  assay <- assay %||% DefaultAssay(object = object)
  DefaultAssay(object = object) <- assay
  features <- features %||% VariableFeatures(object = object)
  ## Why reverse???
  features <- rev(x = unique(x = features))
  disp.max <- disp.max %||% ifelse(test = slot == "scale.data", 
                                   yes = 2.5, no = 6)
  possible.features <- rownames(x = GetAssayData(object = object, 
                                                 slot = slot))
  if (any(!features %in% possible.features)) {
    bad.features <- features[!features %in% possible.features]
    features <- features[features %in% possible.features]
    if (length(x = features) == 0) {
      stop("No requested features found in the ", slot, 
           " slot for the ", assay, " assay.")
    }
    warning("The following features were omitted as they were not found in the ", 
            slot, " slot for the ", assay, " assay: ", paste(bad.features, 
                                                             collapse = ", "))
  }

  if (!is.null(additional.group.sort.by)) {
    if (any(!additional.group.sort.by %in% additional.group.by)) {
      bad.sorts <- additional.group.sort.by[!additional.group.sort.by %in% additional.group.by]
      additional.group.sort.by <- additional.group.sort.by[additional.group.sort.by %in% additional.group.by]
      if (length(x = bad.sorts) > 0) {
        warning("The following additional sorts were omitted as they were not a subset of additional.group.by : ", 
                paste(bad.sorts, collapse = ", "))
      }
    }
  }

  data <- as.data.frame(x = as.matrix(x = t(x = GetAssayData(object = object, 
                                                             slot = slot)[features, cells, drop = FALSE])))

  object <- suppressMessages(expr = StashIdent(object = object, 
                                               save.name = "ident"))
  group.by <- group.by %||% "ident"
  groups.use <- object[[c(group.by, additional.group.by[!additional.group.by %in% group.by])]][cells, , drop = FALSE]
  plots <- list()
  for (i in group.by) {
    data.group <- data
    if (!is_null(additional.group.by)) {
      additional.group.use <- additional.group.by[additional.group.by!=i]  
      if (!is_null(additional.group.sort.by)){
        additional.sort.use = additional.group.sort.by[additional.group.sort.by != i]  
      } else {
        additional.sort.use = NULL
      }
    } else {
      additional.group.use = NULL
      additional.sort.use = NULL
    }

    group.use <- groups.use[, c(i, additional.group.use), drop = FALSE]

    for(colname in colnames(group.use)){
      if (!is.factor(x = group.use[[colname]])) {
        group.use[[colname]] <- factor(x = group.use[[colname]])
      }  
    }

    if (draw.lines) {
      lines.width <- lines.width %||% ceiling(x = nrow(x = data.group) * 
                                                0.0025)
      placeholder.cells <- sapply(X = 1:(length(x = levels(x = group.use[[i]])) * 
                                           lines.width), FUN = function(x) {
                                             return(Seurat:::RandomName(length = 20))
                                           })
      placeholder.groups <- data.frame(rep(x = levels(x = group.use[[i]]), times = lines.width))
      group.levels <- list()
      group.levels[[i]] = levels(x = group.use[[i]])
      for (j in additional.group.use) {
        group.levels[[j]] <- levels(x = group.use[[j]])
        placeholder.groups[[j]] = NA
      }

      colnames(placeholder.groups) <- colnames(group.use)
      rownames(placeholder.groups) <- placeholder.cells

      group.use <- sapply(group.use, as.vector)
      rownames(x = group.use) <- cells

      group.use <- rbind(group.use, placeholder.groups)

      for (j in names(group.levels)) {
        group.use[[j]] <- factor(x = group.use[[j]], levels = group.levels[[j]])
      }

      na.data.group <- matrix(data = NA, nrow = length(x = placeholder.cells), 
                              ncol = ncol(x = data.group), dimnames = list(placeholder.cells, 
                                                                           colnames(x = data.group)))
      data.group <- rbind(data.group, na.data.group)
    }

    order_expr <- paste0('order(', paste(c(i, additional.sort.use), collapse=','), ')')
    group.use = with(group.use, group.use[eval(parse(text=order_expr)), , drop=F])

    plot <- Seurat:::SingleRasterMap(data = data.group, raster = raster, 
                                     disp.min = disp.min, disp.max = disp.max, feature.order = features, 
                                     cell.order = rownames(x = group.use), group.by = group.use[[i]])

    if (group.bar) {
      pbuild <- ggplot_build(plot = plot)
      group.use2 <- group.use
      cols <- list()
      na.group <- Seurat:::RandomName(length = 20)
      for (colname in rev(x = colnames(group.use2))) {
        if (colname == i) {
          colid = paste0('Identity (', colname, ')')
        } else {
          colid = colname
        }

        # Default
        cols[[colname]] <- c(scales::hue_pal()(length(x = levels(x = group.use[[colname]]))))  

        #Overwrite if better value is provided
        if (!is_null(cols.use[[colname]])) {
          req_length = length(x = levels(group.use))
          if (length(cols.use[[colname]]) < req_length){
            warning("Cannot use provided colors for ", colname, " since there aren't enough colors.")
          } else {
            if (!is_null(names(cols.use[[colname]]))) {
              if (all(levels(group.use[[colname]]) %in% names(cols.use[[colname]]))) {
                cols[[colname]] <- as.vector(cols.use[[colname]][levels(group.use[[colname]])])
              } else {
                warning("Cannot use provided colors for ", colname, " since all levels (", paste(levels(group.use[[colname]]), collapse=","), ") are not represented.")
              }
            } else {
              cols[[colname]] <- as.vector(cols.use[[colname]])[c(1:length(x = levels(x = group.use[[colname]])))]
            }
          }
        }

        # Add white if there's lines
        if (draw.lines) {
          levels(x = group.use2[[colname]]) <- c(levels(x = group.use2[[colname]]), na.group)  
          group.use2[placeholder.cells, colname] <- na.group
          cols[[colname]] <- c(cols[[colname]], "#FFFFFF")
        }
        names(x = cols[[colname]]) <- levels(x = group.use2[[colname]])

        y.range <- diff(x = pbuild$layout$panel_params[[1]]$y.range)
        y.pos <- max(pbuild$layout$panel_params[[1]]$y.range) + y.range * 0.015
        y.max <- y.pos + group.bar.height * y.range
        pbuild$layout$panel_params[[1]]$y.range <- c(pbuild$layout$panel_params[[1]]$y.range[1], y.max)

        plot <- suppressMessages(plot + 
                                   annotation_raster(raster = t(x = cols[[colname]][group.use2[[colname]]]),  xmin = -Inf, xmax = Inf, ymin = y.pos, ymax = y.max) + 
                                   annotation_custom(grob = grid::textGrob(label = colid, hjust = 0, gp = gpar(cex = 0.75)), ymin = mean(c(y.pos, y.max)), ymax = mean(c(y.pos, y.max)), xmin = Inf, xmax = Inf) +
                                   coord_cartesian(ylim = c(0, y.max), clip = "off")) 

        if ((colname == i) && label) {
          x.max <- max(pbuild$layout$panel_params[[1]]$x.range)
          x.divs <- pbuild$layout$panel_params[[1]]$x.major %||% pbuild$layout$panel_params[[1]]$x$break_positions()
          group.use$x <- x.divs
          label.x.pos <- tapply(X = group.use$x, INDEX = group.use[[colname]],
                                FUN = median) * x.max
          label.x.pos <- data.frame(group = names(x = label.x.pos), 
                                    label.x.pos)
          plot <- plot + geom_text(stat = "identity", 
                                   data = label.x.pos, aes_string(label = "group", 
                                                                  x = "label.x.pos"), y = y.max + y.max * 
                                     0.03 * 0.5, angle = angle, hjust = hjust, 
                                   size = size)
          plot <- suppressMessages(plot + coord_cartesian(ylim = c(0, 
                                                                   y.max + y.max * 0.002 * max(nchar(x = levels(x = group.use[[colname]]))) * 
                                                                     size), clip = "off"))
        }
      }
    }
    plot <- plot + theme(line = element_blank())
    plots[[i]] <- plot
  }
  if (combine) {
    plots <- CombinePlots(plots = plots)
  }
  return(plots)
}

###DEGs endothelial

Idents(str_endo) <- str_endo@meta.data$cell_type_str

endo_markers <- FindAllMarkers(str_endo, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)

top_endo_markers <- endo_markers %>% group_by(cluster) %>% top_n(10, wt = avg_logFC)

DoMultiBarHeatmap(str_endo, features = top_endo_markers$gene, group.by = "cell_type_str", additional.group.by = "tissue_type",additional.group.sort.by = "tissue_type", cols.use = list(tissue_type = use_colors), draw.lines = F) +
  scale_fill_viridis()
ggsave2("SuppFig6B.png", path = "../results", width = 30, height = 40, units = "cm")
#ggsave2("HeatMap_Endo.pdf", path = "output/fig3", width = 30, height = 40, units = "cm")
#ggsave2("HeatMap_Endo.emf", path = "output/fig3", width = 30, height = 40, units = "cm")

###DEGs fibroblastic

Idents(str_fibro) <- str_fibro@meta.data$cell_type_str

fibro_markers <- FindAllMarkers(str_fibro, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)

top_fibro_markers <- fibro_markers %>% group_by(cluster) %>% top_n(10, wt = avg_logFC)

DoMultiBarHeatmap(str_fibro, features = top_fibro_markers$gene, group.by = "cell_type_str", additional.group.by = "tissue_type",additional.group.sort.by = "tissue_type", cols.use = list(tissue_type = use_colors), draw.lines = F) +
  scale_fill_viridis()
ggsave2("Fig6C.png", path = "../results", width = 30, height = 40, units = "cm")
#ggsave2("HeatMap_Fibro.pdf", path = "output/fig3", width = 30, height = 40, units = "cm")
#ggsave2("HeatMap_Fibro.emf", path = "output/fig3", width = 30, height = 40, units = "cm")

###progeny scores

str_fibro2 <- subset(str_anno, subset = cell_type_str %in% c("Fibroblast1",
                                                             "Fibroblast2",
                                                             "Myofibroblast1",
                                                             "Myofibroblast2",
                                                             "Smoothmuscle1",
                                                             "Smoothmuscle2"))

progeny_scores <- as.data.frame(t(GetAssayData(str_fibro2, assay = "progeny", slot = "scale.data")))
progeny_scores$cell_id <- rownames(progeny_scores)
progeny_scores <- gather(progeny_scores, Pathway, Activity, -cell_id)

cells_clusters <- FetchData(str_anno, c("cell_type_str"))
cells_clusters$cell_id <- rownames(cells_clusters)

progeny_scores <- inner_join(progeny_scores, cells_clusters)

summarized_progeny_scores <- progeny_scores %>% 
  group_by(Pathway, cell_type_str) %>% 
  summarise(avg = mean(Activity), std = sd(Activity)) %>%
  pivot_wider(id_cols = Pathway, names_from = cell_type_str, values_from = avg) %>%
  column_to_rownames("Pathway") %>%
  as.matrix()

pdf("../results/Fig3D.pdf", width = 7, height = 10)
heatmap.2(summarized_progeny_scores, trace = "none", density.info = "none", col = bluered(100), margins = c(10,10))
dev.off()

#############################################immune analyses

### load libraries
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(ggplot2)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(stringr)
library(cowplot)
library(scales)
library(readr)
library(progeny)
library(gplots)
library(tibble)
library(grid)
library(rlang)

theme_set(theme_cowplot())

use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467",
  Alveolar_Macrophages1 = "#6bAEd6",
  Alveolar_Macrophages2 = "#3182BD",
  Alveolar_Macrophages3 = "#08519C",
  CD14_Macrophages1= "#fff500",
  CD14_Macrophages2= "#FE9929",
  CD14_Macrophages3= "#EC7014",
  CD14_Macrophages4= "#CC4C02",
  CD14_Macrophages5= "#8C2D04",
  Macrophages_Proliferating= "#E31A1C",
  Monocytes= "#FA9FB5",
  Myeloid_Dendritic= "#DD3497",
  Plasmacytoid_Dendritic= "#7A0177",
  T_conv1= "#c2e699",
  T_conv2= "#78c679",
  T_reg= "#006837",
  T_CD8_1= "#bcbddc",
  T_CD8_2= "#9e9ac8",
  T_CD8_3= "#807dba",
  T_CD8_Proliferating= "#6a51a3",
  NK_cells= "#4a1486",
  B_cells= "#969696",
  Plasma= "#636363",
  Mast= "#252525")

#load data

#imm_anno <- readRDS("seurat_objects/imm_anno.RDS")

imm_anno@meta.data$cell_type_imm <- ordered(imm_anno@meta.data$cell_type_imm, levels = c("Alveolar_Macrophages1",
                                                                                         "Alveolar_Macrophages2",
                                                                                         "Alveolar_Macrophages3",
                                                                                         "CD14_Macrophages1",
                                                                                         "CD14_Macrophages2",
                                                                                         "CD14_Macrophages3",
                                                                                         "CD14_Macrophages4",
                                                                                         "CD14_Macrophages5",
                                                                                         "Macrophages_Proliferating",
                                                                                         "Monocytes",
                                                                                         "Myeloid_Dendritic",
                                                                                         "Plasmacytoid_Dendritic",
                                                                                         "Mast",
                                                                                         "T_conv1",
                                                                                         "T_conv2",
                                                                                         "T_reg",
                                                                                         "T_CD8_1",
                                                                                         "T_CD8_2",
                                                                                         "T_CD8_3",
                                                                                         "T_CD8_Proliferating",
                                                                                         "NK_cells",
                                                                                         "B_cells",
                                                                                         "Plasma"))

DimPlot(imm_anno, group.by = "tissue_type", cols = use_colors)
#ggsave2("DimPlot_imm_Normal_Tumor.pdf", path = "output/fig4", width = 15, height = 15, units = "cm")
#ggsave2("DimPlot_imm_Normal_Tumor.png", path = "output/fig4", width = 35, height = 15, units = "cm")

DimPlot(imm_anno, group.by = "patient_id", cols = use_colors, pt.size = 0.5)
#ggsave2("DimPlot_imm_patients.pdf", path = "output/fig4", width = 30, height = 15, units = "cm")
ggsave2("SuppFig1C_imm_patients.png", path = "../results", width = 15, height = 15, units = "cm")

DimPlot(imm_anno, group.by = "cell_type_imm", split.by = "tissue_type", cols = use_colors, pt.size = 0.5)
#ggsave2("DimPlot_imm_celltype.pdf", path = "output/fig4", width = 35, height = 15, units = "cm")
ggsave2("Fig4A_umap.png", path = "../results", width = 35, height = 15, units = "cm")

DotPlot(imm_anno, features = c("CD68", "LYZ", "FABP4", "MARCO", "LGMN", "CSF1R", "CD14", "S100A12", "FCN1", "CD1C", "FCER1A", "LILRA4", "IL3RA", "KIT", "GATA2", "CD3E", "CD4", "FOXP3", "IL2RA", "CD8A", "NKG7", "KLRD1", "MS4A1", "CD79A", "JCHAIN", "IGKC", "MKI67"), group.by = "cell_type_imm") + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  coord_flip() + 
  scale_color_viridis()
ggsave2("Fig4B.pdf", path = "../results", width = 20, height = 20, units = "cm")

###subsetting

imm_lympho <- subset(imm_anno, subset = cell_type_imm %in% c("T_conv1",
                                                             "T_conv2",
                                                             "T_reg",
                                                             "T_CD8_1",
                                                             "T_CD8_2",
                                                             "T_CD8_3",
                                                             "T_CD8_Proliferating",
                                                             "NK_cells",
                                                             "B_cells",
                                                             "Plasma"))

imm_lympho <- ScaleData(imm_lympho)

imm_myelo <- subset(imm_anno, subset = cell_type_imm %in% c("Alveolar_Macrophages1",
                                                            "Alveolar_Macrophages2",
                                                            "Alveolar_Macrophages3",
                                                            "CD14_Macrophages1",
                                                            "CD14_Macrophages2",
                                                            "CD14_Macrophages3",
                                                            "CD14_Macrophages4",
                                                            "CD14_Macrophages5",
                                                            "Macrophages_Proliferating",
                                                            "Monocytes",
                                                            "Myeloid_Dendritic",
                                                            "Plasmacytoid_Dendritic",
                                                            "Mast"))

imm_myelo <- ScaleData(imm_myelo)

lympho_counts <- FetchData(imm_lympho, vars = c("tissue_type", "cell_type_imm", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal")))

lympho_counts_tbl <- lympho_counts %>%
  dplyr::count(cell_type_imm, patient_id, tissue_type)
write_csv(lympho_counts_tbl, path = "../results/SuppTable4.csv")

myelo_counts <- FetchData(imm_myelo, vars = c("tissue_type", "cell_type_imm", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal"))) 

myelo_counts_tbl <- myelo_counts %>%
  dplyr::count(cell_type_imm, patient_id, tissue_type)
write_csv(myelo_counts_tbl, path = "../results/SuppTable3.csv")

ggplot(data = lympho_counts, aes(x = tissue_type, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig4A_barplot_lymphoid.pdf", path = "../results", width = 20, height = 5, units = "cm")

ggplot(data = myelo_counts, aes(x = tissue_type, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig4A_barplot_myeloid.pdf", path = "../results", width = 20, height = 5, units = "cm")

lympho_counts %>%
  filter(tissue_type == "Tumor") %>%
  ggplot(aes(x = sample_id, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig4A_barplot_lymphoid_per_patient.pdf", path = "../results", width = 30, height = 30, units = "cm")

myelo_counts %>%
  filter(tissue_type == "Tumor") %>%
  ggplot(aes(x = sample_id, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("Fig4A_barplot_myeloid_per_patient.pdf", path = "../results", width = 30, height = 30, units = "cm")

lympho_counts %>%
  filter(tissue_type == "Normal") %>%
  ggplot(aes(x = sample_id, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("SuppFig7A_lymphoid.pdf", path = "../results", width = 30, height = 30, units = "cm")

myelo_counts %>%
  filter(tissue_type == "Normal") %>%
  ggplot(aes(x = sample_id, fill = cell_type_imm)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = use_colors) +
  coord_flip() +
  scale_y_reverse()
ggsave2("SuppFig7A_myeloid.pdf", path = "../results", width = 30, height = 30, units = "cm")

###DEGs macrophages

imm_macro <- subset(imm_anno, subset = cell_type_imm %in% c("Alveolar_Macrophages1",
                                                            "Alveolar_Macrophages2",
                                                            "Alveolar_Macrophages3",
                                                            "CD14_Macrophages1",
                                                            "CD14_Macrophages2",
                                                            "CD14_Macrophages3",
                                                            "CD14_Macrophages4",
                                                            "CD14_Macrophages5",
                                                            "Macrophages_Proliferating"))

imm_macro <- ScaleData(imm_macro)

Idents(imm_macro) <- imm_macro@meta.data$cell_type_imm

macro_markers <- FindAllMarkers(imm_macro, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)

top_macro_markers <- macro_markers %>% group_by(cluster) %>% top_n(10, wt = avg_logFC)

DoMultiBarHeatmap(imm_macro, features = top_macro_markers$gene, group.by = "cell_type_imm", additional.group.by = "tissue_type",additional.group.sort.by = "tissue_type", cols.use = list(tissue_type = use_colors), draw.lines = F) +
  scale_fill_viridis()
#ggsave2("HeatMap_Macro.pdf", path = "output/fig4", width = 30, height = 40, units = "cm")
ggsave2("SuppFig7B.png", path = "../results", width = 30, height = 40, units = "cm")
#ggsave2("HeatMap_Macro.emf", path = "output/fig4", width = 30, height = 40, units = "cm")

###DEGs T cells

imm_T <- subset(imm_anno, subset = cell_type_imm %in% c("T_conv1",
                                                        "T_conv2",
                                                        "T_reg",
                                                        "T_CD8_1",
                                                        "T_CD8_2",
                                                        "T_CD8_3",
                                                        "T_CD8_Proliferating"))

imm_T <- ScaleData(imm_T)

Idents(imm_T) <- imm_T@meta.data$cell_type_imm

markers_T <- FindAllMarkers(imm_T, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)

top_markers_T <- markers_T %>% group_by(cluster) %>% top_n(10, wt = avg_logFC)

DoMultiBarHeatmap(imm_T, features = top_markers_T$gene, group.by = "cell_type_imm", additional.group.by = "tissue_type",additional.group.sort.by = "tissue_type", cols.use = list(tissue_type = use_colors), draw.lines = F) +
  scale_fill_viridis()
#ggsave2("HeatMap_T.pdf", path = "output/fig4", width = 30, height = 35, units = "cm")
ggsave2("SuppFig7C.png", path = "../results", width = 30, height = 35, units = "cm")
#ggsave2("HeatMap_T.emf", path = "output/fig4", width = 30, height = 35, units = "cm")

###selected hallmark signatures and M1vsM2 signatures

m1m2_pws <- read_lines("../data/data/CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_UP.gmt") %>%
  lapply(str_split, "\\t") %>% 
  unlist(recursive = F) %>% 
  lapply(function(x) setNames(list(x[-c(1:2)]), x[1])) %>% 
  unlist(recursive = F)

m1m2_pws <- append(m1m2_pws, read_lines("../data/data/CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_DN.gmt") %>%
                     lapply(str_split, "\\t") %>% 
                     unlist(recursive = F) %>% 
                     lapply(function(x) setNames(list(x[-c(1:2)]), x[1])) %>% 
                     unlist(recursive = F))

imm_anno <- AddModuleScore(object = imm_anno, features = m1m2_pws, name = c("m1up", "m1dn"), nbin = 12)

VlnPlot(imm_anno, features = c("HALLMARK_INFLAMMATORY_RESPONSE31",
                               "HALLMARK_ALLOGRAFT_REJECTION46",
                               "HALLMARK_INTERFERON_GAMMA_RESPONSE19",
                               "HALLMARK_TNFA_SIGNALING_VIA_NFKB1",
                               "m1up1", 
                               "m1dn2"), 
        group.by = "cell_type_imm", pt.size = 0, ncol = 3, idents = c("Alveolar_Macrophages1",
                                                                      "Alveolar_Macrophages2",
                                                                      "Alveolar_Macrophages3",
                                                                      "CD14_Macrophages1",
                                                                      "CD14_Macrophages2",
                                                                      "CD14_Macrophages3",
                                                                      "CD14_Macrophages4",
                                                                      "CD14_Macrophages5",
                                                                      "Macrophages_Proliferating"), cols = use_colors)
ggsave2("Fig4C.pdf", path = "../results", width = 30, height = 20, units = "cm")

###T cell signatures

T_exhausted <- read_excel("../data/data/CD8_T_cells_exhausted.xlsx", skip = 1)
cytotoxicity <- c("PRF1", "IFNG", "GNLY", "NKG7", "GZMB", "GZMA", "GZMH", "KLRK1", "KLRB1", "KLRD1", "CTSW", "CST7")

T_cell_markers <- list(T_exhausted$GeneSymbol, cytotoxicity)

imm_T <- AddModuleScore(imm_T, features = T_cell_markers, name = c("exhaustion", "cytotoxicity"), nbin = 12)

VlnPlot(imm_T, features = c("cytotoxicity2", "exhaustion1"), pt.size = 0, group.by = "cell_type_imm", cols = use_colors, idents = c("T_CD8_1", "T_CD8_2", "T_CD8_3", "T_CD8_Proliferating"), ncol = 1)
ggsave2("Fig4F.pdf", path = "../results", width = 10, height = 20, units = "cm")

#############################################correlation cell types

### load libraries
library(ggplot2)
library(Seurat)
library(dplyr)
library(tidyr)
library(cowplot)
library(gplots)
library(RColorBrewer)
library(plotrix)
library(corrplot)
library(pcaMethods)
library(readr)

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467",
  lepidic = "#0B775E",
  acinar = "#74A089",
  mucinuous (papillary) = "#E2D200",
  (micro)papillary = "#CEAB07",
  solid = "#B40F20",
  sarcomatoid = "#5B1A18")

###load data and subsetting

#epi_anno <- readRDS("seurat_objects/epi_anno.RDS")

epi_tumor <- SubsetData(SubsetData(epi_anno, subset.name = "cluster_type", accept.value = "Tumor"), subset.name = "tissue_type", accept.value = "Tumor")
epi_tumor <- ScaleData(epi_tumor)

epi_pca <- epi_tumor
epi_pca <- RunPCA(epi_pca)

#imm_anno <- readRDS("seurat_objects/imm_anno.RDS")

imm_anno@meta.data$cell_type_imm <- ordered(imm_anno@meta.data$cell_type_imm, levels = c("Alveolar_Macrophages1",
                                                                                         "Alveolar_Macrophages2",
                                                                                         "Alveolar_Macrophages3",
                                                                                         "CD14_Macrophages1",
                                                                                         "CD14_Macrophages2",
                                                                                         "CD14_Macrophages3",
                                                                                         "CD14_Macrophages4",
                                                                                         "CD14_Macrophages5",
                                                                                         "Macrophages_Proliferating",
                                                                                         "Monocytes",
                                                                                         "Myeloid_Dendritic",
                                                                                         "Plasmacytoid_Dendritic",
                                                                                         "Mast",
                                                                                         "T_conv1",
                                                                                         "T_conv2",
                                                                                         "T_reg",
                                                                                         "T_CD8_1",
                                                                                         "T_CD8_2",
                                                                                         "T_CD8_3",
                                                                                         "T_CD8_Proliferating",
                                                                                         "NK_cells",
                                                                                         "B_cells",
                                                                                         "Plasma"))

imm_lympho <- subset(imm_anno, subset = cell_type_imm %in% c("T_conv1",
                                                             "T_conv2",
                                                             "T_reg",
                                                             "T_CD8_1",
                                                             "T_CD8_2",
                                                             "T_CD8_3",
                                                             "T_CD8_Proliferating",
                                                             "NK_cells",
                                                             "B_cells",
                                                             "Plasma"))

imm_myelo <- subset(imm_anno, subset = cell_type_imm %in% c("Alveolar_Macrophages1",
                                                            "Alveolar_Macrophages2",
                                                            "Alveolar_Macrophages3",
                                                            "CD14_Macrophages1",
                                                            "CD14_Macrophages2",
                                                            "CD14_Macrophages3",
                                                            "CD14_Macrophages4",
                                                            "CD14_Macrophages5",
                                                            "Macrophages_Proliferating",
                                                            "Monocytes",
                                                            "Myeloid_Dendritic",
                                                            "Plasmacytoid_Dendritic",
                                                            "Mast"))

lympho_counts <- FetchData(imm_lympho, vars = c("tissue_type", "cell_type_imm", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal")))

myelo_counts <- FetchData(imm_myelo, vars = c("tissue_type", "cell_type_imm", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal"))) 

#str_anno <- readRDS("seurat_objects/str_anno.RDS")

str_anno@meta.data$cell_type_str <- factor(str_anno@meta.data$cell_type_str, levels = c("Endothelial1",
                                                                                        "Endothelial2",
                                                                                        "Endothelial3",
                                                                                        "Endothelial4",
                                                                                        "Endothelial5",
                                                                                        "Endothelial6",
                                                                                        "Endothelial7",
                                                                                        "Lymphaticendothelial",
                                                                                        "Fibroblast1",
                                                                                        "Fibroblast2",
                                                                                        "Myofibroblast1",
                                                                                        "Myofibroblast2",
                                                                                        "Smoothmuscle1",
                                                                                        "Smoothmuscle2",
                                                                                        "Mesothelial"))

str_endo <- subset(str_anno, subset = cell_type_str %in% c("Endothelial1",
                                                           "Endothelial2",
                                                           "Endothelial3",
                                                           "Endothelial4",
                                                           "Endothelial5",
                                                           "Endothelial6",
                                                           "Endothelial7",
                                                           "Lymphaticendothelial"))

str_fibro <- subset(str_anno, subset = cell_type_str %in% c("Fibroblast1",
                                                            "Fibroblast2",
                                                            "Myofibroblast1",
                                                            "Myofibroblast2",
                                                            "Smoothmuscle1",
                                                            "Smoothmuscle2",
                                                            "Mesothelial"))

endo_counts <- FetchData(str_endo, vars = c("tissue_type", "cell_type_str", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal")))

fibro_counts <- FetchData(str_fibro, vars = c("tissue_type", "cell_type_str", "sample_id", "patient_id")) %>%  
  mutate(tissue_type = factor(tissue_type, levels = c("Tumor", "Normal"))) 

###count immune and stromal cells

myelo_counts_rel <- myelo_counts %>%
  filter(tissue_type == "Tumor") %>%
  dplyr::count(cell_type_imm, patient_id) %>%
  group_by(patient_id) %>% 
  mutate(n_rel = n/sum(n))

myelo_counts_rel <- myelo_counts_rel %>%
  pivot_wider(id_cols = patient_id, names_from = cell_type_imm, values_from = n_rel)

lympho_counts_rel <- lympho_counts %>%
  filter(tissue_type == "Tumor") %>%
  dplyr::count(cell_type_imm, patient_id) %>%
  group_by(patient_id) %>% 
  mutate(n_rel = n/sum(n)) 

lympho_counts_rel <- lympho_counts_rel %>%
  pivot_wider(id_cols = patient_id, names_from = cell_type_imm, values_from = n_rel)

fibro_counts_rel <- fibro_counts %>%
  filter(tissue_type == "Tumor") %>%
  dplyr::count(cell_type_str, patient_id) %>%
  group_by(patient_id) %>% 
  mutate(n_rel = n/sum(n))

fibro_counts_rel <- fibro_counts_rel %>%
  pivot_wider(id_cols = patient_id, names_from = cell_type_str, values_from = n_rel)

endo_counts_rel <- endo_counts %>%
  filter(tissue_type == "Tumor") %>%
  dplyr::count(cell_type_str, patient_id) %>%
  group_by(patient_id) %>% 
  mutate(n_rel = n/sum(n))

endo_counts_rel <- endo_counts_rel %>%
  pivot_wider(id_cols = patient_id, names_from = cell_type_str, values_from = n_rel)

cell_counts_rel <- full_join(myelo_counts_rel, lympho_counts_rel, by = "patient_id")
cell_counts_rel <- full_join(cell_counts_rel, endo_counts_rel, by = "patient_id")
cell_counts_rel <- full_join(cell_counts_rel, fibro_counts_rel, by = "patient_id")

cell_counts_rel[is.na(cell_counts_rel)] <- 0

###cell proportion statistics

test <- cell_counts_rel %>%
  mutate(pattern = ifelse(patient_id %in% c("p032", "p018", "p019", "p024", "p031", "p033"), "N3MC", "CP2E")) %>%
  mutate(T_CD8_exhausted = T_CD8_1 + T_CD8_2)

#for Kim et al. dataset
#test <- cell_counts_rel %>%
#  mutate(pattern = ifelse(patient_id %in% c("P0009", "P0018", "P0020", "P0028"), "CP2E", "N3MC")) %>%
#  mutate(T_CD8_exhausted = T_CD8_1 + T_CD8_2)

ggplot(test, aes(x = pattern, y = CD14_Macrophages1)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = CD14_Macrophages1~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_CD14_Macrophages1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = CD14_Macrophages2)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = CD14_Macrophages2~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_CD14_Macrophages2.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = Myeloid_Dendritic)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = Myeloid_Dendritic~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_Myeloid_Dendritic.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = Plasmacytoid_Dendritic)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = Plasmacytoid_Dendritic~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_Plasmacytoid_Dendritic.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = Myofibroblast1)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = Myofibroblast1~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_Myofibroblast1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = Myofibroblast2)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = Myofibroblast2~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_Myofibroblast2.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = T_conv1)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = T_conv1~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_T_conv1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = NK_cells)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = NK_cells~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_NK_cells.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(test, aes(x = pattern, y = T_CD8_exhausted)) +
  geom_boxplot() +
  ggtitle(paste0("p = ", wilcox.test(formula = T_CD8_exhausted~pattern, data = test, alternative = "two.sided", paired = F)$p.value))
ggsave2("SuppFig8A_Group1vs2_T_CD8_exhausted.pdf", path = "../results", width = 10, height = 10, units = "cm")

###select only cell types that occured in at least three patients

celltype_selected <- cell_counts_rel
celltype_selected[celltype_selected == 0] <- NA
celltype_selected <- cell_counts_rel[, which(colMeans(!is.na(celltype_selected)) >= 0.3)]

cell_counts_rel_selected <- cell_counts_rel[, colnames(cell_counts_rel) %in% colnames(celltype_selected)]

cell_counts_rel_mtrx <- as.matrix(cell_counts_rel_selected[,-1])
rownames(cell_counts_rel_mtrx) <- cell_counts_rel_selected$patient_id

###principal component analysis

cell_counts_pca <- as.data.frame(scores(pca(cell_counts_rel_mtrx, nPcs = 10)))
cell_counts_pca$patient_id <- rownames(cell_counts_pca)

histo_data <- unique(FetchData(epi_tumor, c("patient_id", "histo_subtype")))
cell_counts_pca <- full_join(cell_counts_pca, histo_data, by = "patient_id")

ggplot(cell_counts_pca, aes(x = PC1, y = PC2, color = histo_subtype, label = patient_id)) +
  geom_point(size = 2) +
  geom_text(aes(label = patient_id), color = "black", vjust = 1.5, size = 2) +
  scale_color_manual(values = use_colors) +
  scale_x_continuous(limits = c(-1, 1.2)) +
  scale_y_continuous(limits = c(-1.2, 0.8))
ggsave2("Fig5A.pdf", path = "../results", width = 15, height = 9, units = "cm")

###add mean tumor cell signature score

epi_pca <- AddModuleScore(epi_pca, features = list(c("SCGB3A1", "PIGR", "NAPSA", "C4BPA", "SCGB3A2", "HLA-DRA", "CD74", "ADGRF5", "C16orf89", "FOLR1", "SELENBP1", "HLA-DRB1", "ID4", "MGP", "AQP3", "CA2", "LHX9", "HLA-DPB1", "FMO5", "GKN2", "C5", "MUC1", "NPC2", "RNASE1", "PIK3C2G", "SFTA2", "SLC34A2", "HLA-DPA1", "FGFR3", "PGC")), name = "tumor_signature")
epi_pca <- AddModuleScore(epi_pca, features = list(c("DSG2", "CAMK2N1", "FAM3C", "KRT7", "IFI27", "SLC2A1", "MARCKS", "PLAU", "AHNAK2", "PERP", "S100A4", "KRT19", "COL6A1", "UACA", "COL17A1", "CDA", "TPM2", "S100A16", "KRT8", "PRSS23", "DST", "LAMC2", "S100P", "PRSS3", "LAMA3", "DSP", "ITGA3", "MDK", "FAM83A", "ITGB4")), name = "anti_tumor_signature")

epi_tumor_data <- FetchData(epi_pca, c(colnames(epi_pca@meta.data)))

epi_tumor_signature <- epi_tumor_data %>%
  group_by(patient_id) %>% 
  mutate(diff_mean = mean(tumor_signature1)) %>%
  mutate(undiff_mean = mean(anti_tumor_signature1)) %>%
  dplyr::select(patient_id, diff_mean, undiff_mean) %>%
  unique()

###heatmap of normalized immune and stromal cell counts, and mean tumor cell signatures

cell_counts_rel_heatmap <- full_join(cell_counts_rel_selected, cell_counts_pca, by = "patient_id")
cell_counts_rel_heatmap <- full_join(cell_counts_rel_heatmap, epi_tumor_signature, by = "patient_id")

cell_counts_rel_heatmap_ordered <- cell_counts_rel_heatmap[order(cell_counts_rel_heatmap$PC1),]

cell_counts_rel_heatmap1 <- cell_counts_rel_heatmap_ordered[,2:ncol(celltype_selected)]
cell_counts_rel_heatmap2 <- cell_counts_rel_heatmap_ordered[,c("diff_mean", "undiff_mean")]

rownames(cell_counts_rel_heatmap1) <- cell_counts_rel_heatmap_ordered$patient_id
rownames(cell_counts_rel_heatmap2) <- cell_counts_rel_heatmap_ordered$patient_id

cell_counts_rel_heatmap1 <- scale(cell_counts_rel_heatmap1)

pdf("../results/Fig5B_immune_and_stromal_cells.pdf", width = 7, height = 7)
heatmap.2(as.matrix(t(cell_counts_rel_heatmap1)), trace = "none", density.info = "none", scale = "none", symbreaks = F, col = hcl.colors(100, palette = "orrd", rev = T), Rowv = T, Colv = F, dendrogram = "row", margins = c(5,10))
dev.off()

pdf("../results/Fig5B_tumor_cell_signatures.pdf", width = 7, height = 5)
heatmap.2(as.matrix(t(cell_counts_rel_heatmap2)), trace = "none", density.info = "none", scale = "none", col = bluered(100), Rowv = F, Colv = F, dendrogram = "none")
dev.off()

###correlation plots

# matrix of the p-value of the correlation (from http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram)
cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], method = "spearman", ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

p.mat <- cor.mtest(cell_counts_rel_mtrx)

test <- cor(cell_counts_rel_mtrx, method = "spearman")

png("../results/SuppFig8A.png", width = 30, height = 30, units = "cm", res = 600)
corrplot(test, p.mat = p.mat, sig.level = 0.05, insig = "blank", type = "full", tl.col = "black")
dev.off()

#pdf("output/fig5/corrplot.pdf", width = 7, height = 7)
#corrplot(test, p.mat = p.mat, sig.level = 0.05, insig = "blank", type = "full", tl.col = "black")
#dev.off()

### correlation network plot

test_net <- test

for (i in 1:nrow(test_net)){
  for (j in 1:ncol(test_net)){
    if(p.mat[i,j] > 0.05){
      test_net[i,j] <- 0
    }
    if(test_net[i,j] == 1){
      test_net[i,j] <- 0
    }
  }
}

library(qgraph)

pdf("../results/Fig5C.pdf", width = 6, height = 6)
qgraph(test_net, layout = "spring", threshold = 0.7, labels = colnames(test_net), label.scale = F, label.cex = 0.7, theme = "colorblind", vsize = 3, color = "grey", borders = F)
dev.off()

#####scatter plots

cell_counts_rel <- mutate(cell_counts_rel, T_CD8_exhausted = T_CD8_1+T_CD8_2)

###plots stromal

ggplot(cell_counts_rel, aes(x = Myofibroblast1, y = Myofibroblast2)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.8)) +
  scale_y_continuous(limits = c(0,0.8)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$Myofibroblast1, cell_counts_rel$Myofibroblast2, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$Myofibroblast1, cell_counts_rel$Myofibroblast2, method = "spearman")$estimate))
ggsave2("Fig3E.pdf", path = "../results", width = 10, height = 9, units = "cm")

###plots immune

ggplot(cell_counts_rel, aes(x = CD14_Macrophages1, y = CD14_Macrophages2)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.6)) +
  scale_y_continuous(limits = c(0,0.8)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$CD14_Macrophages1, cell_counts_rel$CD14_Macrophages2, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$CD14_Macrophages1, cell_counts_rel$CD14_Macrophages2, method = "spearman")$estimate))
ggsave2("Fig4D_Macrophages1_2.pdf", path = "../results", width = 10, height = 9, units = "cm")

ggplot(cell_counts_rel, aes(x = CD14_Macrophages2, y = Plasmacytoid_Dendritic)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.8)) +
  scale_y_continuous(limits = c(0,0.06)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$CD14_Macrophages2, cell_counts_rel$Plasmacytoid_Dendritic, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$CD14_Macrophages2, cell_counts_rel$Plasmacytoid_Dendritic, method = "spearman")$estimate))
ggsave2("Fig4D_Macrophages2_pDC.pdf", path = "../results", width = 10, height = 9, units = "cm")

ggplot(cell_counts_rel, aes(x = T_reg, y = T_CD8_exhausted)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.5)) +
  scale_y_continuous(limits = c(0,0.85)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$T_reg, cell_counts_rel$T_CD8_exhausted, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$T_reg, cell_counts_rel$T_CD8_exhausted, method = "spearman")$estimate))
ggsave2("Fig4G_T_reg_exhausted.pdf", path = "../results", width = 10, height = 9, units = "cm")

ggplot(cell_counts_rel, aes(x = CD14_Macrophages2, y = T_CD8_exhausted)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.8)) +
  scale_y_continuous(limits = c(0,0.85)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$CD14_Macrophages2, cell_counts_rel$T_CD8_exhausted, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$CD14_Macrophages2, cell_counts_rel$T_CD8_exhausted, method = "spearman")$estimate))
ggsave2("Fig4G_Macrophages_exhausted.pdf", path = "../results", width = 10, height = 9, units = "cm")

ggplot(cell_counts_rel, aes(x = Plasmacytoid_Dendritic, y = T_CD8_exhausted)) +
  geom_point(aes(color = patient_id), size = 4) +
  scale_x_continuous(limits = c(0,0.06)) +
  scale_y_continuous(limits = c(0,0.85)) +
  scale_color_manual(values = use_colors) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(cell_counts_rel$Plasmacytoid_Dendritic, cell_counts_rel$T_CD8_exhausted, method = "spearman")$p.value, 
                 "\nr = ",cor.test(cell_counts_rel$Plasmacytoid_Dendritic, cell_counts_rel$T_CD8_exhausted, method = "spearman")$estimate))
ggsave2("Fig4G_pDC_exhausted.pdf", path = "../results", width = 10, height = 9, units = "cm")

###compare cell counts with IHC quantification

CD123 <- read_excel("../data/IHC_quantification/Quantification_CD123.xlsx")
CD123_scrnaseq <- cell_counts_rel %>% select(Plasmacytoid_Dendritic)
CD123 <- full_join(CD123, CD123_scrnaseq, by = "patient_id")
CD123 <- pivot_longer(CD123, cols = starts_with("image"), names_to = "image_nr", values_to = "count")
CD123_average <- CD123 %>% group_by(patient_id) %>% mutate(mean = mean(count)) %>% mutate(sd = sd(count)) %>% select(c(patient_id, Plasmacytoid_Dendritic, mean, sd)) %>% unique()
ggplot() +
  geom_jitter(data = CD123, mapping = aes(x = Plasmacytoid_Dendritic, y = count, color = patient_id), size = 0.5, alpha = 0.3) +
  geom_smooth(data = CD123_average, mapping = aes(x = Plasmacytoid_Dendritic, y = mean, color = patient_id), method = "lm", se = FALSE, size = 1, color = "black") +
  geom_errorbar(data = CD123_average, aes(x = Plasmacytoid_Dendritic, ymin=mean-sd, ymax=mean+sd, color = patient_id), size = 0.5, width = 0.002) +
  geom_point(data = CD123_average, aes(x = Plasmacytoid_Dendritic, y=mean, color = patient_id), shape=15, size = 3) +
  scale_color_manual(values = use_colors) +
  ggtitle(paste0("p = ",cor.test(CD123_average$mean, CD123_average$Plasmacytoid_Dendritic, method = "pearson")$p.value, 
                 "\nr = ",cor.test(CD123_average$mean, CD123_average$Plasmacytoid_Dendritic, method = "pearson")$estimate))
ggsave2("Fig4E_pDC_CD123.pdf", path = "../results", width = 10, height = 9, units = "cm")

CTHRC1 <- read_excel("../data/IHC_quantification/Quantification_CTHRC1.xlsx")
CTHRC1_scrnaseq <- cell_counts_rel %>% select(Myofibroblast2)
CTHRC1 <- full_join(CTHRC1, CTHRC1_scrnaseq, by = "patient_id")
CTHRC1 <- pivot_longer(CTHRC1, cols = starts_with("image"), names_to = "image_nr", values_to = "count")
CTHRC1_average <- CTHRC1 %>% group_by(patient_id) %>% mutate(mean = mean(count)) %>% mutate(sd = sd(count)) %>% select(c(patient_id, Myofibroblast2, mean, sd)) %>% unique()
ggplot() +
  geom_jitter(data = CTHRC1, mapping = aes(x = Myofibroblast2, y = count, color = patient_id), size = 0.5, alpha = 0.3) +
  geom_smooth(data = CTHRC1_average, mapping = aes(x = Myofibroblast2, y = mean, color = patient_id), method = "lm", se = FALSE, size = 1, color = "black") +
  geom_errorbar(data = CTHRC1_average, aes(x = Myofibroblast2, ymin=mean-sd, ymax=mean+sd, color = patient_id), size = 0.5) +
  geom_point(data = CTHRC1_average, aes(x = Myofibroblast2, y=mean, color = patient_id), shape=15, size = 3) +
  scale_color_manual(values = use_colors) +
  ggtitle(paste0("p = ",cor.test(CTHRC1_average$mean, CTHRC1_average$Myofibroblast2, method = "pearson")$p.value, 
                 "\nr = ",cor.test(CTHRC1_average$mean, CTHRC1_average$Myofibroblast2, method = "pearson")$estimate))
ggsave2("Fig3F_Myofibroblast2_CTHRC1.pdf", path = "../results", width = 10, height = 9, units = "cm")

CXCL9 <- read_excel("../data/IHC_quantification/Quantification_CXCL9.xlsx")
CXCL9_scrnaseq <- cell_counts_rel %>% select(CD14_Macrophages2)
CXCL9 <- full_join(CXCL9, CXCL9_scrnaseq, by = "patient_id")
CXCL9 <- pivot_longer(CXCL9, cols = starts_with("image"), names_to = "image_nr", values_to = "count")
CXCL9_average <- CXCL9 %>% group_by(patient_id) %>% mutate(mean = mean(count)) %>% mutate(sd = sd(count)) %>% select(c(patient_id, CD14_Macrophages2, mean, sd)) %>% unique()
ggplot() +
  geom_jitter(data = CXCL9, mapping = aes(x = CD14_Macrophages2, y = count, color = patient_id), size = 0.5, alpha = 0.3, width = 0.02) +
  geom_smooth(data = CXCL9_average, mapping = aes(x = CD14_Macrophages2, y = mean, color = patient_id), method = "lm", se = FALSE, size = 1, color = "black") +
  geom_errorbar(data = CXCL9_average, aes(x = CD14_Macrophages2, ymin=mean-sd, ymax=mean+sd, color = patient_id), size = 0.5, width = 0.03) +
  geom_point(data = CXCL9_average, aes(x = CD14_Macrophages2, y=mean, color = patient_id), shape=15, size = 3) +
  scale_color_manual(values = use_colors) +
  ggtitle(paste0("p = ",cor.test(CXCL9_average$mean, CXCL9_average$CD14_Macrophages2, method = "pearson")$p.value, 
                 "\nr = ",cor.test(CXCL9_average$mean, CXCL9_average$CD14_Macrophages2, method = "pearson")$estimate))
ggsave2("Fig4E_CD14_Macrophages_2_CXCL9.pdf", path = "../results", width = 10, height = 9, units = "cm")

########################################TCGA dataset

library(TCGAbiolinks)
library(GSVA)
library(SummarizedExperiment)
library(Seurat)
library(biomaRt)
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
library(corrplot)
library(stringr)
library(survival)
library(survminer)
library(plotrix)
library(readr)
library(HDF5Array)

theme_set(theme_cowplot())

###load gene expression data

luad_se <- loadHDF5SummarizedExperiment(dir="../data/TCGA LUAD")
luad_mut <- read_csv("../data/data/TCGA_LUAD_mutations.csv")

###replace gene symbols, calculate log2-transformed matrix

luad_data <- as.data.frame(assay(luad_se))

luad_data$hgnc_symbol <- luad_se@rowRanges$external_gene_name

luad_data <- luad_data[!(duplicated(luad_data$hgnc_symbol)),]
rownames(luad_data) <- luad_data$hgnc_symbol

luad_data <- luad_data[,1:(ncol(luad_data)-1)]
luad_data <- as.matrix(luad_data)

log_luad_data <- log2(luad_data+1)

###calculate cell cluster marker genes

#imm_anno <- readRDS("seurat_objects/imm_anno.RDS")

imm_anno@meta.data$cell_type_imm <- ordered(imm_anno@meta.data$cell_type_imm, levels = c("Alveolar_Macrophages1",
                                                                                         "Alveolar_Macrophages2",
                                                                                         "Alveolar_Macrophages3",
                                                                                         "CD14_Macrophages1",
                                                                                         "CD14_Macrophages2",
                                                                                         "CD14_Macrophages3",
                                                                                         "CD14_Macrophages4",
                                                                                         "CD14_Macrophages5",
                                                                                         "Macrophages_Proliferating",
                                                                                         "Monocytes",
                                                                                         "Myeloid_Dendritic",
                                                                                         "Plasmacytoid_Dendritic",
                                                                                         "Mast",
                                                                                         "T_conv1",
                                                                                         "T_conv2",
                                                                                         "T_reg",
                                                                                         "T_CD8_1",
                                                                                         "T_CD8_2",
                                                                                         "T_CD8_3",
                                                                                         "T_CD8_Proliferating",
                                                                                         "NK_cells",
                                                                                         "B_cells",
                                                                                         "Plasma"))

imm_lympho <- subset(imm_anno, subset = cell_type_imm %in% c("T_conv1",
                                                             "T_conv2",
                                                             "T_reg",
                                                             "T_CD8_1",
                                                             "T_CD8_2",
                                                             "T_CD8_3",
                                                             "T_CD8_Proliferating",
                                                             "NK_cells",
                                                             "B_cells",
                                                             "Plasma"))
imm_lympho <- ScaleData(imm_lympho)

imm_myelo <- subset(imm_anno, subset = cell_type_imm %in% c("Alveolar_Macrophages1",
                                                            "Alveolar_Macrophages2",
                                                            "Alveolar_Macrophages3",
                                                            "CD14_Macrophages1",
                                                            "CD14_Macrophages2",
                                                            "CD14_Macrophages3",
                                                            "CD14_Macrophages4",
                                                            "CD14_Macrophages5",
                                                            "Macrophages_Proliferating",
                                                            "Monocytes",
                                                            "Myeloid_Dendritic",
                                                            "Plasmacytoid_Dendritic",
                                                            "Mast"))
imm_myelo <- ScaleData(imm_myelo)

#str_anno <- readRDS("seurat_objects/str_anno.RDS")

str_anno@meta.data$cell_type_str <- factor(str_anno@meta.data$cell_type_str, levels = c("Endothelial1",
                                                                                        "Endothelial2",
                                                                                        "Endothelial3",
                                                                                        "Endothelial4",
                                                                                        "Endothelial5",
                                                                                        "Endothelial6",
                                                                                        "Endothelial7",
                                                                                        "Lymphaticendothelial",
                                                                                        "Fibroblast1",
                                                                                        "Fibroblast2",
                                                                                        "Myofibroblast1",
                                                                                        "Myofibroblast2",
                                                                                        "Smoothmuscle1",
                                                                                        "Smoothmuscle2",
                                                                                        "Mesothelial"))

str_endo <- subset(str_anno, subset = cell_type_str %in% c("Endothelial1",
                                                           "Endothelial2",
                                                           "Endothelial3",
                                                           "Endothelial4",
                                                           "Endothelial5",
                                                           "Endothelial6",
                                                           "Endothelial7",
                                                           "Lymphaticendothelial"))
str_endo <- ScaleData(str_endo)

str_fibro <- subset(str_anno, subset = cell_type_str %in% c("Fibroblast1",
                                                            "Fibroblast2",
                                                            "Myofibroblast1",
                                                            "Myofibroblast2",
                                                            "Smoothmuscle1",
                                                            "Smoothmuscle2",
                                                            "Mesothelial"))
str_fibro <- ScaleData(str_fibro)

Idents(imm_myelo) <- imm_myelo@meta.data$cell_type_imm
Idents(imm_lympho) <- imm_lympho@meta.data$cell_type_imm
Idents(str_endo) <- str_endo@meta.data$cell_type_str
Idents(str_fibro) <- str_fibro@meta.data$cell_type_str

myelo_markers <- FindAllMarkers(imm_myelo, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)
write_csv(myelo_markers, path = "../results/SuppTable9.csv")

lympho_markers <- FindAllMarkers(imm_lympho, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)
write_csv(lympho_markers, path = "../results/SuppTable10.csv")

endo_markers <- FindAllMarkers(str_endo, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)
write_csv(endo_markers, path = "../results/SuppTable7.csv")

fibro_markers <- FindAllMarkers(str_fibro, only.pos = T, min.pct = 0.25, min.diff.pct = 0.25)
write_csv(fibro_markers, path = "../results/SuppTable8.csv")

clusters_myelo <- as.character(unique(myelo_markers$cluster))
clusters_lympho <- as.character(unique(lympho_markers$cluster))
clusters_endo <- as.character(unique(endo_markers$cluster))
clusters_fibro <- as.character(unique(fibro_markers$cluster))

markers <- list()

for (i in seq_along(clusters_myelo)){
  markers[paste0(clusters_myelo[i])] <- myelo_markers %>% filter(cluster == clusters_myelo[i]) %>% dplyr::select(gene)
}

for (i in seq_along(clusters_lympho)){
  markers[paste0(clusters_lympho[i])] <- lympho_markers %>% filter(cluster == clusters_lympho[i]) %>% dplyr::select(gene)
}

for (i in seq_along(clusters_endo)){
  markers[paste0(clusters_endo[i])] <- endo_markers %>% filter(cluster == clusters_endo[i]) %>% dplyr::select(gene)
}

for (i in seq_along(clusters_fibro)){
  markers[paste0(clusters_fibro[i])] <- fibro_markers %>% filter(cluster == clusters_fibro[i]) %>% dplyr::select(gene)
}

### genes from PCA of tumor epithelial cells

markers[["PC1_pos"]] <- c("SCGB3A1", "PIGR", "NAPSA", "C4BPA", "SCGB3A2", "HLA-DRA", "CD74", "ADGRF5", "C16orf89", "FOLR1", "SELENBP1", "HLA-DRB1", "ID4", "MGP", "AQP3", "CA2", "LHX9", "HLA-DPB1", "FMO5", "GKN2", "C5", "MUC1", "NPC2", "RNASE1", "PIK3C2G", "SFTA2", "SLC34A2", "HLA-DPA1", "FGFR3", "PGC")
markers[["PC1_neg"]] <- c("DSG2", "CAMK2N1", "FAM3C", "KRT7", "IFI27", "SLC2A1", "MARCKS", "PLAU", "AHNAK2", "PERP", "S100A4", "KRT19", "COL6A1", "UACA", "COL17A1", "CDA", "TPM2", "S100A16", "KRT8", "PRSS23", "DST", "LAMC2", "S100P", "PRSS3", "LAMA3", "DSP", "ITGA3", "MDK", "FAM83A", "ITGB4")

###short signatures

for (i in seq_along(clusters_myelo)){
  markers[paste0(clusters_myelo[i], "_short")] <- myelo_markers %>% filter(cluster == clusters_myelo[i]) %>% top_n(5, wt = avg_logFC) %>% dplyr::select(gene)
}

for (i in seq_along(clusters_lympho)){
  markers[paste0(clusters_lympho[i], "_short")] <- lympho_markers %>% filter(cluster == clusters_lympho[i]) %>% top_n(5, wt = avg_logFC)%>% dplyr::select(gene)
}

for (i in seq_along(clusters_endo)){
  markers[paste0(clusters_endo[i], "_short")] <- endo_markers %>% filter(cluster == clusters_endo[i]) %>% top_n(5, wt = avg_logFC)%>% dplyr::select(gene)
}

for (i in seq_along(clusters_fibro)){
  markers[paste0(clusters_fibro[i], "_short")] <- fibro_markers %>% filter(cluster == clusters_fibro[i]) %>% top_n(5, wt = avg_logFC)%>% dplyr::select(gene)
}

markers[["PC1_pos_short"]] <- c("SCGB3A1", "PIGR", "NAPSA", "C4BPA", "SCGB3A2")
markers[["PC1_neg_short"]] <- c("DSG2", "CAMK2N1", "FAM3C", "KRT7", "IFI27")

###complete pattern signatures

#long
markers[["N3MC_complete"]] <- unique(c(markers[["CD14_Macrophages1"]], markers[["Myofibroblast1"]], markers[["NK_cells"]], markers[["Myeloid_Dendritic"]], markers[["T_conv1"]], markers[["PC1_pos"]]))
markers[["CP2E_complete"]] <- unique(c(markers[["CD14_Macrophages2"]], markers[["Myofibroblast2"]], markers[["Plasmacytoid_Dendritic"]], markers[["T_CD8_1"]], markers[["T_CD8_2"]], markers[["PC1_neg"]]))

#short
markers[["N3MC_complete_short"]] <- unique(c(markers[["CD14_Macrophages1_short"]], markers[["Myofibroblast1_short"]], markers[["NK_cells_short"]], markers[["Myeloid_Dendritic_short"]], markers[["T_conv1_short"]], markers[["PC1_pos_short"]]))
markers[["CP2E_complete_short"]] <- unique(c(markers[["CD14_Macrophages2_short"]], markers[["Myofibroblast2_short"]], markers[["Plasmacytoid_Dendritic_short"]], markers[["T_CD8_1_short"]], markers[["T_CD8_2_short"]], markers[["PC1_neg_short"]]))

###simple pattern signatures

#long
markers[["N3MC_simple"]] <- unique(c(markers[["CD14_Macrophages1"]], markers[["Myofibroblast1"]]))
markers[["CP2E_simple"]] <- unique(c(markers[["Myofibroblast2"]], markers[["PC1_neg"]]))

#short
markers[["N3MC_simple_short"]] <- unique(c(markers[["CD14_Macrophages1_short"]], markers[["Myofibroblast1_short"]]))
markers[["CP2E_simple_short"]] <- unique(c(markers[["Myofibroblast2_short"]], markers[["PC1_neg_short"]]))

###save gene sets
write_csv(as_data_frame(markers["N3MC_complete"]), path = "../results/SuppTable11_1.csv")
write_csv(as_data_frame(markers["CP2E_complete"]), path = "../results/SuppTable11_2.csv")
write_csv(as_data_frame(markers["N3MC_simple_short"]), path = "../results/SuppTable11_3.csv")
write_csv(as_data_frame(markers["CP2E_simple_short"]), path = "../results/SuppTable11_4.csv")

###single-sample Gene set Enrichment Analysis (ssGSEA)

luad_gsva <- gsva(expr = log_luad_data, gset.idx.list = markers, method = "ssgsea", verbose = TRUE, kcdf = "Gaussian")

luad_gsva <- as.data.frame(t(luad_gsva))

###calculate ratios

luad_gsva <- mutate(luad_gsva, PC1_ratio = PC1_pos / PC1_neg)
luad_gsva <- mutate(luad_gsva, Myofibro_ratio = Myofibroblast1 / Myofibroblast2)
luad_gsva <- mutate(luad_gsva, Macro_ratio = CD14_Macrophages1 / CD14_Macrophages2)
luad_gsva <- mutate(luad_gsva, Pattern_complete_ratio = N3MC_complete / CP2E_complete)
luad_gsva <- mutate(luad_gsva, Pattern_simple_ratio = N3MC_simple / CP2E_simple)

luad_gsva <- mutate(luad_gsva, PC1_ratio_short = PC1_pos_short / PC1_neg_short)
luad_gsva <- mutate(luad_gsva, Myofibro_ratio_short = Myofibroblast1_short / Myofibroblast2_short)
luad_gsva <- mutate(luad_gsva, Macro_ratio_short = CD14_Macrophages1_short / CD14_Macrophages2_short)
luad_gsva <- mutate(luad_gsva, Pattern_complete_short_ratio = N3MC_complete_short / CP2E_complete_short)
luad_gsva <- mutate(luad_gsva, Pattern_simple_short_ratio = N3MC_simple_short / CP2E_simple_short)

###add follow-up data

luad_gsva$days_to_death <- luad_se$days_to_death
luad_gsva$vital_status <- luad_se$vital_status
luad_gsva$days_to_last_follow_up <- luad_se$days_to_last_follow_up

luad_gsva$overall_survival <- ifelse(luad_gsva$vital_status == "Dead", luad_gsva$days_to_death, luad_gsva$days_to_last_follow_up)

###binarize data

luad_gsva <- mutate(luad_gsva, vital_status_bin = ifelse(vital_status == "Alive", 0, 1))

ES_names <- colnames(luad_gsva[1:98])

for (i in seq_along(ES_names)) {
  luad_gsva[,paste0(ES_names[i], "_bin")] <- ifelse(luad_gsva[,ES_names[i]] > median(luad_gsva[,ES_names[i]]), 1, 0)
}

###correlation plot all clusters

# matrix of the p-value of the correlation (from http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram)
cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], method = "spearman", ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

p.mat <- cor.mtest(as.matrix(luad_gsva[,1:40]))

corr_tcga <- cor(as.matrix(luad_gsva[,1:40]), method = "spearman")

png(file = "../results/SuppFig8B.png", width = 30, height = 30, units = "cm", res = 600)
corrplot(corr_tcga, p.mat = p.mat, sig.level = 0.05, insig = "blank", type = "full", tl.col = "black")
dev.off()

#pdf(file = "output/fig5/corrplot_tcga.pdf", width = 7, height = 7)
#corrplot(corr_tcga, p.mat = p.mat, sig.level = 0.05, insig = "blank", type = "full", tl.col = "black")
#dev.off()

###correlation plot for selected cell clusters

### PC1 ~ Macrophages

ggplot(luad_gsva, aes(y = PC1_pos, x = CD14_Macrophages1)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_pos, luad_gsva$CD14_Macrophages1, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_pos, luad_gsva$CD14_Macrophages1, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_pos_CD14_Macrophages1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_pos, x = CD14_Macrophages2)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_pos, luad_gsva$CD14_Macrophages2, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_pos, luad_gsva$CD14_Macrophages2, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_pos_CD14_Macrophages2.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_neg, x = CD14_Macrophages1)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_neg, luad_gsva$CD14_Macrophages1, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_neg, luad_gsva$CD14_Macrophages1, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_neg_CD14_Macrophages1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_neg, x = CD14_Macrophages2)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_neg, luad_gsva$CD14_Macrophages2, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_neg, luad_gsva$CD14_Macrophages2, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_neg_CD14_Macrophages2.pdf", path = "../results", width = 10, height = 10, units = "cm")

### PC1 ~ Myofibroblasts

ggplot(luad_gsva, aes(y = PC1_pos, x = Myofibroblast1)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_pos, luad_gsva$Myofibroblast1, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_pos, luad_gsva$Myofibroblast1, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_pos_Myofibroblast1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_pos, x = Myofibroblast2)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_pos, luad_gsva$Myofibroblast2, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_pos, luad_gsva$Myofibroblast2, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_pos_Myofibroblast2.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_neg, x = Myofibroblast1)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_neg, luad_gsva$Myofibroblast1, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_neg, luad_gsva$Myofibroblast1, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_neg_Myofibroblast1.pdf", path = "../results", width = 10, height = 10, units = "cm")

ggplot(luad_gsva, aes(y = PC1_neg, x = Myofibroblast2)) +
  geom_point(color = "grey", size = 0.1) +
  geom_smooth(method = "lm", se = FALSE, size = 1, color = "black") +
  ggtitle(paste0("p = ",cor.test(luad_gsva$PC1_neg, luad_gsva$Myofibroblast2, method = "spearman")$p.value, 
                 "\nrho = ",cor.test(luad_gsva$PC1_neg, luad_gsva$Myofibroblast2, method = "spearman")$estimate))
ggsave2("Fig5E_PC1_neg_Myofibroblast2.pdf", path = "../results", width = 10, height = 10, units = "cm")

###Survival analysis

# Kaplan-Meyer and Log-rank

for (i in c(104:201)) {
  ggsurvplot(survfit(Surv(as.numeric(overall_survival/30.4375), vital_status_bin)~luad_gsva[[i]], data = luad_gsva), pval = T, pval.method = T, title = paste0(colnames(luad_gsva)[i]), xlim = c(0,120), break.x.by = 30, palette = c("Red", "Blue"))
  ggsave2(paste0("Fig5F_SuppFig11_", colnames(luad_gsva)[i],".pdf"), path = "../results", width = 8, height = 10, units = "cm")
}

# Cox regression

summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Macro_ratio_bin, data = luad_gsva)) 
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Myofibro_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ PC1_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Pattern_complete_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Pattern_simple_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Pattern_complete_short_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Pattern_simple_short_ratio_bin, data = luad_gsva))
summary(coxph(Surv(as.numeric(overall_survival), vital_status_bin) ~ Macro_ratio_bin + Myofibro_ratio_bin + PC1_ratio_bin, data = luad_gsva))

###add mutation data

#genes from nNGMv2 Panel
mutations <- c("ALK", "BRAF", "CTNNB1", "EGFR", "ERBB2", "FGFR1", "FGFR2", "FGFR3", "FGFR4", "HRAS", "IDH1", "KEAP1", "KRAS", "MAP2K1", "MET", "NRAS", "NTRK1", "NTRK2", "NTRK3", "PIK3CA", "PTEN", "RET", "ROS1", "STK11", "TP53")

patients_mut <- list()

for (i in seq_along(mutations)){
  patients_mut[[i]] <- filter(luad_mut, Hugo_Symbol == mutations[i] & Variant_Classification != "Silent")$Tumor_Sample_Barcode %>% as.vector() %>% substr(1,12)
}

names(patients_mut) <- mutations

luad_gsva_mut <- luad_gsva

luad_gsva_mut[, mutations] <- 0
luad_gsva_mut$samples <- rownames(luad_gsva_mut) %>% substr(1,12)

luad_gsva_mut$Pattern_ratio_complete_bin <- ifelse(luad_gsva_mut$Pattern_complete_ratio_bin == 1, "N3MC_complete", "CP2E_complete")
luad_gsva_mut$Pattern_ratio_simple_short_bin <- ifelse(luad_gsva_mut$Pattern_simple_short_ratio_bin == 1, "N3MC_simple_short", "CP2E_simple_short")

luad_gsva_mut$TMB <- as.character(luad_se$paper_Nonsilent.Mutations.per.Mb) %>% gsub(pattern = ",", replacement = ".", fixed = T) %>% as.numeric()

for(i in seq_along(mutations)){
  luad_gsva_mut[, mutations[i]] <- ifelse(luad_gsva_mut$sample %in% patients_mut[[i]], "mut", "wt")
}

luad_gsva_mut <- filter(luad_gsva_mut, samples %in% substr(luad_mut$Tumor_Sample_Barcode, 1, 12))

#Pattern_ratio_bin

for(i in seq_along(mutations)){
  mut_counts <- luad_gsva_mut %>% group_by(Pattern_ratio_complete_bin) %>% dplyr::count(eval(parse(text = mutations[i]))) %>% pivot_wider(id_cols = Pattern_ratio_complete_bin, names_from = "eval(parse(text = mutations[i]))", values_from = n)

  mut_counts[is.na(mut_counts)] <- 0

  chi <- chisq.test(mut_counts[,2:3])

  ggplot(luad_gsva_mut) + 
    geom_bar(aes(x = Pattern_ratio_complete_bin, fill = eval(parse(text = mutations[i]))), position = "fill") +
    theme(legend.title = element_blank()) +
    ggtitle(paste0(mutations[i], " X虏 = ", chi$statistic, " p = ", chi$p.value))
  ggsave2(paste0("Fig5G_mut_", mutations[[i]], ".pdf"), path = "../results", width = 20, height = 10, units = "cm")
}

#for(i in seq_along(mutations)){
#  mut_counts <- luad_gsva_mut %>% group_by(Pattern_ratio_simple_short_bin) %>% dplyr::count(eval(parse(text = mutations[i]))) %>% pivot_wider(id_cols = Pattern_ratio_simple_short_bin, names_from = "eval(parse(text = mutations[i]))", values_from = n)
#  
#  mut_counts[is.na(mut_counts)] <- 0
#  
#  chi <- chisq.test(mut_counts[,2:3])
#  
#  ggplot(luad_gsva_mut) + 
#    geom_bar(aes(x = Pattern_ratio_simple_short_bin, fill = eval(parse(text = mutations[i]))), position = "fill") +
#    theme(legend.title = element_blank()) +
#    ggtitle(paste0(mutations[i], " X虏 = ", chi$statistic, " p = ", chi$p.value))
#  ggsave2(paste0("Fig5G_mut_simple_short_", mutations[[i]], ".pdf"), path = "../results", width = 20, height = 10, units = "cm")
#}

###Tumor mutational burden

ggplot(luad_gsva_mut) +
  geom_boxplot(aes(x = Pattern_ratio_complete_bin, y = TMB)) +
  ggtitle(paste0("p = ", t.test(formula = TMB~Pattern_ratio_complete_bin, data = luad_gsva_mut, alternative = "two.sided", paired = F)$p.value))
ggsave2("Fig5G_TMB.pdf", path = "../results", width = 20, height = 10, units = "cm")

#ggplot(luad_gsva_mut) +
#  geom_boxplot(aes(x = Pattern_ratio_simple_short_bin, y = TMB)) +
#  ggtitle(paste0("p = ", t.test(formula = TMB~Pattern_ratio_simple_short_bin, data = luad_gsva_mut, alternative = "two.sided", paired = F)$p.value))
#ggsave2("Fig5G_TMB_simple_short.pdf", path = "../results", width = 20, height = 10, units = "cm")

###add fusion data

fusions.data <- data.frame(luad_se$paper_Fusions, luad_se$barcode)
fusions.data$barcode <- fusions.data$luad_se.barcode

luad_gsva_mut$barcode <- rownames(luad_gsva_mut)

luad_gsva_fus <- inner_join(fusions.data, luad_gsva_mut, by = "barcode")

luad_gsva_fus$alk_fusions <- ifelse(grepl("ALK", luad_gsva_fus$luad_se.paper_Fusions) == T, "fus", "wt")
luad_gsva_fus$ret_fusions <- ifelse(grepl("RET", luad_gsva_fus$luad_se.paper_Fusions) == T, "fus", "wt")
luad_gsva_fus$ros_fusions <- ifelse(grepl("ROS", luad_gsva_fus$luad_se.paper_Fusions) == T, "fus", "wt")

fusions <- c("alk_fusions", "ret_fusions", "ros_fusions")

#Pattern_ratio_bin

for(i in seq_along(fusions)){
  fus_counts <- luad_gsva_fus %>% group_by(Pattern_ratio_complete_bin) %>% dplyr::count(eval(parse(text = fusions[i]))) %>% pivot_wider(id_cols = Pattern_ratio_complete_bin, names_from = "eval(parse(text = fusions[i]))", values_from = n)

  fus_counts[is.na(fus_counts)] <- 0

  chi <- chisq.test(fus_counts[,2:3])

  ggplot(luad_gsva_fus) + 
    geom_bar(aes(x = Pattern_ratio_complete_bin, fill = eval(parse(text = fusions[i]))), position = "fill") +
    theme(legend.title = element_blank()) +
    ggtitle(paste0(fusions[i], " X? = ", chi$statistic, " p = ", chi$p.value))
  ggsave2(paste0("Fig5G_fus_", fusions[[i]], ".pdf"), path = "../results", width = 20, height = 10, units = "cm")      
}

#for(i in seq_along(fusions)){
#  fus_counts <- luad_gsva_fus %>% group_by(Pattern_ratio_simple_short_bin) %>% dplyr::count(eval(parse(text = fusions[i]))) %>% pivot_wider(id_cols = Pattern_ratio_simple_short_bin, names_from = "eval(parse(text = fusions[i]))", values_from = n)
#  
#  fus_counts[is.na(fus_counts)] <- 0
#  
#  chi <- chisq.test(fus_counts[,2:3])
#  
#  ggplot(luad_gsva_fus) + 
#    geom_bar(aes(x = Pattern_ratio_simple_short_bin, fill = eval(parse(text = fusions[i]))), position = "fill") +
#    theme(legend.title = element_blank()) +
#    ggtitle(paste0(fusions[i], " X? = ", chi$statistic, " p = ", chi$p.value))
#  ggsave2(paste0("Fig5_fus_simple_short_", fusions[[i]], ".pdf"), path = "../results", width = 20, height = 10, units = "cm")      
#}

print(nrow(luad_gsva))
print(nrow(luad_gsva_mut))
luad_gsva_mut %>% filter(TMB != "NA") %>% nrow() %>% print()

Reference

https://codeocean.com/capsule/8321305/tree/v1
P. Bischoff et al., “Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma,” Oncogene 2021 4050, vol. 40, no. 50, pp. 6748–6758, Oct. 2021, doi: 10.1038/s41388-021-02054-3.

Introduction

Somethimes we need to transform gene ID of MicroArray data to gene symbol. For example, "AFFX-LysX-3_at", we don’t known what it means until we transform it.

Code exmaple

1. Check the platform

   !Series_platform_id     "GPL15207"

   This text was presented in my data, so I known my data were produced from GPL15207.

2. install GEOquery package (Code is not fixed )

   if (!requireNamespace("BiocManager", quietly = TRUE))
   install.packages("BiocManager")
   BiocManager::install("GEOquery")

3. Transform
   I want to trasform rownames(mt) to gene symbol

   #get annotation file
   GPL15207 <- getGEO("GPL15207", destdir = ".")
   tran_dt <- Table(GPL15207)

   We need two column, ID and Gene symbol.
   ID means raw ID, such like AFFX-LysX-3_at
   Gene symbol was gene symbol that we need

   tran_dt <- tran_dt[,c("ID","Gene Symbol")]

   Tramsform (by your way)

   
   mt <- mt[rownames(mt) %in% tran_dt$ID,]

#### References
The manual of getGEO function
#### Donation



    

        
        
        
        
        
        
        
        
    

地贫公益片

Introduction

  Technological and computational advances in genomics and interactomics have made it possible to identify how disease mutations perturb protein–protein interaction (PPI) networks within human cells.
  A PPI in which there is significantenrichment in interface mutations in one or the other of the two protein-binding partners across individuals will be defined as an oncoPPI. For each gene gi and its PPI interfaces, we assume that the observed number of mutations for a given interface follows a binomial distribution, binomial (T,P_gi ), in which T is the total number of mutations observed in one gene and P_gi is the estimated mutation rate for the region of interest under the null hypothesis that the region was not recurrently mutated. Using length (g_i) to represent the length of the protein product of gene _gi, for each interface, we computed the P value—the probability of observing >k mutations around this interface out of T total mutations observed in this gene—using the following equation
  in which _gi = length of interface/ length(_gi).
Finally, they set the minimal P value across all the interfaces in a specific protein as the representative P value of its coding gene gi, denoted P(P_gi). The significance of each PPI is defined as the product of P values of the two proteins (gene products). All P values were adjusted for multiple testing using the Bonferroni correction.
  Authors found that somatic missense mutations are also significantly enriched in PPI interfaces compared to noninterfaces in 10,861 tumor exomes.

Code explanation

Read requred file

# read protein length
uniprot_length <- read.table("UniProt_len.txt",colClasses=c("character","numeric"))

psites = list()
data <- read.table('PPI_sites.txt',sep="\t",colClasses="character")  
for(i in 1:nrow(data)){
    psites[[data[i,1]]] <- unlist(strsplit(data[i,2],","))  # psites information
}

cancer_types <- sapply(list.files('HQinterfaces_mutations/'),function(x) unlist(strsplit(x,split="_"))[1])

for each protein, run binomial test

for(cancer in cancer_types){
    Mutations <- read.table(paste0('/data/Data/TCGA_Somatic_Mutations/Mutect2/TCGA.', cancer, '.mutect.NonSilent.maf'),
        skip=1,header=TRUE,colClasses="character",sep="\t") # total mutation
    PPI_mutations <- read.table(paste0('HQinterfaces_mutations/', cancer, '_HQinterfaces_mutations.txt'),colClasses="character")  # PPI mutation

    all_uniprot <- unique(PPI_mutations[,3])
    uniprot_gene_mapping <- rep(NA,length(all_uniprot))
    names(uniprot_gene_mapping) <- all_uniprot
    uniprot_pvalues <- rep(1,length(all_uniprot))
    names(uniprot_pvalues) <- all_uniprot
    for(uniprot in all_uniprot){
        uniprot_gene_mapping[uniprot] <- PPI_mutations[match(uniprot,PPI_mutations[,3]),2]
        mutation_number <- sum(Mutations[,68]==uniprot)
        PPI_mutation_number <- sum(PPI_mutations[,3]==uniprot)
        gene_len <- uniprot_length[match(uniprot,uniprot_length[,1]),2]
        bind_len <- length(psites[[uniprot]])
        if(PPI_mutation_number > 0 & mutation_number > 0){
            p <- binom.test(x=PPI_mutation_number,n=mutation_number,p=bind_len/gene_len,alternative="greater")$p.value
            uniprot_pvalues[uniprot] <- p
        }
    }
    write.table(cbind(uniprot_gene_mapping,uniprot_pvalues),
                file=paste("./uniprot_pvalues/",cancer,"_uniprot_pvalues.txt",sep=""),
                col.names=FALSE,
                sep="\t",
                quote=FALSE)
}

cancer <- cancer_types[1]
Mutations <- read.table(paste0('/data/Data/TCGA_Somatic_Mutations/Mutect2/TCGA.', cancer, '.mutect.NonSilent.maf'),
        skip=1,header=TRUE,colClasses="character",sep="\t")
PPI_mutations <- read.table(paste0('HQinterfaces_mutations/', cancer, '_HQinterfaces_mutations.txt'),colClasses="character")
#Pan-cancer
for(cancer in cancer_types[2:33]){
    Mutations <- rbind(Mutations,
                       read.table(paste0('/data/Data/TCGA_Somatic_Mutations/Mutect2/TCGA.', cancer, '.mutect.NonSilent.maf'),
                       skip=1,header=TRUE,colClasses="character",sep="\t"))
    PPI_mutations <- rbind(PPI_mutations,
                           read.table(paste0('HQinterfaces_mutations/', cancer, '_HQinterfaces_mutations.txt'),
                            colClasses="character"))
}

all_uniprot <- unique(PPI_mutations[,3])
uniprot_gene_mapping <- rep(NA,length(all_uniprot))
names(uniprot_gene_mapping) <- all_uniprot
uniprot_pvalues <- rep(1,length(all_uniprot))
names(uniprot_pvalues) <- all_uniprot
# for each protein, calculate pvalue of PPI mutation
for(uniprot in all_uniprot){
    uniprot_gene_mapping[uniprot] <- PPI_mutations[match(uniprot,PPI_mutations[,3]),2]
    mutation_number <- sum(Mutations[,68]==uniprot)
    PPI_mutation_number <- sum(PPI_mutations[,3]==uniprot)
    gene_len <- uniprot_length[match(uniprot,uniprot_length[,1]),2]
    bind_len <- length(psites[[uniprot]])
    if(PPI_mutation_number > 0 & mutation_number > 0){
    # binominal test
        p <- binom.test(x=PPI_mutation_number,n=mutation_number,p=bind_len/gene_len,alternative="greater")$p.value
        uniprot_pvalues[uniprot] <- p
    }
}
write.table(cbind(uniprot_gene_mapping,uniprot_pvalues),
            file=paste("./uniprot_pvalues/Pan-cancer_uniprot_pvalues.txt",sep=""),
            col.names=FALSE,
            sep="\t",
            quote=FALSE)

Interpretation of results

Authors inspected 1,750,987 missense somatic mutations from 10,861 tumor exomes across 33 cancer types from The Cancer Genome Atlas (TCGA) in the interface regions of
121,575 PPIs. They found a significantly higher somatic mutation burden at PPI interfaces compared to noninterfaces across all 33 cancer types.

References

[1] F Cheng, Zhao J , Wang Y , et al. Comprehensive characterization of protein–protein interactions perturbed by disease mutations[J]. Nature Genetics, 2021, 53(3):1-12.

Original

https://github.com/ChengF-Lab/oncoPPIs

Donation

Introduction

 NicheNet, a method that predicts ligand–target links between interacting cells by combining their expression data with prior knowledge on signaling and gene regulatory networks.
 Authors applied NicheNet analysis to uncover the diverse external stimulus for LAMP3+ cDCs from different origins.

Code explanation

Load the required library.

## /data2/csj/tools/R-3.6.3/bin/R
.libPaths("/data2/csj/tools/Rlib3.6.3")

library(nichenetr)
library(tidyverse)

They first calculated two groups of differently expressed genes (DEGs): cDC1-derived LAMP3+ cDCs versus cDC1s (log2FC > 0.2 and adjusted Pvalue < 0.05) and cDC2-derived LAMP3+ cDCs versus cDC2s (log2FC > 0.15 and adjusted P-value < 0.05).
 
ligand_target_matrix.rds file contains a Ligand-target model, which denotes the prior potential that a particular ligand might regulate the expression of a specific target gene.

ligand_target_matrix = readRDS("ligand_target_matrix.rds")
ligand_target_matrix[1:5,1:5] # target genes in rows, ligands in columns

h5ad <- readRDS("/data2/csj/Pan_Myeloid/A20191105/NicheNet_analysis/cDC_for_NicheNet.rds")

DC1_sig <- readRDS("/data2/csj/Pan_Myeloid/A20191105/NicheNet_analysis/cDC1_maturation_sig.rds")
DC2_sig <- readRDS("/data2/csj/Pan_Myeloid/A20191105/NicheNet_analysis/cDC2_maturation_sig.rds")

DC1_tmp <- DC1_sig[which(DC1_sig$Adj_pval < 0.05 & DC1_sig$log2FC > 0.2 ),]
DC2_tmp <- DC2_sig[which(DC2_sig$Adj_pval < 0.05 & DC2_sig$log2FC > 0.15),]

Then, they compared the two lists of DEGs and generated three groups of maturation associated genes: cDC1-specific, cDC2-specific and shared maturation.

overlapped_genes <- intersect(DC1_tmp$gene,DC2_tmp$gene)
DC1_specific_vector <- DC1_tmp[!(DC1_tmp$gene %in% overlapped_genes),]
DC2_specific_vector <- DC2_tmp[!(DC2_tmp$gene %in% overlapped_genes),]
library(dplyr)

Finally, the top 100 genes ordered by log2FC from each maturation signature were used as gene sets of interest for NicheNet analysis.

DC1_specific_df <- DC1_specific_vector %>%  top_n(n = 100, wt = log2FC)
DC2_specific_df <- DC2_specific_vector %>%  top_n(n = 100, wt = log2FC)
common_sig_df <- DC1_tmp[DC1_tmp$gene %in% overlapped_genes, ] %>%  top_n(n = 100, wt = log2FC)

DC1_specific <- as.vector(DC1_specific_df$gene)
DC2_specific <- as.vector(DC2_specific_df$gene)
common_sig <- as.vector(common_sig_df$gene)

df_genes <- data.frame("cDC1-specific" = DC1_specific, "cDC2-specific" = DC2_specific, "shared" = common_sig)
write.csv(df_genes,"cDC maturation sigantures.csv", quote = F, row.names = F, col.names = F)

DC2_specific <- DC2_specific[DC2_specific!='IGLV2-14']
DC2_specific <- DC2_specific[DC2_specific!="IGKV3-20"]
length(common_sig)
length(DC1_specific)
length(DC2_specific)

NicheNet analysis.
First, Define expressed genes in sender and receiver cell populations. Authors selected expressed genes as receptor gene because they only have DC cells. Then they checker the number of expressed genes, which should a 'reasonable' number of total expressed genes in a cell type, e.g. between 5000-10000 (and not 500 or 20000).

expression = h5ad$expression
sample_info = h5ad$metadata
sample_info$cell = rownames(sample_info)

### Define expressed genes in sender and receiver cell populations
# here we do not consider sender population as we only have DC cells

cDC3_cDC1_ids = sample_info %>% filter(MajorCluster == 'cDC3-cDC1') %>% pull(cell)
## selected expressed genes (We will consider a gene to be expressed when it is expressed in at least 10% of cells in one cluster.)
expressed_genes_receiver = expression[,cDC3_cDC1_ids] %>% apply(1,function(x){sum(x>0)/length(x)})  %>% .[. >= 0.1] %>% names()

# Check the number of expressed genes: should be a 'reasonable' number of total expressed genes in a cell type, e.g. between 5000-10000 (and not 500 or 20000)
length(expressed_genes_receiver)
## [1] 6351

Define the gene set of interst and a backgroud of genes. All expressed genes in cDC1- or cDC2-derived LAMP3+ cDCs were used as background of genes.

geneset_oi = DC1_specific %>% .[. %in% rownames(ligand_target_matrix)] 

length(geneset_oi)
## [1] "SERPINE1" "TGFBI"    "MMP10"    "LAMC2"    "P4HA2"    "PDPN"

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
length(background_expressed_genes)
head(background_expressed_genes)
## [1] "RPS11"   "ELMO2"   "PNMA1"   "MMP2"    "TMEM216" "ERCC5"

Define a set of potential ligands

### Define a set of potential ligands
lr_network = readRDS("lr_network.rds")  #Putative ligand-receptor links were gathered from NicheNet’s ligand-receptor data sources.

# If wanted, users can remove ligand-receptor interactions that were predicted based on protein-protein interactions and only keep ligand-receptor interactions that are described in curated databases. To do this: uncomment following line of code:
# lr_network = lr_network %>% filter(database != "ppi_prediction_go" & database != "ppi_prediction")

ligands = lr_network %>% pull(from) %>% unique()
expressed_ligands = ligands ## here we used all ligands intersect(ligands,expressed_genes_sender)

receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,expressed_genes_receiver)

lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
head(lr_network_expressed)

potential_ligands = lr_network_expressed %>% pull(from) %>% unique()

Perform NicheNet ligand activity analysis on the gene set of interest

ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand)

Infer target genes of top-ranked ligands and visualize in a heatmap

active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()
active_ligand_target_links_df <- na.omit(active_ligand_target_links_df)

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)

order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = active_ligand_target_links_df$target %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()

cDC1_specific = vis_ligand_target %>% make_heatmap_ggplot("Pro cDC1 maturation ligands","cDC1-specific maturation signature", color = "red",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "red", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))

cDC1_specific

Similarly, they analyzed cDC2 and common maturation signature

##########cDC2 maturation signature
cDC3_cDC2_ids = sample_info %>% filter(MajorCluster == 'cDC3-cDC2') %>% pull(cell)
expressed_genes_receiver = expression[,cDC3_cDC2_ids] %>% apply(1,function(x){sum(x>0)/length(x)})  %>% .[. >= 0.1] %>% names()

# Check the number of expressed genes: should be a 'reasonable' number of total expressed genes in a cell type, e.g. between 5000-10000 (and not 500 or 20000)
length(expressed_genes_receiver)
## [1] 6351

### Define the gene set of interest and a background of genes
geneset_oi = DC2_specific %>% .[. %in% rownames(ligand_target_matrix)] 
length(geneset_oi)

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
length(background_expressed_genes)

### Define a set of potential ligands
lr_network = readRDS("lr_network.rds")

ligands = lr_network %>% pull(from) %>% unique()
expressed_ligands = ligands ## here we used all ligands intersect(ligands,expressed_genes_sender)

receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,expressed_genes_receiver)

lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
head(lr_network_expressed)

potential_ligands = lr_network_expressed %>% pull(from) %>% unique()

### Perform NicheNet鈥檚 ligand activity analysis on the gene set of interest
ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand)

### Infer target genes of top-ranked ligands and visualize in a heatmap
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)

order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = active_ligand_target_links_df$target %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets[order_targets %in% rownames(active_ligand_target_links)],order_ligands] %>% t()

cDC2_specific = vis_ligand_target %>% make_heatmap_ggplot("Pro cDC2 maturation ligands","cDC2-specific maturation signature", color = "red",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "red", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))

cDC2_specific

##########common maturation signature
cDC3_cDC2_ids = sample_info %>% filter(MajorCluster == 'cDC3-cDC2') %>% pull(cell)
cDC3_cDC1_ids = sample_info %>% filter(MajorCluster == 'cDC3-cDC1') %>% pull(cell)

expressed_genes_receiver = expression[,c(cDC3_cDC2_ids,cDC3_cDC1_ids)] %>% apply(1,function(x){sum(x>0)/length(x)})  %>% .[. >= 0.1] %>% names()

# Check the number of expressed genes: should be a 'reasonable' number of total expressed genes in a cell type, e.g. between 5000-10000 (and not 500 or 20000)
length(expressed_genes_receiver)
## [1] 6351

### Define the gene set of interest and a background of genes
geneset_oi = common_sig %>% .[. %in% rownames(ligand_target_matrix)] 
length(geneset_oi)

background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
length(background_expressed_genes)

### Define a set of potential ligands
lr_network = readRDS("lr_network.rds")

ligands = lr_network %>% pull(from) %>% unique()
expressed_ligands = ligands ## here we used all ligands intersect(ligands,expressed_genes_sender)

receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,expressed_genes_receiver)

lr_network_expressed = lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) 
head(lr_network_expressed)

potential_ligands = lr_network_expressed %>% pull(from) %>% unique()

### Perform NicheNet ligand activity analysis on the gene set of interest
ligand_activities = predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands)
best_upstream_ligands = ligand_activities %>% top_n(20, pearson) %>% arrange(-pearson) %>% pull(test_ligand)

### Infer target genes of top-ranked ligands and visualize in a heatmap
active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)

order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = active_ligand_target_links_df$target %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets[order_targets %in% rownames(active_ligand_target_links)],order_ligands] %>% t()

common_res = vis_ligand_target %>% make_heatmap_ggplot("Pro cDC maturation ligands","cDC common maturation signature", color = "red",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "red", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))

common_res

pdf(file="Common_maturation.pdf",width=6.99, height=4.96)
common_res
dev.off()

pdf(file="cDC2_specific_maturation.pdf",width=4.28, height=4.96)
cDC2_specific
dev.off()

pdf(file="cDC1_specific_maturation.pdf",width=5.36, height=4.96)
cDC1_specific
dev.off()

Heatmaps showing potential ligands driving the maturation of cDC1s.

Interpretation of results

Authors identified specific ligands for each cDC subset. Of the ligands regulating the maturation of cDC1s, IL-4 and IL12B were significantly enriched, consistent
with previous studies (Maier et al., 2020), whereas IL-15, a cytokine reported to induce the conversion of monocytes tomature DCs (Saikh et al., 2001), was identified from cDC2-specific maturation signatures.

References

Cheng S , Li Z , Gao R , et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[J]. Cell, 2021, 184(3):792-809.e23.

Original

https://github.com/Sijin-ZhangLab/PanMyeloid
https://github.com/saeyslab/nichenetr

Donation

Introduction

RUVseq can conduct a differential expression (DE) analysis that controls for “unwanted variation”, e.g., batch, library preparation, and other nuisance effects, using the between-sample normalization methods proposed. They call this approach RUVSeq for remove unwanted variation from RNA-Seq data
If no genes are known a priori not to be influenced by the covariates of interest, one can obtain a set of “in-silico empirical” negative controls, e.g., least significantly DE genes based on a first-pass DE analysis performed prior to RUVg normalization.

Code exmaple

Installation of RUVseq and edgeR from bioconductor

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# The following initializes usage of Bioc devel
BiocManager::install(version='devel')

BiocManager::install("RUVSeq")
BiocManager::install("edgeR")

Load needed packge and data

library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)

zfGenes data: one gene in a row, one sample in a column
Filter out non-expressed genes, by requiring more than 5 reads in at least two samples for each gene.

filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]

Defined group and saved data in SeqExpressionSet

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
phenoData = data.frame(x, row.names=colnames(filtered)))

Firstly we run DEG, then we get empirical control gene by considering all but the top 5000 genes as ranked by edgeR p-values.

design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]

Now we can estimate the factors of unwanted variation

set2 <- RUVg(set, empirical, k=1)
pData(set2)
## x W_1
## Ctl1 Ctl -0.10879677
## Ctl3 Ctl 0.23066424
## Ctl5 Ctl 0.19926266
## Trt9 Trt 0.07672121
## Trt11 Trt -0.83540924
## Trt13 Trt 0.43755790

Now , We can put W_1 as the factors of unwanted variation combining with group factor into the negative binomial GLM model of edgeR. Then the DE genes without batch effect came out.

design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)  # DE gene order by pvalue

If someone wants to use DESeq2 instead of edgeR, the codes are as follows.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts(set1),
colData = pData(set1),
design = ~ W_1 + x)
dds <- DESeq(dds)
res <- results(dds)
res

References

http://www.bioconductor.org/packages/devel/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.pdf

Introduction

Whether it is single-cell RNA sequencing or bulk RNA sequencing, we need to perform differential gene expression analysis and GSVA pathway enrichment analysis. The following code is from a single-cell article

Code explanation

Load the required library and function. parse_h5ad is used to read anndata from the workflow of Scanpy.As we can see, expression data and meta data are returned.

## load data
.libPaths("/data2/csj/tools/Rlib")
library(Seurat)
library(dplyr)
library(monocle)
options(stringsAsFactors=FALSE)
library(reticulate)

parse_h5ad <- function(adata){
    require(reticulate)
    ad <- import("anndata", convert = FALSE)
    ada <- ad$read_h5ad(adata)
    meta <- py_to_r(ada$obs)
    if(class(ada$raw$X)[1] == "scipy.sparse.csr.csr_matrix"){
        exp <- t(py_to_r(ada$raw$X$toarray()))
    }
    else{
        exp <- t(py_to_r(ada$raw$X))
    }
    rownames(exp) <- rownames(py_to_r(ada$raw$var))
    colnames(exp) <- rownames(meta)
    return(
        list(
        metadata = meta,
        expression = exp
        )
    )
}

Authors calculated log fold change, p value and adjusted p value. "pseudocount.use" was used to avoid zero log. Wilcox.test was used to make test and then p value was adjusted by BH method.

############## different expression
DEGene <- function(h5ad,Cluster1,Cluster2){
    ## h5ad: the result of function parse_h5ad
    ## Cluster1 : a vector of Cluster
    ## Cluster2 : a vector of Cluster

    pseudocount.use =1

    cell_name_1 <- rownames(h5ad$metadata[h5ad$metadata$MajorCluster %in% Cluster1,])
    cell_name_2 <- rownames(h5ad$metadata[h5ad$metadata$MajorCluster %in% Cluster2,])

    Expression_1 <- h5ad$expression[,cell_name_1]
    Expression_2 <- h5ad$expression[,cell_name_2]

    ## log FC
    mean_c1 <- as.data.frame(rowMeans(as.matrix(Expression_1)))
    colnames(mean_c1) <- "mean_c1"
    mean_c2 <- as.data.frame(rowMeans(as.matrix(Expression_2)))
    colnames(mean_c2) <- "mean_c2"
    log2fc <- data.frame(log2fc = log2(mean_c1$mean_c1 + pseudocount.use) - log2(mean_c2$mean_c2 + pseudocount.use))
    rownames(log2fc) <- rownames(mean_c1)
    log2fc$gene <- rownames(log2fc)

    ## wilcox test
    group.info <- data.frame(row.names = c(cell_name_1, cell_name_2))
    group.info[cell_name_1, "group"] <- "Group1"
    group.info[cell_name_2, "group"] <- "Group2"
    group.info[, "group"] <- factor(x = group.info[, "group"])
    data.use <- h5ad$expression[, rownames(x = group.info), drop = FALSE]

    p_val <- sapply(
        X = 1:nrow(x = data.use),
        FUN = function(x) {
        return(wilcox.test(data.use[x, ] ~ group.info[, "group"])$p.value)
    })

    ## BH correction
    adj_p_val <- p.adjust(p_val, method="BH")

    ## DE table
    result <- data.frame(gene=log2fc$gene, log2FC=log2fc$log2fc, Pvalue=p_val, Adj_pval=adj_p_val)
    return(result)
}

The clusterProfiler package was used to run GSVA analysis.

##################################################### GSVA
runGSVA <- function(h5ad,Cluster1,Cluster2,kcdf="Gaussian",AdjPvalueCutoff=0.05){
    ## h5ad: the result of function parse_h5ad
    ## Cluster1 : Cluster1
    ## Cluster2 : Cluster2
    ## kcdf: kcdf="Gaussian" for continuous and 'Poisson for integer counts'

    require(GSVA)
    require(GSEABase)
    require(GSVAdata)
    require(clusterProfiler)
    data(c2BroadSets)
    library(limma)

    expression <- h5ad$expression[,rownames(h5ad$metadata[h5ad$metadata$MajorCluster %in% c(Cluster1,Cluster2),])]

    ## change gene symbol to geneid
    gene_entrezid <- bitr(rownames(expression), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")
    expression_filt <- expression[gene_entrezid$SYMBOL,]
    rownames(expression_filt) <- gene_entrezid$ENTREZID
    expression_filt <- as.matrix(expression_filt)

    res_gsva <- gsva(expression_filt, c2BroadSets, parallel.sz=10,kcdf=kcdf) 

    annotation_col = data.frame(CellType = factor(h5ad$metadata[colnames(res_gsva),]$MajorCluster))

    rownames(annotation_col) = colnames(res_gsva)

    ## using limma to conduct DE analysis
    f <- factor(annotation_col$CellType)
    design <- model.matrix(~0+f)
    colnames(design) <- c("C1","C2")
    rownames(design) <- colnames(res_gsva)

    fit <- lmFit(res_gsva, design)
    cont.matrix=makeContrasts('C1-C2',levels = design)
    fit2=contrasts.fit(fit,cont.matrix)
    fit2=eBayes(fit2)
    gs <- topTable(fit2,adjust='BH', number=Inf, p.value=AdjPvalueCutoff)
    gs$cluster <- ifelse(gs$t > 0 , "C1", "C2")
    return(gs)
}

Bar plot showing different pathways enriched in C1QC+ macrophage and SPP1+ macrophage in lung cancer scored per cell by gene set variation analysis
(GSVA). t values are calculated with limma regression

Interpretation of results

Although there is not plot function of codes, it is easy to run it in bar plot function.
Figure shows that the SPP1+ TAMs showed preferential expression of genes involved in angiogenesis as a result of validation.

References

Cheng S , Li Z , Gao R , et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[J]. Cell, 2021, 184(3):792-809.e23.

Original

https://github.com/Sijin-ZhangLab/PanMyeloid

Introduction

Somethimes we need to replace multiple contents of vector in R , like I want to replace a to c and b to d in x vector . Obviously you can use gsub to finish it by two times : gsub (a , c , x ) ; gsub (b , d , x ) . Actually, we can do it in a more readable way .

Code exmaple

> library(plyr)
> x <- c("a", "b", "f")
> mapvalues(x, c("a", "b"),c("c", "d"))
[1] "c" "d" "f"

Here the old names were changed into new names at one time .

References

The manual of mapvalues function

Original

https://github.com/Sijin-ZhangLab/PanMyeloid

Introduction

By performing a pan-cancer analysis of single myeloid cells , authors found that some mutations were correlated with the fractions of myeloid subset.

Code explanation

Load the required library.

## /share/tools/R/R-3.6.0/bin/R
.libPaths("/data2/csj/tools/Rlib3.6")

suppressPackageStartupMessages({
library("reshape2")
library("plyr")
library("dplyr")
library("ggpubr")
library("data.table")
library("MASS")
library("R.utils")
library("nlme")
#library("MuMIn")
library("forestplot")
library("glmmLasso")
#library("lmmen")
})

Read TMB data, mutation data, cell fraction data. And then merge them by patient name

options(stringsAsFactors=F)
## read TMB
tmb = read.table("/data2/csj/Pan_Myeloid/A20191105/mutation/TMB.merge.txt", header=F, check.names=F, stringsAsFactors=F)
colnames(tmb) = c("patient","cancerType","TMB")
tmb$patient = gsub("-",".",tmb$patient)
tmb$patient = gsub("BC","BRCA",tmb$patient)
tmb$patient = gsub("PACA","PAAD",tmb$patient)
old.name <- c("BC", "CHOL", "CRC", "ESCA", "HCC", "LUNG", "FTC", "OV", "PACA", "RC", "THCA", "UCEC")
new.name <- c("BRCA", "CHOL", "CRC", "ESCA", "HCC", "LUNG", "FTC", "OV", "PAAD", "RC", "THCA", "UCEC")
tmb$cancerType <- plyr::mapvalues(x = tmb$cancerType, from = old.name, to = new.name)
# change cancer type name

## read cellInfo
cellInfo = read.csv("/data2/csj/Pan_Myeloid/A20191105/scanorama_integrate/metadata.csv")
cellInfo = cellInfo[cellInfo$Cluster != 'Myeloid',]

cellInfo = cellInfo[cellInfo$tissue=="T",]
cellInfo = cellInfo[cellInfo$source=="ZhangLab",] 
current.cluster.ids <- c("ESCA", "L", "MM", "OV", "PACA", "RC", "THCA", "UCEC")
new.cluster.ids <- c("ESCA", "LYM", "MM", "OV", "PAAD", "RC", "THCA", "UCEC")
cellInfo$cancer <- plyr::mapvalues(x = cellInfo$cancer, from = current.cluster.ids, to = new.cluster.ids)
cellInfo <- cellInfo[c('cancer','patient','Cluster','MajorCluster','library_id')]
cellInfo$patient <- paste0(cellInfo$cancer,".",cellInfo$patient)

cellInfo$subcluster <- unlist(lapply(strsplit(cellInfo$MajorCluster,"_"),function(x){x[2]}))
old.cluster.name <- c("Mono", "Mast", "Macro", "pDC", "cDC3", "cDC2", "cDC1", "Monolike")
new.cluster.name <- c("Mo/Mq", "Mast", "Mo/Mq", "pDC", "cDC3", "cDC2", "cDC1", "Mo/Mq")
cellInfo$subcluster <- plyr::mapvalues(x = cellInfo$subcluster, from = old.cluster.name, to = new.cluster.name)
colnames(cellInfo) <- c('cancerType','patient','Majorcluster','Cluster','library_id','meta.cluster')

# subset, calculate and then merge tables
cells = cellInfo[,c("patient","meta.cluster")] %>% table %>% as.data.frame()
cells = cells %>% group_by(patient) %>% dplyr::mutate(Percent=100*Freq/sum(Freq))
cells$patient = as.character(cells$patient)
cells$meta.cluster = as.character(cells$meta.cluster)

overlap = intersect(unique(tmb$patient), unique(cells$patient))
length(overlap)  ## 35
tmb = tmb[tmb$patient %in% overlap, ]
cellInfo = cellInfo[cellInfo$patient %in% overlap, ]
cells = cells[cells$patient %in% overlap,]

dat = merge(cells, tmb[,c("patient","cancerType","TMB")], by="patient")

# Spec. Mut.
## load files, make the dataframe
mut = read.table("/data2/csj/Pan_Myeloid/A20191105/mutation/mut_mtx.flt.fmt.txt", sep="\t", header=T, check.names=F, stringsAsFactors=F)
mut$Gene = gsub("-","_", mut$Gene)
colnames(mut) = gsub("-",".",colnames(mut))
colnames(mut) = gsub("PACA","PAAD",colnames(mut))
colnames(mut) = gsub("BC","BRCA",colnames(mut))
rownames(mut) = mut$Gene
#
flag = mut[,overlap]
flag = ifelse(flag==".", 0, 1)
mut$Num2 = rowSums(flag)
mut$Num2r = round(mut$Num2/length(overlap), 2)
#
census = fread("/data2/csj/Pan_Myeloid/A20191105/mutation/Census_all_Wed_Oct_30_07_30_28_2019.GRCh37_COSMICv90.tsv", sep="\t", header=T, stringsAsFactors=F)
census$<code>Gene Symbol</code> = gsub("-","_", census$<code>Gene Symbol</code>)

## choose gene used
# set1: Num2r>=0.1 && in Census table
table(mut$Num2)
genes = mut[mut$Num2r >= 0.1,"Gene"]
genes = genes[genes %in% census$<code>Gene Symbol</code>]
genes = genes[!grepl("[;,]",genes)]
length(genes)
print(genes)

#reshape
mut = mut[genes, overlap]
mut = as.data.frame(t(ifelse(mut==".", 0, 1)))

dat2 = dcast(cells, patient ~ meta.cluster, value.var="Percent")
dat2[is.na(dat2)] = 0
dat2 = merge(unique(dat[,c("patient","cancerType")]), dat2, by="patient")
clusters = colnames(dat2)[3:ncol(dat2)]

mut$patient = rownames(mut)
dat2 = merge(dat2, mut, by="patient")
dat2 = merge(dat2, tmb[,c("patient","TMB")], by="patient")
dat2$cancerType = factor(dat2$cancerType)
## used save lm result
allRes <<- as.data.frame(matrix(NA, nrow=0, ncol=5))
colnames(allRes) = c("Gene", "Cluster", "Estimate", "EstimateCI", "Pvalue")

Now we can model cell proportions and mutations. If we want a lasso model, we have to select a suitable lambda firstly. Here, authors used three folds cross validation to choose lambda, and then use that lambda to build lasso model.

## define functions
## lasso variables selection, lambda choosed using 3-fold cv
runLasso = function(par.dat2, par.tarClu, par.genes){
  # https://stats.stackexchange.com/questions/70753/correct-estimation-of-arguments-for-glmmlasso-function/134559
  set.seed(123)
  N = nrow(par.dat2)
  ind = sample(N,N)
  lambda = seq(500, 0, by=-5)
  family = gaussian(link="identity")
  #
  par.dat2$y = par.dat2[,par.tarClu[1]]
  if (length(par.tarClu) > 1){
    for(i in par.tarClu[2:length(par.tarClu)]){
       par.dat2$y = par.dat2$y + par.dat2[,i]
    }
  }
  formula = as.formula(paste("y ~", paste(c("TMB",par.genes), collapse='+')))
  #
  kk = 3
  nk = floor(N/kk)
  Devianz_ma = matrix(Inf,ncol=kk,nrow=length(lambda))
  ## first fit good starting model
  PQL = glmmPQL(y~1, random=~1|cancerType, family=family, data=par.dat2, verbose=F)
  Delta.start = c(as.numeric(PQL$coef$fixed), rep(0,length(par.genes)+1), as.numeric(t(PQL$coef$random$cancerType)))
  Q.start = as.numeric(VarCorr(PQL)[1,1])
  ## train and test
  for(j in 1:length(lambda)){
    for (i in 1:kk){
      if (i < kk){
        indi = ind[(i-1)*nk+(1:nk)]
      }else{
        indi = ind[((i-1)*nk+1):N]
      }

      par.dat2.train = par.dat2[-indi,]
      par.dat2.test = par.dat2[indi,]

      glm2 = try(glmmLasso(formula, rnd=list(cancerType=~1), family=family, data=par.dat2.train, lambda=lambda[j], 
                  switch.NR=F, final.re=TRUE, control=list(start=Delta.start, q_start=Q.start)), silent=TRUE) 

      if(class(glm2)!="try-error"){  
        y.hat = predict(glm2, par.dat2.test)    
        Devianz_ma[j,i] = sum(family$dev.resids(par.dat2.test$y, y.hat, wt=rep(1,length(y.hat))))
      }
    }
    #print(paste("Iteration ", j, sep=""))
    #print(sum(Devianz_ma[j,]))
  }
  ## get best lambda
  Devianz_vec = apply(Devianz_ma,1,sum)
  opt2 = which.min(Devianz_vec)
  ## final lasso regression
  glm2_final = glmmLasso(formula, rnd=list(cancerType=~1), family=family, data=par.dat2, lambda=lambda[opt2],
                         switch.NR=F, final.re=TRUE, control=list(start=Delta.start,q_start=Q.start))

  return(glm2_final)
}

Function showing lasso model . Content includes coefficient ,p value , and confident interval .

## show LASSO model info
showLasso = function(par.dat2, par.tarClu, par.lm){
  sum = summary(par.lm)

  tab = as.data.frame(sum$coefficients)
  tab = tab[tab$Estimate!=0,]
  tab = tab[!grepl("^(TMB|\\(Intercept\\))$",rownames(tab)),]
  #
  if (nrow(tab)==0){
    print ("No good model")
    return(NULL)
  }
  print(sum)
  #
  tab$Gene = rownames(tab)
  tab = tab[order(tab$Estimate, decreasing=T), ]
  #
  tab$upper = tab$Estimate + 2*tab$StdErr
  tab$lower = tab$Estimate - 2*tab$StdErr
  tab$sig = ifelse(tab$p.value <0.05, "*", " ")
  #tab$sig = ifelse(tab$p.value <0.01, "**", tab$sig)
  tab$Pval = sprintf("%s%s", tab$p.value, tab$sig)
  tab$Info = sprintf("%.2f (%.2f - %.2f)%s", tab$Estimate, tab$lower, tab$upper, tab$sig)
  tab$Estimate = sprintf("%.2f", tab$Estimate)
  tab$StdErr = sprintf("%.2f", tab$StdErr)

  # forestplot
  print(
    forestplot(rbind(c(c("Gene", NA, "Esitimate (95% CI)")), cbind(tab$Gene, NA, tab$Info)),
      c(NA,tab$Estimate), c(NA,tab$lower), c(NA,tab$upper),
      graph.pos=2, xlab="Esitimate", is.summary=c(T, rep(F,nrow(tab))), boxsize=0.15, 
      hrzl_lines=T, new_page=T,
      col=fpColors(box=c("#1c61b6","black"), lines="#1c61b6", zero="gray50"), 
      lwd.ci=2, ci.vertices=T, ci.vertices.height = 0.05,
      txt_gp=fpTxtGp(label=gpar(cex=1), ticks=gpar(cex=0.7), xlab=gpar(cex=1), title=gpar(cex=1)),
      xticks=seq(from=floor(min(c(0,tab$lower))), length.out=10,
                  by=round((ceiling(max(c(0,tab$upper)))-floor(min(c(0,tab$lower))))/10, 2))
    )
  )
  # save and return
  tab2save(par.dat2, par.tarClu, tab, allRes)
  return(tab$Gene)
}

boxplot showing the comparsison of cell proportion between cancer patients with or without mutations

## boxplot for selected genes
loop.box = function(par.dat2, par.tarClu, par.genes){
  for (i in par.genes){
    p = ggboxplot(par.dat2, x=i, y=paste(tarClu,collapse="+"), color=i, add="jitter", palette="jco", outlier.shape=NA) + 
      stat_compare_means(comparisons=list(c("0","1"))) + 
      theme(legend.position="right", plot.title=element_text(size=9, vjust = -1)) +
      ylab("Percent") + xlab("Mutation") 
    print(p)
  }
}

Save information from lasso model to a data frame. And the column names of data frame are Gene, Cluster, Estimate, EstimateCI, Pvalue.

## save res
tab2save = function(par2.dat2, par2.tarClu, par2.tab, par2.allRes){
  # cid
  cid = par2.tarClu
  if (length(par2.tarClu) > 1){
    cid = paste(par2.tarClu, collapse="+")
  }
  # test.data
  test.dat = rowSums(par2.dat2[,par2.tarClu,drop=F])
  # calculate
  tmp = par2.tab
  tmp$Cluster = cid
  tmp$EstimateCI = gsub("\\*$","",tmp$Info)
  tmp$EstimateCI = sub(" - ", " ~ ",tmp$EstimateCI)
  tmp$Pvalue = NA
  for (i in tmp$Gene){
    t = wilcox.test(test.dat[ par2.dat2[,i]==0 ], test.dat[ par2.dat2[,i]==1])
    tmp[tmp$Gene==i,"Pvalue"] = sprintf("%.3f",t$p.value)
  }
  # save
  this.res = tmp[,c("Gene", "Cluster", "Estimate", "EstimateCI", "Pvalue")]
  allRes <<- unique(rbind(par2.allRes, this.res))
}

Finally we can use the defined functions for different cell cluster.

## Lasso
tarClu = c("cDC1")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
print(lm.genes)
#loop.box(dat2, "cDC3", "MUC4")

tarClu = c("cDC2")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
print(lm.genes)

tarClu = c("cDC3")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
pdf("cDC3_MUC4.pdf",width=3,height=3)
loop.box(dat2, tarClu, "MUC4")
dev.off()

### plot with cancer type
# ggboxplot(dat2, x='MUC4', y=paste('cDC3',collapse="+"), add="jitter", palette="jco", outlier.shape=NA,add.params = list(color = 'cancerType'))+ 
      # stat_compare_means(comparisons=list(c("0","1"))) + 
      # theme(legend.position="right", plot.title=element_text(size=9, vjust = -1)) +
      # ylab("The proportion of LAMP3+ DC") + xlab("MUC4 Mutation")

tarClu = c("pDC")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
print(lm.genes)

tarClu = c("Mast")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
print(lm.genes)

tarClu = c("Mo/Mq")
lm.lasso = runLasso(dat2, tarClu, genes)
lm.genes = showLasso(dat2, tarClu, lm.lasso)
print(lm.genes)

Save total result

## Final save result

allRes$Estimate.p = as.numeric(allRes$Estimate)
allRes$Estimate.p = ifelse(allRes$Estimate.p > 5, 5, allRes$Estimate.p)
allRes$Estimate.p = ifelse(allRes$Estimate.p < -5, -5, allRes$Estimate.p)
#
allRes$EstimateCI.p = sub(" ","\n", allRes$EstimateCI)
allRes$EstimateCI.p = gsub(" ","", allRes$EstimateCI.p)
allRes$EstimateCI.p = gsub("(\\d\\.\\d)\\d","\\1", allRes$EstimateCI.p)
#
allRes$sig = ifelse(allRes$Pvalue<0.05, "P<0.05", NA)
allRes$sig = factor(allRes$sig, levels=c("P<0.05"))
sig.color = c("#EE2C2C")
names(sig.color) = c("P<0.05")
#

Heatmap to show coefficient of genes in lasso

pdf("mutation_myeloid_subset_new.pdf", width=5, height=2)
ggplot(allRes, aes(x=Gene, y=Cluster, fill=Estimate.p)) +
                geom_tile(size=3) +
                geom_tile(aes(color=sig),size=1.2, alpha=0) +
                geom_text(aes(label=Estimate), angle=0, color="black", size=3) +
                scale_fill_gradient2( high="#EEAEEE", mid="#FFFFFF", low="#1C86EE", midpoint=0, limits=c(-5,5), name="Coef") +
                scale_colour_manual(name="significance",values=sig.color) +
                theme_bw() + 
                theme(axis.text.x=element_text(size=8, vjust=1, hjust=0, angle=300)) +
                xlab('') + ylab("cellType")+theme(axis.text.x = element_text(colour="black"), axis.text.y = element_text(colour="black"))

dev.off()               

The following results are the same, the gene is changed from mutated genes to a pathways

########################################################################################################3
## pathway level regression
# Spec. Mut.
## load files, make the dataframe
mut = read.table("/data2/csj/Pan_Myeloid/A20191105/mutation/hallmark_mtx.txt", sep="\t", header=T, check.names=F, stringsAsFactors=F)
colnames(mut) = gsub("-",".",colnames(mut))
colnames(mut) = gsub("PACA","PAAD",colnames(mut))
colnames(mut) = gsub("BC","BRCA",colnames(mut))
colnames(mut) = gsub("PATH","Gene",colnames(mut))
rownames(mut) = mut$Gene
#
flag = mut[,overlap]
flag = ifelse(flag==".", 0, 1)
mut$Num2 = rowSums(flag)
mut$Num2r = round(mut$Num2/length(overlap), 2)

## choose pathway used
# set1: Num2r>=0.1 
table(mut$Num2)
pathways = mut[mut$Num2r >= 0.1,"Gene"]

length(pathways)
print(pathways)

#reshape
mut = mut[pathways, overlap]
mut = as.data.frame(t(ifelse(mut==".", 0, 1)))

dat2 = dcast(cells, patient ~ meta.cluster, value.var="Percent")
dat2[is.na(dat2)] = 0
dat2 = merge(unique(dat[,c("patient","cancerType")]), dat2, by="patient")
clusters = colnames(dat2)[3:ncol(dat2)]

mut$patient = rownames(mut)
dat2 = merge(dat2, mut, by="patient")
dat2 = merge(dat2, tmb[,c("patient","TMB")], by="patient")
dat2$cancerType = factor(dat2$cancerType)

## used save lm result
allRes <<- as.data.frame(matrix(NA, nrow=0, ncol=5))
colnames(allRes) = c("Gene", "Cluster", "Estimate", "EstimateCI", "Pvalue") 

## Lasso
tarClu = c("cDC1")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)

tarClu = c("cDC2")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)

tarClu = c("cDC3")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)
pdf("Path1.pdf",width=4,height=3)
loop.box(dat2, tarClu, "HALLMARK_APOPTOSIS")
dev.off()
pdf("Path2.pdf",width=5,height=3)
loop.box(dat2, tarClu, "HALLMARK_INFLAMMATORY_RESPONSE")
dev.off()

### plot with cancer type
# ggboxplot(dat2, x='HALLMARK_APOPTOSIS', y=paste('cDC3',collapse="+"), add="jitter", palette="jco", outlier.shape=NA,add.params = list(color = 'cancerType'))+ 
      # stat_compare_means(comparisons=list(c("0","1"))) + 
      # theme(legend.position="right", plot.title=element_text(size=9, vjust = -1)) +
      # ylab("The proportion of LAMP3+ DC") + xlab("Apoptosis Mutation")

# ggboxplot(dat2, x='HALLMARK_INFLAMMATORY_RESPONSE', y=paste('cDC3',collapse="+"), add="jitter", palette="jco", outlier.shape=NA,add.params = list(color = 'cancerType'))+ 
      # stat_compare_means(comparisons=list(c("0","1"))) + 
      # theme(legend.position="right", plot.title=element_text(size=9, vjust = -1)) +
      # ylab("The proportion of LAMP3+ DC") + xlab("Inflammatory Response Mutation")

tarClu = c("pDC")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)

tarClu = c("Mast")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)

tarClu = c("Mo/Mq")
lm.lasso = runLasso(dat2, tarClu, pathways)
lm.pathways = showLasso(dat2, tarClu, lm.lasso)
print(lm.pathways)

## Final save result

allRes$Estimate.p = as.numeric(allRes$Estimate)
allRes$Estimate.p = ifelse(allRes$Estimate.p > 5, 5, allRes$Estimate.p)
allRes$Estimate.p = ifelse(allRes$Estimate.p < -5, -5, allRes$Estimate.p)
#
allRes$EstimateCI.p = sub(" ","\n", allRes$EstimateCI)
allRes$EstimateCI.p = gsub(" ","", allRes$EstimateCI.p)
allRes$EstimateCI.p = gsub("(\\d\\.\\d)\\d","\\1", allRes$EstimateCI.p)
#
allRes$sig = ifelse(allRes$Pvalue<0.05, "P<0.05", NA)
allRes$sig = factor(allRes$sig, levels=c("P<0.05"))
sig.color = c("#EE2C2C")
names(sig.color) = c("P<0.05")
#
allRes$Gene = gsub("HALLMARK_","",allRes$Gene)

pdf("pathway_myeloid_subset_new.pdf",width=5,height=2)
ggplot(allRes, aes(x=Gene, y=Cluster, fill=Estimate.p)) +
                geom_tile(size=3) +
                geom_tile(aes(color=sig),size=1.2, alpha=0) +
                geom_text(aes(label=Estimate), angle=0, color="black", size=3) +
                scale_fill_gradient2( high="#EEAEEE", mid="#FFFFFF", low="#1C86EE", midpoint=0, limits=c(-5,5), name="Coef") +
                scale_colour_manual(name="significance",values=sig.color) +
                theme_bw() + 
                theme(axis.text.x=element_text(size=8, vjust=1, hjust=0, angle=340)) +
                xlab('') + ylab("cellType")+theme(axis.text.x = element_text(colour="black"), axis.text.y = element_text(colour="black"))

dev.off()

Heatmap showing the coefficients of genes with somatic mutations significantly associated with the fractions of myeloid subset. Lasso regression.

Boxplot showing the comparison of LAMP3+ cDCs fraction between cancer patients with or without mutations in MUC4. Unpaired two-sided Wilcoxon test.

Heatmap showing the coefficients of pathways with somatic mutations significantly associated with the fractions of myeloid subset. Lasso regression.

Boxplot showing the comparison of LAMP3+ cDCs fraction between cancer patients with or without mutations in apoptosis and inflammatory response
pathways. Unpaired two-sided Wilcoxon test.

Interpretation of results

Authors identified several genes with somatic mutations
correlated with myeloid cell subset proportions. Of
note, they identified a negative association between the proportion of LAMP3+ cDCs and mutations in MUC4, a
gene known to be associated with tumor progression via modulating tumor cell apoptosis.

References

Cheng S , Li Z , Gao R , et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[J]. Cell, 2021, 184(3):792-809.e23.

Original

https://github.com/Sijin-ZhangLab/PanMyeloid

Introduction

  In this paper, by performing a pan-cancer analysis of single myeloid cells from 210 patients across 15 human cancer types, authors found that tumor-derived mast cells showed low expression of TNF in multiple cancer types. Apart from TNF, mast-derived VEGFA is a key gene associated with tumor angiogenesis.So they use ROGUE to assessing the herterogeneity of mast cells cluster in cancer.
  ROGUE： Often, it is not even clear whether a given cluster is uniform in unsupervised scRNA-seq data analyses. Here, we proposed the concept of cluster purity and introduced a conceptually novel statistic, named ROGUE, to examine whether a given cluster is a pure cell population.

Code explanation

Load the required library and tool function. Single cell data in this research were processed for dimension reduction and unsupervised clustering by following the workflow in Scanpy, so the result of Scanpy(h5ad files) were read to run ROGUE.

.libPaths("/data2/csj/tools/Rlib3.6")
suppressMessages(library(ROGUE))
suppressMessages(library(ggplot2))
suppressMessages(library(tidyverse))

library(dplyr)
library(monocle)
options(stringsAsFactors=FALSE)
library(reticulate)

parse_h5ad <- function(adata){
    require(reticulate)
    ad <- import("anndata", convert = FALSE)
    ada <- ad$read_h5ad(adata)
    meta <- py_to_r(ada$obs)
    if(class(ada$raw$X)[1] == "scipy.sparse.csr.csr_matrix"){
        exp <- t(py_to_r(ada$raw$X$toarray()))
    }
    else{
        exp <- t(py_to_r(ada$raw$X))
    }
    rownames(exp) <- rownames(py_to_r(ada$raw$var))
    colnames(exp) <- rownames(meta)
    return(
        list(
        metadata = meta,
        expression = exp
        )
    )
}

read h5ad data. Extracted expresion matrix, cluster name, sample name from h5ad

### using the merged dataset of A20191105
h5ad_label <- parse_h5ad("/data2/csj/Pan_Myeloid/A20191105/scanorama_integrate/All_Myeloid.h5ad")
cancers <- c("ESCA","L", "MM", "OV", "PACA", "RC", "THCA", "UCEC")
meta_merged <- h5ad_label$metadata
meta_merged <- meta_merged[meta_merged$cancer %in% cancers,]
meta_merged <- meta_merged[rownames(meta_merged) %in% rownames(meta),]

expr <- expr[,rownames(meta_merged)]
meta <- meta_merged
meta$M1 <- unlist(lapply(as.vector(meta$MajorCluster),function(x){unlist(strsplit(x,"_"))[2]}))

current.cluster.ids <- c("Monolike","Macro","cDC2","Mast","cDC1","pDC","cDC3", "Mono")
new.cluster.ids <- c("Mono/Mq","Mono/Mq","cDC2","Mast","cDC1","pDC","cDC3", "Mono/Mq")
meta$M_cluster <- plyr::mapvalues(x = meta$M1, from = current.cluster.ids, to = new.cluster.ids)

Run rogue

rogue.res <- rogue(expr, labels = meta$M_cluster, samples = meta$patient, platform = "UMI", span = 0.6)
rogue.boxplot(rogue.res)

Authors used a major cell classes this time.

current.cluster.ids <- c("Monolike","Macro","cDC2","Mast","cDC1","pDC","cDC3", "Mono")
new.cluster.ids <- c("Mono/Mq","Mono/Mq","cDC","Mast","cDC","pDC","cDC", "Mono/Mq")
meta$MM_cluster <- plyr::mapvalues(x = meta$M1, from = current.cluster.ids, to = new.cluster.ids)

And run rogue again

res <- rogue(expr, labels = meta$MM_cluster, samples = meta$patient, platform = "UMI", span = 0.6)
rogue.boxplot(res)

saveRDS(rogue.res,"Sub_ROGUE.rds")
saveRDS(res,"/data2/csj/Pan_Myeloid/Maj_ROGUE.rds")

Plot figure to compare ROGUE index between different clusters.

res <- readRDS("/data2/csj/Pan_Myeloid/Sub_ROGUE.rds")
res$patient <- rownames(res)
df <- melt(res)
df <- df[df$variable %in% c('pDC','cDC1','cDC2','cDC3'),]
# df$variable <- factor(df$variable, levels=c("Mono/Mq","cDC1","cDC2","cDC3","pDC","Mast"))
df$variable <- factor(df$variable, levels=c("pDC","cDC1","cDC2","cDC3"))

ggboxplot(df, x = "variable", y = "value",color = "variable",palette ='Paired',
 add = c("mean_se", "dotplot"))+theme(legend.position='none')+ylab("ROGUE index")+xlab("")

ggboxplot(df, x = "variable", y = "value",color = "variable",palette ='Paired',
          ylab = "ROGUE index", 
          add = "jitter",                              # Add jittered points
          add.params = list(size = 1.5, jitter = 0.2)  # Point size and the amount of jittering
          )+theme(legend.position='none')

Use different color panels, and stat_compare_means function to make test

# Color panel -----------------
  c54 <- c("#1A62AF","#26933B","#F3951B","#E71638")   
 my_comparisons <- list( c("cDC2", "pDC"), c("cDC2", "cDC1"), c("cDC2", "cDC3") ) 
  p <- ggboxplot(df, x = "variable", y = "value",color = "variable", palette = c54,
          ylab = "ROGUE index", 
          add = "jitter",                              # Add jittered points
          add.params = list(size = 1.5, jitter = 0.2)  # Point size and the amount of jittering
          )+theme(legend.position='none')+stat_compare_means(comparisons = my_comparisons)+theme(
        legend.position = "null",
        axis.text.x = element_text(
          size = 12,
          angle = 45,
          hjust = 1
        ))
    pdf("/data2/csj/Pan_Myeloid/A20191105/ROGUE_cDC.pdf",width=3.41,height=3.56)
    p
    dev.off()

library(ggbeeswarm)
ggplot(df, aes_string(x = "variable", y = "value")) +
      theme_bw() + xlab("") + ggtitle("") + ylab("ROGUE index")+
      theme(
        legend.position = "null",
        plot.title = element_text(
          size = 16,
          face = "bold",
          hjust = 0.5
        ),
        text = element_text(size = 10),
        plot.margin = unit(c(1, 1, 1, 1), "char"),
        axis.text.x = element_text(
          size = 12,
          angle = 45,
          hjust = 1
        ),
        axis.text.y = element_text(size = 12),
        axis.title = element_text(size = 15)
      )+ geom_boxplot(
          outlier.alpha = 0.5,
          outlier.color = "#a3a3a3",
          outlier.size = 2
        ) + geom_quasirandom(
          cex = 2,
          aes_string(color="variable"),
          width = 0.25,
          alpha = 1
        ) + scale_color_manual(values = c54)

res <- readRDS("/data2/csj/Pan_Myeloid/Maj_ROGUE.rds")
res$patient <- rownames(res)
df <- melt(res)

my_comparisons <- list( c("Mast", "Mono/Mq"), c("Mast", "cDC"), c("Mast", "pDC") )
p <- ggboxplot(df, x = "variable", y = "value",color = "variable", palette = c('dodgerblue2','green4','#E31A1C','#6A3D9A'),
          ylab = "ROGUE index", 
          add = "jitter",                              # Add jittered points
          add.params = list(size = 1.5, jitter = 0.2)  # Point size and the amount of jittering
          )+theme(legend.position='none')+stat_compare_means(comparisons = my_comparisons)+theme(
        legend.position = "null",
        axis.text.x = element_text(
          size = 12,
          angle = 45,
          hjust = 1
        ))
pdf("/data2/csj/Pan_Myeloid/A20191105/ROGUE_major.pdf",width=3.41,height=3.56)
p
dev.off()

# Color panel -----------------
  c54 <- c('dodgerblue2','green4','#E31A1C','#6A3D9A','#FF7F00',
         '#FB9A99','#CAB2D6','khaki2','deeppink1','blue1',      
         'steelblue4','green1','yellow4','yellow3','forestgreen',
         'red2','orange','cornflowerblue', 'magenta','darkolivegreen4',
         'indianred1','tan4','darkblue','mediumorchid1','firebrick4',
         'yellowgreen','lightsalmon','tan3','tan1','darkgray',
         'wheat4','#DDAD4B','chartreuse','seagreen1','moccasin',
         'mediumvioletred','seagreen','cadetblue1','darkolivegreen1','tan2',
         'tomato3','#7CE3D8','gainsboro','gold1','skyblue2',
         'palegreen2','#FDBF6F','gray70','darkorange4','orchid1',
         'darkturquoise','maroon','brown','black')        
library(ggbeeswarm)
ggplot(df, aes_string(x = "variable", y = "value")) +
      theme_bw() + xlab("") + ggtitle("") + ylab("ROGUE index")+
      theme(
        legend.position = "null",
        plot.title = element_text(
          size = 16,
          face = "bold",
          hjust = 0.5
        ),
        text = element_text(size = 10),
        plot.margin = unit(c(1, 1, 1, 1), "char"),
        axis.text.x = element_text(
          size = 12,
          angle = 45,
          hjust = 1
        ),
        axis.text.y = element_text(size = 12),
        axis.title = element_text(size = 15)
      )+ geom_boxplot(
          outlier.alpha = 0.5,
          outlier.color = "#a3a3a3",
          outlier.size = 1
        ) + geom_quasirandom(
          cex = 1,
          aes_string(color="variable"),
          width = 0.25,
          alpha = 1
        ) + scale_color_manual(values = c54)

Interpretation of results

Mast cells showed the highest homogeneity among the four major lineages, so is is possibly a cluster that cannot be subdivided. Then authors directly compared the expression patterns of TNF and
VEGFA in mast cells and observed both mutually exclusive
and co-expression pattern across different cancer types

References

Cheng S , Li Z , Gao R , et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[J]. Cell, 2021, 184(3):792-809.e23.
ROGUE package: https://github.com/PaulingLiu/ROGUE

Original

https://github.com/Sijin-ZhangLab/PanMyeloid

# set a color for each cancer type
library(RColorBrewer)
cancer_color <- data.frame("type"=sort(unique(readRDS("CCLE_heterogeneity_Rfiles/CCLE_metadata.RDS") $cancer_type_trunc)), "color"= c(brewer.pal(12, "Set3")[c(1:6,8,7,10:11)], "maroon","gray93", "yellow2", "goldenrod1", "slateblue2"),stringsAsFactors = F)

References

Kinker G S , Greenwald A C , Tal R , et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity[J]. Nature Genetics, 2020, 52(11):1208-1218.

Original code

https://github.com/broadinstitute/single_cell_classification

Introduction

In this paper, single-cell 10X RNA sequencing was performed on a mixture of a variety of known cell lines to study the interaction mode between polyclones. Here we introduce the single-cell sequencing gene expression cell classification operation. However, the article uses the known inherent SNP for classification, and the gene expression classification is used for comparison with the SNP classification.

代码讲解

Read in the seurat obj that have passed the QC, and see the previous blog for the filtering code , see the previous blog for the filtering codehttp://www.thecodesearch.com/2021/01/28/%e5%8d%95%e7%bb%86%e8%83%9erna%e6%b5%8b%e5%ba%8fqc%e5%b7%b2%e7%9f%a5%e6%a0%b7%e6%9c%acsnp/

Set the parameters according to the number of sequenced cell lines, n_pcs is the number of principal components in principal component analysis, and clust_res is the resolution parameter of running FindCluster function

n_cls <- nlevels(seuObj)
n_pcs <- 2*n_cls
if(n_cls > 25-param_range & n_cls < 25+param_range) {
  clust_res <- 1
} else if(n_cls > 50-param_range & n_cls < 50+param_range) {
  clust_res <- 2
} else if(n_cls > 100-param_range & n_cls < 100+param_range) {
  clust_res <- 4
} else {
  stop('Not implemented')
}

Select the cells that have been quality-controlled according to the markers in the previous QC

#use just good singlets
cq <- Seurat::FetchData(seuObj, vars = c('cell_quality'))
seuObj <- seuObj[, which(cq$cell_quality == 'normal')]

Perform cell sorting

seuObj <- Seurat::ScaleData(object = seuObj, vars.to.regress = vtr)
seuObj <- Seurat::FindVariableFeatures(object = seuObj,
                                       nfeatures = 1e5,
                                       selection.method = 'vst')
seuObj <- Seurat::RunPCA(object = seuObj,
                         features = Seurat::VariableFeatures(seuObj),
                         seed.use = 1,
                         npcs = n_pcs,
                         verbose = FALSE)

seuObj <- Seurat::RunTSNE(object = seuObj,
                          dims = 1:n_pcs,
                          check_duplicates = FALSE,
                          seed.use = 1,
                          perplexity = 30,
                          verbose = FALSE) 

#cluster cells
seuObj <- Seurat::FindNeighbors(seuObj, reduction = 'pca',
                                dims = 1:n_pcs,
                                k.param = 10, 
                                force.recalc = TRUE,
                                verbose = FALSE)
seuObj <- Seurat::FindClusters(seuObj, resolution = clust_res, verbose = FALSE)

References

Kinker G S , Greenwald A C , Tal R , et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity[J]. Nature Genetics, 2020, 52(11):1208-1218.

Original code

https://github.com/broadinstitute/single_cell_classification