Tomoliverberry · Tomoliverberry · Sep 16, 2021
diff --git a/code.R b/code.R
@@ -5,13 +5,15 @@ devtools::install_github("GreenleafLab/ArchR", ref="master", repos = BiocManager
 writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")
 Sys.which("make")
 
-setwd("C:/Users/Tom Berry/Desktop/single_cell/")
-SAMPLES <- c("510", "611", "993")
+setwd("C:/Users/Tom/Desktop/single_cell/")
+
 library(ArchR)
 BiocManager::install("BSgenome.Hsapiens.UCSC.hg38", force = TRUE)
 BiocManager::install("Matrix", force = TRUE)
 addArchRGenome("hg38")
 
+SAMPLES <- c("510", "611", "993")
+
 Arrowfiles <- createArrowFiles(
     inputFiles = c(paste0("C:/Users/Tom/Desktop/single_cell/FastFile-3aEnrFqAw78jBfpr/", SAMPLES[1], "_FC_fragments.tsv.gz"), 
                   paste0("C:/Users/Tom/Desktop/single_cell/FastFile-3aEnrFqAw78jBfpr/", SAMPLES[2], "_FC_fragments.tsv.gz"),
@@ -23,15 +25,17 @@ Arrowfiles <- createArrowFiles(
 
 
 project1 <- ArchRProject(
-    ArrowFiles = Arrowfiles, outputDirectory = "singlecell", copyArrows = TRUE
+    ArrowFiles = Arrowfiles, outputDirectory = "C:/Users/Tom/Desktop/single_cell", copyArrows = TRUE
 )
 
 
-saveArchRProject(project, outputDirectory = "save_proj_1", load = FALSE)
+saveArchRProject(project1, outputDirectory = "C:/Users/Tom/Desktop/single_cell", load = FALSE)
 
 #Remove Doublets
 doubletscore <- addDoubletScores(input = project1, k = 10, knnMethod = "UMAP", LSIMethod = 1)
+
 project2 <- filterDoublets(doubletscore)
+
 #Dimentionality Reduction with Iterative LSI
 project2 <- addIterativeLSI(
     ArchRProj = project2,
@@ -78,16 +82,19 @@ project2 <- addClusters(
 
 #Cluster counts batch corrected 
 cluster_counts_harmony <- as.data.frame(t(as.data.frame(as.vector((table(project2$Clusters_Harmony))))))
-rownames(cluster_counts) <- NULL
-colnames(cluster_counts) <- names(table(project2$Clusters_Harmony)) 
+rownames(cluster_counts_harmony) <- NULL
+colnames(cluster_counts_harmony) <- names(table(project2$Clusters_Harmony)) 
 
 cM_nobatchcorrection <- confusionMatrix(paste0(project2$Clusters_no_batch_correction), paste0(project2$Sample))
 cM <- confusionMatrix(paste0(project2$Clusters_Harmony), paste0(project2$Sample))
 
 #RNAseq data read in
+
 RNAseq <- readRDS("seurat.pfc.final.rds")
 RNAseq$cellIDs <- gsub('FC-', '', RNAseq$cellIDs)
 
+#unconstrained integration
+
 project2 <- addGeneIntegrationMatrix(
      ArchRProj = project2,
      useMatrix = "GeneScoreMatrix",
@@ -101,24 +108,30 @@ project2 <- addGeneIntegrationMatrix(
      nameScore = "predictedScore_Un"
  )
 
+
+
 cM <- as.matrix(confusionMatrix(project2$Clusters_Harmony, project2$predictedGroup_Un))
 preClust <- colnames(cM)[apply(cM, 1, which.max)]
 cbind(preClust, rownames(cM))
 
+clust_cM <- pheatmap::pheatmap(
+     mat = as.matrix(cM),
+     color = paletteContinuous("whiteBlue"),
+     border_color = "black", display_numbers = TRUE, number_format = "%.0f"
+ ) 
+
 unique(unique(project2$predictedGroup_Un))
 ExN <- paste0(c("ExN-1","ExN-3","ExN-5","ExN-2","ExN-4","ExN-6"), collapse = "|")
 InN <- paste0(c("InN-1","InN-2","InN-3","InN-4"), collapse = "|")
 RG <- paste0(c("RG-1","RG-2"), collapse = "|")
-MG <- paste0("MG", collapse = "|")
-Else <- paste0(c("IP","CycPro"), collapse="|")
+MG <- paste0("MG")
+
 
 clustExN <- rownames(cM)[grep(ExN, preClust)]
 clustInN <- rownames(cM)[grep(InN, preClust)]
 clustRG <- rownames(cM)[grep(RG, preClust)]
 clustMG <- rownames(cM)[grep(MG, preClust)]
-clustElse <- rownames(cM)[grep(Else, preClust)]
 
-RNA <- RNAseq[grep(clust, RNAseq$cellIDs)]
 RNAcells_ExN <- colnames(RNAseq)[grep(pattern = ExN, x = RNAseq$cellIDs)]
 RNAcells_InN <- colnames(RNAseq)[grep(pattern = InN, x = RNAseq$cellIDs)]
 RNAcells_RG <- colnames(RNAseq)[grep(pattern = RG, x = RNAseq$cellIDs)]
@@ -149,7 +162,11 @@ project2 <- addGeneIntegrationMatrix(
     groupRNA = "cellIDs",
     nameCell = "predictedCell_Co",
     nameGroup = "predictedGroup_Co",
-    nameScore = "predictedScore_Co"
+    nameScore = "predictedScore_Co",
+	dimsToUse = 1:26,
+	k.score = 26,
+	npcs = 26,
+	dims = 1:26
 )
 
 pal <- paletteDiscrete(values = RNAseq$cellIDs)
@@ -191,23 +208,276 @@ cM2 <- confusionMatrix(project3$Clusters_Harmony, project3$predictedGroup)
 labelOld <- rownames(cM2)
 labelNew <- colnames(cM2)[apply(cM2, 1, which.max)]
 
-project3$Clusters_Harmony <- mapLabels(project3$Clusters_Harmony, newLabels = labelNew, oldLabels = labelOld)
+project3$Clusters_RNAmapped <- mapLabels(project3$Clusters_Harmony, newLabels = labelNew, oldLabels = labelOld)
 
-p1 <- plotEmbedding(project3, colorBy = "cellColData", name = "Clusters_Harmony")
+p1 <- plotEmbedding(project3, colorBy = "cellColData", name = "Clusters_RNAmapped")
 
 #making pseudo bulk replicates
 
-project4 <- addGroupCoverages(ArchRProj = project3, groupBy = "Clusters_Harmony")
+project4 <- addGroupCoverages(ArchRProj = project3, groupBy = "Clusters_RNAmapped")
 
 #Calling peaks with tile matrix
 
-project4@projectMetadata$outputDirectory <- "C:/Users/Tom/Desktop/single_cell/save_proj_1"
+project4@projectMetadata$outputDirectory <- "C:/Users/Tom/Desktop/single_cell"
 
 project4 <- addReproduciblePeakSet(
     ArchRProj = project4, 
-    groupBy = "Clusters_Harmony",
+    groupBy = "Clusters_RNAmapped",
     peakMethod = "Tiles",
-    method = "p",
-    cutOff = 0.01,
-    extendSummits = 500
+    method = "p"
+)
+
+getPeakSet(project4)
+
+------------------------------------------------------------------------------------------------------------------------------
+TABLES AND FIGURES
+------------------------------------------------------------------------------------------------------------------------------
+
+#PRE-CLUSTERING
+
+#Exported plots as PNGs through RStudio
+
+#Read in pre-filter quality control data generated by ArchR
+
+`510-Pre-Filter-Metadata` <- readRDS("~/QualityControl/510/510-Pre-Filter-Metadata.rds")
+`611-Pre-Filter-Metadata` <- readRDS("~/QualityControl/611/611-Pre-Filter-Metadata.rds")
+`993-Pre-Filter-Metadata` <- readRDS("~/QualityControl/993/993-Pre-Filter-Metadata.rds")
+
+#Table 1: Preclustering QC summary for donors---------------------------------------------------------------------------------
+
+table_precluster_A <- data.frame(
+"Sample" = "510",
+"Cells_Passed_QC" = sum(`510-Pre-Filter-Metadata`$Keep),
+"Cells_Failed_QC" = sum(`510-Pre-Filter-Metadata`$Keep == 0),
+"Total_Frags" = sum(`510-Pre-Filter-Metadata`$nFrags),
+"Median_Frags" = median(`510-Pre-Filter-Metadata`$nFrags[`510-Pre-Filter-Metadata`$Keep ==1]),
+"Median_TSS_Enrichment" = median(`510-Pre-Filter-Metadata`$TSSEnrichment[`510-Pre-Filter-Metadata`$Keep == 1])
 )
+
+newrow_precluster_A <- data.frame(
+"Sample" = "611",
+"Cells_Passed_QC" = sum(`611-Pre-Filter-Metadata`$Keep),
+"Cells_Failed_QC" = sum(`611-Pre-Filter-Metadata`$Keep == 0),
+"Total_Frags" = sum(`611-Pre-Filter-Metadata`$nFrags),
+"Median_Frags" = median(`611-Pre-Filter-Metadata`$nFrags[`611-Pre-Filter-Metadata`$Keep ==1]),
+"Median_TSS_Enrichment" = median(`611-Pre-Filter-Metadata`$TSSEnrichment[`611-Pre-Filter-Metadata`$Keep == 1])
+)
+table_precluster_A <- rbind(table_precluster_A, newrow_precluster_A)
+
+newrow_precluster_A <- data.frame(
+"Sample" = "993",
+"Cells_Passed_QC" = sum(`993-Pre-Filter-Metadata`$Keep),
+"Cells_Failed_QC" = sum(`993-Pre-Filter-Metadata`$Keep == 0),
+"Total_Frags" = sum(`993-Pre-Filter-Metadata`$nFrags),
+"Median_Frags" = median(`993-Pre-Filter-Metadata`$nFrags[`993-Pre-Filter-Metadata`$Keep ==1]),
+"Median_TSS_Enrichment" = median(`993-Pre-Filter-Metadata`$TSSEnrichment[`993-Pre-Filter-Metadata`$Keep == 1])
+)
+table_precluster_A <- rbind(table_precluster_A, newrow_precluster_A)
+
+#Table 2: Count table for doublet removal-------------------------------------------------------------------------------------
+
+`510-Doublet-Summary` <- readRDS("C:/Users/Tom/Desktop/single_cell/QualityControl/510/510-Doublet-Summary.rds")
+`611-Doublet-Summary` <- readRDS("C:/Users/Tom/Desktop/single_cell/QualityControl/510/510-Doublet-Summary.rds")
+`993-Doublet-Summary` <- readRDS("C:/Users/Tom/Desktop/single_cell/QualityControl/510/510-Doublet-Summary.rds")
+
+table_doubletremoval_B <- data.frame(
+"Sample" = c("510", "611", "993"),
+"Cell_Count_Pre_Doublet_Removal" = c("8321","5029","2287"),
+"Cell_Count_Post_Doublet_Removal" = c("7629","5029","2235"),
+"Percentage_of_Cells_Removed_per_Donor" = c("8.3%","0%","2.3%")
+)
+
+#Plot 1: TSS Enrichment vs log10(number of unique fragments)------------------------------------------------------------------
+
+df <- getCellColData(project2, select = c("log10(nFrags)", "TSSEnrichment"))
+
+df611 <- df[grep("611", rownames(df)), ]
+df510 <- df[grep("510", rownames(df)), ]
+df993 <- df[grep("993", rownames(df)), ]
+
+plot611_1 <- ggPoint(
+     x = df611[,1], 
+     y = df611[,2], 
+     colorDensity = TRUE,
+     continuousSet = "sambaNight",
+     xlabel = "Log10 Unique Fragments",
+     ylabel = "TSS Enrichment",
+     xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
+     ylim = c(0, quantile(df[,2], probs = 0.99))
+ ) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")
+
+ plot510_1 <- ggPoint(
+    x = df510[,1], 
+    y = df510[,2], 
+    colorDensity = TRUE,
+    continuousSet = "sambaNight",
+    xlabel = "Log10 Unique Fragments",
+    ylabel = "TSS Enrichment",
+    xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
+    ylim = c(0, quantile(df[,2], probs = 0.99))
+) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")
+
+plot993_1 <- ggPoint(
+    x = df993[,1], 
+    y = df993[,2], 
+    colorDensity = TRUE,
+    continuousSet = "sambaNight",
+    xlabel = "Log10 Unique Fragments",
+    ylabel = "TSS Enrichment",
+    xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
+    ylim = c(0, quantile(df[,2], probs = 0.99))
+) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")
+
+#Ridge plot for TSS enrichment per donor--------------------------------------------------------------------------------------
+
+#Batch corrected data used for pre clustering with filtering
+
+Ridge_plot <- plotGroups(
+     ArchRProj = project2,
+     groupBy = "Sample",
+     colorBy = "cellColData",
+     name = "TSSEnrichment",
+     plotAs = "ridges"
+ )
+
+#Fragment size plot-----------------------------------------------------------------------------------------------------------
+
+plot_fragmentsizes <- plotFragmentSizes(project2, groupBy = "Sample")
+
+#POST-CLUSTERING--------------------------------------------------------------------------------------------------------------
+
+#Cell counts per cluster per donor--------------------------------------------------------------------------------------------
+
+table_cell_per_cluster <- as.matrix(confusionMatrix(project2$Sample, project2$Clusters_Harmony))
+
+#Cell counts per cluster------------------------------------------------------------------------------------------------------
+
+#Janitor package will be used for data frame editing
+
+install.packages("janitor")
+library(janitor)
+
+#Make data frame copy of table containing number of cells per cluster per donor
+
+table_cell_cluster_all <- as.data.frame(table_cell_per_cluster)
+
+#Move rownames to column as use of adorn_totals creates a "Total" cell which shifts data
+
+setDT(table_cell_cluster_all, keep.rownames = TRUE)[]
+
+#Use janitor feature to apply new row to data frame containing totals
+
+table_cell_cluster_all %<>% adorn_totals(dat = table_cell_cluster_all, where = "row")
+
+#Remove rows containing donor specific data and column containing "Total" cell
+
+table_cell_cluster_all <- table_cell_cluster_all[-c(1,2,3), ]
+table_cell_cluster_all <- table_cell_cluster_all[,-1]
+
+#UMAPs------------------------------------------------------------------------------------------------------------------------
+#IterativeLSI UMAP------------------------------------------------------------------------------------------------------------
+project2 <- addUMAP(
+    ArchRProj = project2,
+    reducedDims = "IterativeLSI",
+    name = "UMAP",
+    nNeighbors = 30,
+    minDist = 0.5,
+    metric = "cosine"
+)
+plot_UMAP_by_cluster <- plotEmbedding(ArchRProj = project2, colorBy = "cellColData", name = "Clusters_no_batch_correction", embedding = "UMAP")
+plot_UMAP_by_donor <- plotEmbedding(ArchRProj = project2, colorBy = "cellColData", name = "Sample", embedding = "UMAP")
+
+#Harmony UMAP-----------------------------------------------------------------------------------------------------------------
+
+project2 <- addUMAP(
+    ArchRProj = project2,
+    reducedDims = "Harmony",
+    name = "UMAP_Harmony",
+    nNeighbors = 30,
+    minDist = 0.5,
+    metric = "cosine"
+)
+plot_UMAP_by_cluster_harmony <- plotEmbedding(ArchRProj = project2, colorBy = "cellColData", name = "Clusters_Harmony", embedding = "UMAP_Harmony")
+plot_UMAP_by_donor_harmony <- plotEmbedding(ArchRProj = project2, colorBy = "cellColData", name = "Sample", embedding = "UMAP_Harmony")
+
+
+
+
+#Final RNA cell type mapping to atac-seq data
+
+ atac_and_rna_p3 <- plotEmbedding(
+     project4, 
+     colorBy = "cellColData",
+     embedding = "UMAP_Harmony",
+     name = "predictedGroup", 
+     pal = pal
+ )
+
+table_cell_cluster_RNA_ATAC <- as.matrix(confusionMatrix(project4$Sample, project4$Clusters_RNAmapped))
+
+
+------------------------------------------------------------------------------------------------------------------------------
+FINEMAPPING
+------------------------------------------------------------------------------------------------------------------------------
+
+#Load in rsnps package, to obtain hg38 information from hg19 data by accessing OpenSNP and NCBI's dbSNP SNP database
+install.packages("rsnps")
+library(rsnps)
+
+#Load readxl package for PGC3 SNP data
+install.packages("readxl")
+library(readxl)
+
+#Read .xlsx file containing PGC3 data
+SNPs_PGC3 <- read_excel("Copy of PGC3 revised 11-8-21 finemapped SNPs with at least 0.1 posterior probability(1335).xlsx")
+
+#rsnps requires the RSID's to be in a character string. Extract RSID's to new vector
+RSID <- as.list(SNPs_PGC3$rsid)
+
+#Perform NCBI search of all RSID's. Package limits to one SNP per second.
+snp_query <- ncbi_snp_query(RSID)
+
+#Read in all peakcalling data from ArchR for each RNA-seq cell type cluster
+peak_ExN_2 <- readRDS(file = "~/PeakCalls/ExN.2-reproduciblePeaks.gr.rds")
+peak_ExN_3 <- readRDS(file = "~/PeakCalls/ExN.3-reproduciblePeaks.gr.rds")
+peak_ExN_4 <- readRDS(file = "~/PeakCalls/ExN.4-reproduciblePeaks.gr.rds")
+peak_ExN_5 <- readRDS(file = "~/PeakCalls/ExN.5-reproduciblePeaks.gr.rds")
+peak_ExN_6 <- readRDS(file = "~/PeakCalls/ExN.6-reproduciblePeaks.gr.rds")
+peak_InN_1 <- readRDS(file = "~/PeakCalls/InN.1-reproduciblePeaks.gr.rds")
+peak_InN_2 <- readRDS(file = "~/PeakCalls/InN.2-reproduciblePeaks.gr.rds")
+peak_InN_3 <- readRDS(file = "~/PeakCalls/InN.3-reproduciblePeaks.gr.rds")
+peak_RG_1 <- readRDS(file = "~/PeakCalls/RG.1-reproduciblePeaks.gr.rds")
+peak_RG_2 <- readRDS(file = "~/PeakCalls/RG.2-reproduciblePeaks.gr.rds")
+
+
+SNPs_h19_h38 <- SNPs_PGC3
+SNPs_h19_h38 <- cbind(SNPs_h19_h38, hg38_position = snp_query$bp)
+SNPs_h19_h38 <- cbind(SNPs_h19_h38, gene = snp_query$gene)
+
+rsID_in_peak_range <- cbind(SNPs_h19_h38[FALSE,], "Cell_type" = character(), "Peak_Start" = integer(), "Peak_Stop" = integer())
+rsID_not_in_peak_range <- cbind(SNPs_h19_h38[FALSE,], "Cell_type" = character(), "Peak_Start" = integer(), "Peak_Stop" = integer())
+
+data_names <- c("ExN_2", "ExN_3", "ExN_4", "ExN_5", "ExN_6", "InN_1", "InN_2", "InN_3", "RG_1", "RG_2")
+
+ for (k in data_names) {
+     peak_CELLTYPE <- eval(parse(text = paste0("peak_",k,"@ranges@start")))
+     for(i in 1:nrow(SNPs_h19_h38)) {
+         break_check = 0
+         for(j in peak_CELLTYPE) {
+             if(j > SNPs_h19_h38[i, 21]){
+                 #print(SNPs_h19_h38[i, ])
+                 rsID_not_in_peak_range <- add_row(rsID_not_in_peak_range, SNPs_h19_h38[i, ], "Cell_type" = k, "Peak_Start" = j, "Peak_Stop" = (j+500))	
+                 break_check = 1
+                 break
+             } else if (SNPs_h19_h38[i, 21] <= j+500){
+                 rsID_in_peak_range <- add_row(rsID_in_peak_range, SNPs_h19_h38[i, ], "Cell_type" = k, "Peak_Start" = j, "Peak_Stop" = (j+500))	
+                 break_check = 1
+                 break
+             }
+         }#end of j
+         if(break_check == 0){
+             rsID_not_in_peak_range <- add_row(rsID_not_in_peak_range, SNPs_h19_h38[i, ], "Cell_type" = k, "Peak_Start" = j, "Peak_Stop" = (j+500))	
+         }
+     } #end of i
+ }#end k
+-------------------------------------------------------------------------------------------------------------------------------