|
| 1 | +library(Seurat) |
| 2 | +library(ggplot2) |
| 3 | +library(optparse) |
| 4 | +library(dplyr) |
| 5 | +library(stringr) |
| 6 | + |
| 7 | +option_list <- list( |
| 8 | + make_option(c("-w", "--workdir"), type='character', action='store', default=NA, |
| 9 | + help="Path to the working directory"), |
| 10 | + make_option(c("-r", "--rdsfiles"), type='character', action='store', default=NA, |
| 11 | + help="List of rds files for the samples to integrate") |
| 12 | +) |
| 13 | + |
| 14 | +opt <- parse_args(OptionParser(option_list=option_list)) |
| 15 | + |
| 16 | + |
| 17 | +seur_list <- sapply(strsplit(opt$rdsfiles, ',')[[1]], readRDS) |
| 18 | +names(seur_list) <- sapply(seur_list, function(x) x$Sample[1]) |
| 19 | + |
| 20 | + |
| 21 | +setwd(opt$workdir) |
| 22 | + |
| 23 | +obj <- merge(seur_list[[1]], seur_list[2:length(seur_list)], add.cell.ids=names(seur_list)) |
| 24 | + |
| 25 | +temp <- obj |
| 26 | + |
| 27 | +obj <- NormalizeData(obj) |
| 28 | +obj <- FindVariableFeatures(obj) |
| 29 | +obj <- ScaleData(obj) |
| 30 | +obj <- RunPCA(obj) |
| 31 | +obj <- FindNeighbors(obj, dims = 1:30, reduction = "pca") |
| 32 | +obj <- FindClusters(obj, resolution = 0.2, cluster.name = "merged_clusters") |
| 33 | +obj <- RunUMAP(obj, dims = 1:30, reduction = "pca", reduction.name = "umap.merged") |
| 34 | + |
| 35 | +options(future.globals.maxSize = 8000 * 1024^2) |
| 36 | + |
| 37 | +for (method in c('CCAIntegration', 'RPCAIntegration', 'HarmonyIntegration', 'JointPCAIntegration')) { |
| 38 | + try({ |
| 39 | + reduction <- paste0("integrated.", gsub('Integration', '', method) %>% tolower()) |
| 40 | + clusters <- paste0(gsub('Integration', '', method) %>% tolower(), '_clusters') |
| 41 | + umap <- gsub('integrated', 'umap', reduction) |
| 42 | + obj <- IntegrateLayers( |
| 43 | + object = obj, method = method, |
| 44 | + orig.reduction = "pca", new.reduction = reduction, |
| 45 | + verbose = FALSE |
| 46 | + ) |
| 47 | + obj <- FindClusters(obj, resolution = 0.2, cluster.name = clusters) |
| 48 | + obj <- RunUMAP(obj, reduction = reduction, dims = 1:30, reduction.name = umap) |
| 49 | + }) |
| 50 | +} |
| 51 | + |
| 52 | +saveRDS(obj, 'integrated_norm.rds') |
| 53 | + |
| 54 | +groupwidth <- 2/3 |
| 55 | +plotwidth <- 19/3 |
| 56 | +columns <- obj$Sample %>% unique %>% length %>% sqrt %>% ceiling |
| 57 | +rows <- ((obj$Sample %>% unique %>% length)/columns) %>% ceiling |
| 58 | +for (reduction in obj@reductions %>% names %>% grep(pattern='umap', value=T)) { |
| 59 | + name <- gsub(pattern='umap.', replacement='', x=reduction) |
| 60 | + groupcols <- (obj[[paste0(name, '_clusters')]][,1] %>% unique %>% length/20) %>% ceiling |
| 61 | + ggsave(filename=paste0('UMAP_Norm_', reduction, '_Sample.png'), plot=DimPlot(obj, reduction=reduction, group.by='Sample', label=T), units='in', width=7, height=7, dpi=300) |
| 62 | + ggsave(filename=paste0('UMAP_Norm_', reduction, '_Sample-Split.png'), plot=DimPlot(obj, reduction=reduction, group.by=paste0(name, '_clusters'), split.by='Sample', label=T, ncol=columns), units='in', width=(columns*plotwidth)+groupwidth*groupcols, height=rows*7, dpi=300) |
| 63 | + ggsave(filename=paste0('UMAP_Norm_', reduction, '_Cluster.png'), plot=DimPlot(obj, reduction=reduction, group.by=paste0(name, '_clusters'), label=T), units='in', width=plotwidth+groupwidth*groupcols, height=7, dpi=300) |
| 64 | +} |
| 65 | + |
| 66 | +obj <- SCTransform(temp) |
| 67 | +obj <- RunPCA(obj, npcs = 30, verbose = F) |
| 68 | +obj <- FindNeighbors(obj, dims = 1:30, reduction = "pca") |
| 69 | +obj <- FindClusters(obj, resolution = 0.2, cluster.name = "merged_clusters") |
| 70 | +obj <- RunUMAP(obj, dims = 1:30, reduction = "pca", reduction.name = "umap.merged") |
| 71 | + |
| 72 | +for (method in c('CCAIntegration', 'RPCAIntegration', 'HarmonyIntegration', 'JointPCAIntegration')) { |
| 73 | + try({ |
| 74 | + reduction <- paste0("integrated.", gsub('Integration', '', method) %>% tolower()) |
| 75 | + clusters <- paste0(gsub('Integration', '', method) %>% tolower(), '_clusters') |
| 76 | + umap <- gsub('integrated', 'umap', reduction) |
| 77 | + obj <- IntegrateLayers( |
| 78 | + object = obj, method = method, |
| 79 | + orig.reduction = "pca", new.reduction = reduction, |
| 80 | + verbose = FALSE, normalization.method='SCT' |
| 81 | + ) |
| 82 | + obj <- FindNeighbors(obj, reduction = reduction, dims = 1:30) |
| 83 | + obj <- FindClusters(obj, resolution = 0.2, cluster.name = clusters) |
| 84 | + obj <- RunUMAP(obj, reduction = reduction, dims = 1:30, reduction.name = umap) |
| 85 | + }) |
| 86 | +} |
| 87 | + |
| 88 | +saveRDS(obj, 'integrated_sct.rds') |
| 89 | + |
| 90 | +columns <- obj$Sample %>% unique %>% length %>% sqrt %>% ceiling |
| 91 | +rows <- ((obj$Sample %>% unique %>% length)/columns) %>% ceiling |
| 92 | +for (reduction in obj@reductions %>% names %>% grep(pattern='umap', value=T)) { |
| 93 | + name <- gsub(pattern='umap.', replacement='', x=reduction) |
| 94 | + groupcols <- (obj[[paste0(name, '_clusters')]][,1] %>% unique %>% length/20) %>% ceiling |
| 95 | + ggsave(filename=paste0('UMAP_SCT_', reduction, '_Sample.png'), plot=DimPlot(obj, reduction=reduction, group.by='Sample', label=T), units='in', width=7, height=7, dpi=300) |
| 96 | + ggsave(filename=paste0('UMAP_SCT_', reduction, '_Sample-Split.png'), plot=DimPlot(obj, reduction=reduction, group.by=paste0(name, '_clusters'), split.by='Sample', label=T, ncol=columns), units='in', width=(columns*plotwidth)+groupwidth*groupcols, height=rows*7, dpi=300) |
| 97 | + ggsave(filename=paste0('UMAP_SCT_', reduction, '_Cluster.png'), plot=DimPlot(obj, reduction=reduction, group.by=paste0(name, '_clusters'), label=T), units='in', width=plotwidth+groupwidth*groupcols, height=7, dpi=300) |
| 98 | +} |
| 99 | + |
| 100 | +writeLines(capture.output(devtools::session_info()), 'seuratIntegrate_sessionInfo.txt') |
| 101 | + |
0 commit comments