<jupyter_start><jupyter_text>Installation<jupyter_code># if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("chromVAR", version = "3.8")
# BiocManager::install("motifmatchr", version = "3.8")
# BiocManager::install("BSgenome.Hsapiens.UCSC.hg19", version = "3.8")
# BiocManager::install("JASPAR2016", version = "3.8")<jupyter_output><empty_output><jupyter_text>Import packages<jupyter_code>library(chromVAR)
library(motifmatchr)
library(Matrix)
library(SummarizedExperiment)
library(BiocParallel)
library('JASPAR2016')
library(BSgenome.Hsapiens.UCSC.hg19)
packageVersion("chromVAR")
register(MulticoreParam(10))<jupyter_output><empty_output><jupyter_text>Obtain Feature Matrix<jupyter_code>start_time <- Sys.time()
set.seed(2019)
metadata <- read.table('../../input/metadata.tsv',
                         header = TRUE,
                         stringsAsFactors=FALSE,quote="",row.names=1)
head(metadata)
peakfile <- "../../input/atac_v1_pbmc_5k_peaks.bed"
peaks <- getPeaks(peakfile, sort_peaks = TRUE)<jupyter_output>Warning message in getPeaks(peakfile, sort_peaks = TRUE):
“Peaks are not equal width!Use resize(peaks, width = x, fix = "center") to make peaks equal in size, where x is the desired size of the peaks)”<jupyter_text>`width=500` will cause the error: Error in loadFUN(x, seqname, ranges): trying to load regions beyond the boundaries of non-circular sequence "chr17"<jupyter_code>peaks <- resize(peaks, width = 450, fix = "center")
bamfile <- list.files(path = "../../input/sc-bams_nodup/", pattern = "\\.bam$")
length(bamfile)
head(bamfile)
cellnames <- sapply(strsplit(bamfile,'.',fixed = TRUE), "[[", 2)
sum(cellnames == rownames(metadata))
# ix = match(rownames(metadata),cellnames)
# bamfile = bamfile[ix]
# cellnames = cellnames[ix]
# sum(cellnames == rownames(metadata))
head(paste0("../../input/sc-bams_nodup/",bamfile))
fragment_counts <- getCounts(paste0("../../input/sc-bams_nodup/",bamfile), 
                             peaks, 
                             paired = TRUE, 
                             by_rg = FALSE, 
                             format = "bam", 
                             colData = data.frame(cellid = cellnames,
                                                  label = metadata$label))
fragment_counts
fragment_counts <- addGCBias(fragment_counts, genome = BSgenome.Hsapiens.UCSC.hg19)
head(rowData(fragment_counts))
counts_filtered <- filterPeaks(fragment_counts, non_overlapping = TRUE)
bg <- getBackgroundPeaks(counts_filtered)
# Potentially save the bg object
saveRDS(bg, file = "background_peaks_motifs.rds")
motifs <- getJasparMotifs()
motif_ix <- matchMotifs(motifs, counts_filtered, genome = BSgenome.Hsapiens.UCSC.hg19)
dev <- computeDeviations(object = counts_filtered, annotations = motif_ix,
                         background_peaks = bg)
end_time <- Sys.time()
end_time - start_time
df_zscores = dev@assays[[1]]
dim(df_zscores)
df_zscores[1:5,1:5]
sum(sapply(strsplit(colnames(df_zscores),'.',fixed = TRUE), "[[", 2) == rownames(metadata))
colnames(df_zscores) = rownames(metadata)
df_zscores[1:5,1:5]
saveRDS(df_zscores, file = '../../output/feature_matrices/FM_ChromVAR_10xpbmc5k_motifs.rds')<jupyter_output><empty_output><jupyter_text>Downstream Analysis<jupyter_code>variability <- computeVariability(dev)
plotVariability(variability, use_plotly = FALSE)
head(variability)
dim(variability)
dev$label <- as.factor(metadata$label)
tsne_results <- deviationsTsne(dev, threshold = 1.5, perplexity = 10)
tsne_plots <- plotDeviationsTsne(dev, tsne_results, 
                                 sample_column = "label", 
                                 shiny = FALSE)
tsne_plots[[1]]
sessionInfo()
save.image(file = 'ChromVAR_10xpbmc5k_motifs.RData')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1389545573942737, 'fraction_numerical': 0.017073811400052534, 'mean_word_length': 4.132075471698113, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '421640', 'license': 'MIT', 'n_tokens_mistral': 1382, 'n_tokens_neox': 1278, 'n_words': 212, 'path': 'Real_Data/10x_PBMC_5k/run_methods/chromVAR/chromVAR_10xpbmc5k_motifs.ipynb', 'repo_name': 'DaneseAnna/scATAC-benchmarking', 'size': '3807'}