sc141939.R

setwd("/Users/robert/Desktop/bioinformatics/ipf_sc141939//")

library(data.table)
library(Seurat)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(RCurl)

install.packages("tidyseurat")
library(tidyseurat)
library(dplyr)
library(tidyr)
library(purrr)
library(magrittr)
library(ggplot2)
library(Seurat)
library(tidyseurat)
library(celldex)
library(RColorBrewer)

#reading in sc data
list.files("filtered_feature_bc_matrix/")

scdata <- Read10X("filtered_feature_bc_matrix/", gene.column = 1)
head(scdata)[1:5,1:5]

data_10 <- CreateSeuratObject(
  counts = scdata,
  project = "sc141939",
  min.cells = 3,
  min.features = 250 
)

data_10

head(data_10@meta.data)
data_10@assays$RNA[1:5, 1:10]

#quality control
  #mt%
  head(data_10@assays$RNA@counts)[1:6, 1:20]
  grep("^MT-", rownames(data_10@assays$RNA@counts), value = TRUE)  
  grep("^RPL", rownames(data_10@assays$RNA@counts), value = TRUE)  
  
  data_10$percent.MT <- PercentageFeatureSet(data_10,
                                             pattern = "^MT-")
  data_10$ribo <- PercentageFeatureSet(data_10,
                                       pattern = "^RPL")
  data_10$smallribo <- PercentageFeatureSet(data_10,
                                       pattern = "^RPS")
  head(data_10)  
  
plot(data_10$percent.MT, data_10$nFeature_RNA)
plot(filt.data$percent.MT, filt.data$nFeature_RNA)

#percent largest gene
  data_10$Percent.Largest.Gene <- apply(
    data_10@assays$RNA@counts,
    2,
    function(x)(100*max(x))/sum(x)
  )

  tail(data_10$Percent.Largest.Gene)
  
#plotting qc metrics
VlnPlot(data_10,
        features = c("nCount_RNA", "nFeature_RNA", "percent.MT", "Percent.Largest.Gene"), ncol = 4)

FeatureScatter(data_10,feature1 = "nFeature_RNA", feature2 = "Percent.Largest.Gene")

#filtering
filt.data <- subset(
  data_10,
  nFeature_RNA > 500 &
    percent.MT < 10 &
    Percent.Largest.Gene < 20
)

data_10
filt.data

#replot
FeatureScatter(filt.data, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

#normalization
filt.data[["RNA"]]@data[1:8,1:20]

filt.data <- NormalizeData(filt.data, normalization.method = "LogNormalize")
head(filt.data@assays$RNA@data)[1:5,1:5]

gene.expression <- apply(filt.data@assays$RNA@data, 1, mean)
gene.expression <- sort(gene.expression, decreasing = TRUE)

head(gene.expression, n = 50)

#gene selection
data.clr <- filt.data

FindVariableFeatures(
  data.clr, 
  selection.method = "vst",
  nfeatures = 1000
) -> data.clr

data.clr

variance.data <- as_tibble(HVFInfo(data.clr), rownames = "Gene")
head(variance.data)


head(VariableFeatures(data.clr))

variance.data %>%
  mutate(hypervariable=Gene %in% VariableFeatures(data.clr)
  ) -> variance.data

head(variance.data, n=10)

  #plot variance vs mean
variance.data %>%
  ggplot(aes(log(mean),log(variance),color=hypervariable)) +
  geom_point() +
  scale_color_manual(values=c("black","red"))

#scaling
data.clr <- ScaleData(data.clr, features = rownames(data.clr))
#after scaling
head(data.clr@assays$RNA@scale.data[1:10,1:10])

#dimensionality reduction
#PCA
data.clr <- RunPCA(data.clr, 
                   features = VariableFeatures(data.clr),
                   npcs = 25)

DimPlot(data.clr, reduction = "pca")
DimPlot(data.clr,reduction="pca", dims=c(3,4))
DimPlot(data.clr,reduction="pca", dims=c(9,10))

  #elbow plot
ElbowPlot(data.clr) #evens out around 19-20 pcs
  #heatmap
DimHeatmap(data.clr, dims = 1:20, cells = 500)
DimHeatmap(data.clr, dims = 21:25, cells = 500)

#tSNE
RunTSNE(
  data.clr,
  dims = 1:20,
  perplexity = 30
) -> data.clr

pdf("tSNE perplexity 30.pdf", width = 10, height = 10)
DimPlot(data.clr, reduction = "tsne", pt.size = 0.5) + ggtitle("tSNE with perplexity 30")
dev.off()

#UMAP
data.clr <- RunUMAP(data.clr, dims = 1:20)

pdf("UMAP 20 dims.pdf", width = 10, height = 10)
DimPlot(data.clr, reduction = "umap", pt.size = 0.5) + ggtitle("UMAP with 20 dim")
dev.off()

#defining cell clusters
data.clr <- FindNeighbors(data.clr, k.param = 20, dims = 1:20)

data.clr
data.clr@graphs$RNA_snn[1:10, 1:10]

data.clr <- FindClusters(data.clr, resolution = 0.8)

head(data.clr@meta.data)
head(data.clr$seurat_clusters, n = 50)

length(data.clr$seurat_clusters[data.clr$seurat_clusters == 4])
cluster4cells <- (data.clr$seurat_clusters[data.clr$seurat_clusters == 4])
cluster5cells <- (data.clr$seurat_clusters[data.clr$seurat_clusters == 5])

data.clr$seurat_clusters[data.clr$seurat_clusters == 5] <- 4

length(data.clr$seurat_clusters[data.clr$seurat_clusters == 4])

saveRDS(data.clr, "data.clr_BRUSH.rds")
DimPlot(data.clr, reduction = "pca", label = T) + ggtitle("PC1 vs PC2 with Clusters")


view(data.clr$seurat_clusters[])

view(data.clr)

DimPlot(data.clr,reduction="tsne", pt.size = 0.5, label = TRUE, label.size = 4)

DimPlot(data.clr,reduction="umap",pt.size = 0.5, label = TRUE, label.size = 4)

#examining cluster properties
VlnPlot(data.clr, features = "nCount_RNA")
VlnPlot(data.clr, features = "nFeature_RNA")
VlnPlot(data.clr, features = "Percent.Largest.Gene")
VlnPlot(data.clr, features = "percent.MT")
VlnPlot(data.clr, features = c("TPSB2", "TPSAB1", "FABP4"))


#finding markers for clusters
unique(Idents(data.clr))

mark_clu0<-FindMarkers(data.clr, ident.1 = 4, min.pct = 0.25)
head(mark_clu0)
VlnPlot(data.clr, features = c("ZMYND10", "WDR54", "ALDH3B1"))

#FIND ALL MARKERS
markers <- FindAllMarkers(data.clr, only.pos = TRUE, min.pct = 0.25,
                                       logfc.threshold = 0.25,slot = "data",test.use ="wilcox",return.thresh = 0.1)


markers[markers$cluster == 1]
head(markers)
nrow(markers)

list.files()
annotations <- read.csv("annotation.csv")
head(annotations)


annotated_markers <- markers %>%
  left_join(y = unique(annotations[, c("geneID", "description")]), by = c("gene" = "geneID"))

head(annotated_markers)

annot_markers <- annotated_markers

top50RNA <- annot_markers %>%
  group_by(cluster) %>% 
  subset(avg_log2FC > 0.7 & p_val_adj < 0.05) %>%
  top_n(n = 50, wt = avg_log2FC)

head(top50RNA)

top50matrix <- matrix(NA, nrow = 50, ncol = 20)
top50matrix[1:1000] <- top50RNA$gene
top50matrix <- as.data.frame(top50matrix)
colnames(top50matrix) = c(seq(1,20))
(top50matrix)

data.clr$orig.ident[1]

levels <- c("BR002", "BR004", "BR005", "BR006", "BR007", "BR008", "BR011", "BR014", "BR015")
levels2 <- c("BR001", "BR003", "BR009", "BR010", "BR012", "BR013")
select <- subset(data.clr, subset = orig.ident %in% levels)
select2 <- subset(data.clr, subset = orig.ident %in% levels2)

select

uip <- select
uip$seurat_clusters <- NULL
uip

uip <- FindNeighbors(uip, k.param=20, dims = 1:20)
uip <- FindClusters(uip, resolution = 0.8)
head(uip$seurat_clusters, n=50)

pdf("uip umap new.pdf", width = 8, height = 8)
DimPlot(uip, reduction = "umap", pt.size = 0.5, label = T, label.size = 5)
dev.off()

data.clr <- FindNeighbors(data.clr, k.param = 20, dims = 1:20)
data.clr <- FindClusters(data.clr, resolution = 0.8)
head(data.clr$seurat_clusters, n = 50)


pdf("uip vs nonuip reductions.pdf", width = 8, height = 8)
DimPlot(select, reduction = "umap", pt.size = 0.5, label = T, label.size = 5)
DimPlot(select, reduction = "tsne", pt.size = 0.5, label = T, label.size = 5)
DimPlot(select2, reduction = "umap", pt.size = 0.5, label = T, label.size = 5)
DimPlot(select2, reduction = "tsne", pt.size = 0.5, label = T, label.size = 5)
dev.off()



FeaturePlot(data.clr, c("KRT17", "PHLDA2"), order = T)
VlnPlot(data.clr, features = "KRT17")
VlnPlot(data.clr, features = "PTPRC")

write.table(top50matrix, "top50RNAgenes.csv", sep = ",")

data.clr <- ScaleData(data.clr, features = top50RNA$gene)

length(data.clr$seurat_clusters[data.clr$seurat_clusters==5])

pdf("heatmaptop50.pdf", width = 50, height = 50)
DoHeatmap(data.clr, features = c(top50RNA$gene), slot = "scale.data", raster = F, draw.lines = T) + theme(axis.text=element_text(size=4))
dev.off()

FeaturePlot(data.clr, 
            c("STPG4", "KRT86", "CAPN12", "NUP210L", "TUBD1", "CDRT1"), 
            max.cutoff = 'q99', min.cutoff = 'q10', 
            order = T, label = TRUE,
            repel = TRUE)

VlnPlot(object = data.clr, log = F,
        features = c("STPG4", "KRT86", "CAPN12", "NUP210L", "TUBD1", "CDRT1"))

top50matrix[1,1]

view(top50matrix)

view(data.clr)
length(data.clr$orig.ident)

data.clr$orig.ident

data.clr$orig.ident[1]


data.clr$cohort <- "UIP"

data.clr

for (i in 1:10058) {
  if (data.clr$orig.ident[i] == "BR001")
    data.clr$cohort[i] = "nonUIP"
}

library(glue)
#visualizing markers for clusters using umap cluster plot
for (i in 1:20) {
  list = c(" ", " ", " ", " ", " ")
  for (j in 1:5) {
    list[j] = top50matrix[j,i]
  }
  pdf(glue("cluster {i-1} marker plot.pdf", width = 8, height = 8))
  print(FeaturePlot(data.clr,
              features = c(list),
              max.cutoff = 'q99', min.cutoff = 'q10',
              order = T, label = T, 
              repel = T))
  dev.off()
}



top50matrix[1,1]
list

#tsne with cluster label
DimPlot(data.clr, reduction = "tsne", pt.size = 0.5, label = T, label.size = 6)

grep("CD", markers$gene, value = T)


#gene markers for each cluster
#cluster 0
FeaturePlot(object = data.clr,
            features = c("CXCL2"),
            order = TRUE,
            min.cutoff = 'q10', max.cutoff = "q99",
            label = TRUE,
            repel = TRUE)
cxcl2_cells <- FindMarkers(data.clr, 
                           ident.1 = 0,
                           ident.2 = 7)

cxcl2_cells <- cxcl2_cells %>%
  rownames_to_column(var = "gene") %>%
  left_join(y = unique(annotations[, c("geneID", "description")]),
            by = c("gene" = "geneID"))

cxcl2_cells <- cxcl2_cells[, c(1, 3:5,2,6:7)]
cxcl2_cells <- cxcl2_cells %>%
  dplyr::arrange(p_val_adj)

head(cxcl2_cells)
write.csv(cxcl2_cells, "cxcl2_cells.csv")
FeaturePlot(object = data.clr,
            features = c("KRT15", "CXCL1", "CXCL2"),
            order = T,
            min.cutoff = 'q20', max.cutoff = 'q99',
            label = T)

top50RNA$pct.1[1]
biggest = 0
counter = 1
store = 0
for (j in 1:20) {
  biggest = 0
  for (i in counter:counter+49) {
    x = (top50RNA$pct.1[i] - top50RNA$pct.2[i])
    if (x > biggest) {
      biggest = x
      store = (top50RNA$gene[i])
    }
  }
  print(store)
  counter = counter + 50
}

biggest = 0
store = ""
for (i in 51:100) {
  x = top50RNA$pct.1[i] - top50RNA$pct.2[i]
  if (x > biggest) {
    biggest = x
    store = top50RNA$gene[i]
  }
}

print(store)

write.csv(top50RNA, "top50RNA_annotated.csv")

view(top50RNA$description)

view(top50matrix)


library(SingleR)

hpca.se <- celldex::HumanPrimaryCellAtlasData()
table(hpca.se$label.main)
table(hpca.se$label.fine)[1:10]
str(hpca.se)

data.clr
filt.data
DimPlot(data.clr, reduction = "tsne", pt.size = 0.5, label = T, label.size = 5)

blueprint_encode <- BlueprintEncodeData()
table(blueprint_encode$label.main)

immune.se2 <- MonacoImmuneData()
table(immune.se2$label.main)

table(immune.se2$label.fine)

immune.se <- DatabaseImmuneCellExpressionData()
table(immune.se$label.main)

ImmGenData.se3 <- ImmGenData()
table(ImmGenData.se3$label.main)

ref = list(BP=blueprint_encode,HPCA=hpca.se, IMM=immune.se2, IMCel=immune.se, IMGen=ImmGenData.se3)
labels = list(blueprint_encode$label.fine,hpca.se$label.fine, immune.se2$label.fine,immune.se$label.fine, ImmGenData.se3$label.fine)

table(labels[1])[1:10]
DefaultAssay(data.clr) <- "RNA"
counts <- GetAssayData(object = data.clr, assay = "RNA", slot = "data")
head(counts)[1:5,1:5]

pred.hesc <- SingleR(test = counts, ref = ref[1], labels = labels[1])
head(pred.hesc)

names(pred.hesc$orig.results)
head(pred.hesc$orig.results$BP)
table(pred.hesc$orig.results$BP$pruned.labels)

types <- pred.hesc$orig.results$BP$labels
view(types)
data.clr <- AddMetaData(data.clr, metadata=types,
                    col.name=paste0("BP","_type" ) )

head(data.clr@meta.data)

pdf("AnnotCells_OA.pdf", width=6, height=6)
par(mar=c(11,4,3,3),mgp= c(3, 0.5, 0)) #c(bottom, left, top, right)
barplot(sort(table(pred.hesc$labels) ,decreasing = T), main="Cell Annotation OA 10X",xlab = "",ylab="Number of cells",col="lightblue",las=2)
dev.off()

DimPlot(data.clr, group.by='BP_type',  reduction="umap",label = T)
DimPlot(data.clr, group.by='BP_type',  reduction="tsne",label = T)

Idents(data.clr)[1:10]

data.clr$celltypes<-Idents(data.clr)
DimPlot(data.clr, group.by='BP_type',  reduction="umap",label = T)
DimPlot(data.clr, group.by='celltypes',  reduction="tsne",label = T)

data.clr$BP_type

data.clr <- RenameIdents(object = data.clr,
                         "0" = "Adipocytes",
                         "1" = "CD4+ T-cells",
                         "2" = "CD4+ Tcm",
                         "3" = "CD4+ Tem",
                         "4" = "CD8+ T-cells",
                         "5" = "CD8+ Tcm",
                         "6" = "CD8+ Tem",
                         "7" = "Chondrocytes",
                         "8" = "CLP",
                         "9" = "CMP",
                         "10" = "DC",
                         "11" = "Epithelial cells",
                         "12" = "GMP",
                         "13" = "HSC",
                         "14" = "Keratinocytes",
                         "15" = "Macrophages",
                         "16" = "Macrophages M1",
                         "17" = "Megakaryocytes",
                         "18" = "Melanocytes",
                         "19")

table(data.clr$BP_type,data.clr$celltypes)#[1:10,1:10]
theme_set(theme_bw(base_size = 9))
cell_pred<-(table(data.clr$BP_type,data.clr$celltypes))
pt <- as.data.frame(cell_pred)

pt$Var1 <- as.character(pt$Var1)
colores_ID=c("#A6CEE3","#1F78B4", "#B2DF8A","#33A02C","orange","#536771","#2d1e60","#7c9c60",
             "#d2ded5","gray","#CC0000", "#006600", "#000000", "#00CCCC","#cce0cc",brewer.pal(9,"Purples")[c(2,4,6,9)],brewer.pal(8,"Dark2"))

pdf("cluster_celltypes.pdf", width = 8, height = 8)
ggplot(pt%>%arrange(Var1), aes(x = Var1, y = Freq, fill = Var2)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  scale_fill_manual(values = c(colores_ID,"#4952fc","red",colores_ID))+
  theme_classic() +
  theme(legend.title = element_blank())+theme(axis.text.x = element_text(angle = 45,vjust = 1, hjust=1, colour = "black"))+
  ggtitle("OA cells")
dev.off()