<jupyter_start><jupyter_text>Get sequences for datasetThis workflow to get sequences was inspired by [notebooks/functions from this repository](https://github.com/emdann/scATAC_prep).<jupyter_code>%load_ext rpy2.ipython
# Ignore R warning messages
#Note: this can be commented out to get more verbose R output
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)<jupyter_output><empty_output><jupyter_text>Load dataset<jupyter_code>ad = annd.read("data/Marrow_pooled.h5ad")
sc.pl.umap(ad, color='most_abundant_celltype')
gene_ids = ad.var["symbol"]
gene_ids[0:10]<jupyter_output><empty_output><jupyter_text>Use R to extract sequencesMost/all of these are available through R's package manager or Bioconductor. For human genome, use `library(EnsDb.Hsapiens.v86)` instead.<jupyter_code>%%R 
library(Matrix)
library(GenomicRanges)
library(ensembldb)
library(EnsDb.Mmusculus.v79) 
library(Biostrings)
library(tidyr)
library(AnnotationHub)
options(stringsAsFactors = F)
%%R 

genome <- ensembldb:::getGenomeTwoBitFile(EnsDb.Mmusculus.v79)
genome_genes <- ensembldb::genes(EnsDb.Mmusculus.v79)<jupyter_output><empty_output><jupyter_text>Since we are working with symbols, we set the gene names to symbols:<jupyter_code>%%R

names(genome_genes) <- genome_genes$symbol<jupyter_output><empty_output><jupyter_text>Get genes and filter out those that are not included:<jupyter_code>%%R -i gene_ids -o gene_ids

gene_ids <- as.vector(gene_ids)

good_seqname <- as.character(seqnames(genome_genes)) %in% c(as.character(1:19), "MT", "X", "Y") 
genome_genes <- genome_genes[good_seqname]

gene_ids <- gene_ids[gene_ids %in% genome_genes$symbol]<jupyter_output><empty_output><jupyter_text>Define sequence length relative to (canonical) TSS of genes and extract sequences<jupyter_code>%%R -i gene_ids -o promoter_seq

upstream = 500
downstream = 500
# This returns a GRanges object:
promoter_ranges <- GenomicRanges::promoters(genome_genes[gene_ids], 
                                            upstream = upstream, 
                                            downstream = downstream)

promoter_seq <- Biostrings::getSeq(genome, promoter_ranges)
promoter_seq <- tolower(as.character(promoter_seq))
ad = ad[:,ad.var["symbol"].isin(gene_ids)]
ad.var["promoter_seq"] = promoter_seq
sum(["N" in x for x in ad.var["promoter_seq"]])
ad.var
ad.write("data/Marrow_pooled.h5ad")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11622073578595318, 'fraction_numerical': 0.010033444816053512, 'mean_word_length': 5.997076023391813, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 25, 'https://': 1, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '261154', 'license': 'MIT', 'n_tokens_mistral': 898, 'n_tokens_neox': 825, 'n_words': 164, 'path': 'example_notebooks/02_get_sequences.ipynb', 'repo_name': 'jacobhepkema/scanem', 'size': '2392'}