Mercurial > repos > gaelcge > r_signac_galaxy
diff Signac.R @ 0:6e0b320d8b6a draft default tip
"planemo upload commit dc808171975d0012e25bd7b32adc7a5a5c56a145-dirty"
author | gaelcge |
---|---|
date | Tue, 02 Aug 2022 19:11:27 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Signac.R Tue Aug 02 19:11:27 2022 +0000 @@ -0,0 +1,170 @@ +rm(list = ls()) + +library(Signac) +library(Seurat) +library(GenomeInfoDb) +library(EnsDb.Hsapiens.v75) +library(ggplot2) +library(patchwork) + +library(future) +plan("multicore", workers = (availableCores()-1)) +options(future.globals.maxSize = 3000000 * 1024^2) + +setwd("/home/gaelcge/projects/def-jsjoyal/gaelcge/scATACseq/10XData/atac_v1_pbmc_10k_output/") + +counts <- Read10X_h5(filename = "./atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5") + + +metadata <- read.csv( + file = "./atac_v1_pbmc_10k_singlecell.csv", + header = TRUE, + row.names = 1 +) + +metadata <- metadata[colnames(counts),] + +chrom_assay <- CreateChromatinAssay( + counts = counts, + sep = c(":", "-"), + genome = 'hg19', + fragments = './atac_v1_pbmc_10k_fragments.tsv.gz', + min.cells = 10, + min.features = 200 +) + +pbmc <- CreateSeuratObject( + counts = chrom_assay, + assay = "peaks", + meta.data = metadata +) + +setwd("/home/gaelcge/projects/def-jsjoyal/gaelcge/scATACseq/Signac_analysis/atac_pbmc_500_nextgem") + +pbmc[['peaks']] + +granges(pbmc) + +# extract gene annotations from EnsDb +annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75) + +# change to UCSC style since the data was mapped to GRCh38 +seqlevelsStyle(annotations) <- 'UCSC' +genome(annotations) <- "GRCh38" + +# add the gene information to the object +Annotation(pbmc) <- annotations + + +# compute nucleosome signal score per cell +pbmc <- NucleosomeSignal(object = pbmc) + +# compute TSS enrichment score per cell +pbmc <- TSSEnrichment(object = pbmc, fast = FALSE) + +# add blacklist ratio and fraction of reads in peaks +pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100 +pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments + + + +pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low') + +png("Tssplot.png") +TSSPlot(pbmc, group.by = 'high.tss') + NoLegend() +dev.off() + + + +pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4') +png("FragmentHistogram.png") +FragmentHistogram(object = pbmc, group.by = 'nucleosome_group') +dev.off() + +png("VlnPlot_QC.png", width=1000) +VlnPlot( + object = pbmc, + features = c('pct_reads_in_peaks', 'peak_region_fragments', + 'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'), + pt.size = 0.1, + ncol = 5 +) +dev.off() + + +pbmc <- subset( + x = pbmc, + subset = peak_region_fragments > 3000 & + peak_region_fragments < 20000 & + pct_reads_in_peaks > 15 & + blacklist_ratio < 0.05 & + nucleosome_signal < 2 & + TSS.enrichment > 1 +) +pbmc + + +pbmc <- RunTFIDF(pbmc) +pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0') +pbmc <- RunSVD(pbmc) + +png("DepthCor.png") +DepthCor(pbmc) +dev.off() + +pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30) +pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30) +pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3) + + +png("UMAP.png") +DimPlot(object = pbmc, label = TRUE) + NoLegend() +dev.off() + + +gene.activities <- GeneActivity(pbmc) + + +# add the gene activity matrix to the Seurat object as a new assay and normalize it +pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities) +pbmc <- NormalizeData( + object = pbmc, + assay = 'RNA', + normalization.method = 'LogNormalize', + scale.factor = median(pbmc$nCount_RNA) +) + + +DefaultAssay(pbmc) <- 'RNA' + +png("FeaturePlot_knownMarkers.png", width=1000) +FeaturePlot( + object = pbmc, + features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'), + pt.size = 0.1, + max.cutoff = 'q95', + ncol = 3 +) +dev.off() + + +saveRDS(pbmc, paste("Seurat_object.rds", sep=".")) + + + + + + + + + + + + + + + + + + +q("no")