<jupyter_start><jupyter_text>Run BLAST<jupyter_code>file1 = 'path/to/transcriptome1.fasta'
type1 = 'nucl' #or 'prot' if file1 is a proteome
id1 = 'hu' #2-character ID (e.g. 'hu' for human)

file2 = 'path/to/transcriptome2.fasta'
type2 = 'nucl'
id2 = 'mo' #2-character ID (e.g. 'mo' for mouse)
!bash map_genes.sh --tr1 {file1} --t1 {type1} --n1 {id1} --tr2 {file2} --t2 {type2} --n2 {id2}<jupyter_output><empty_output><jupyter_text>Run SAMap<jupyter_code>from samap.mapping import SAMAP
from samap.analysis import get_mapping_scores, GenePairFinder
from samalg import SAM<jupyter_output><empty_output><jupyter_text>SAMap accepts file paths to unprocessed, raw `.h5ad` files. Alternatively, if you already have a processed `SAM` object, you can load them in directly.Prior to running SAMap, you should have run `map_genes.sh`, which expects 2-character identifiers describing each species. For example, human and mouse might get `hu` and `mo` identifiers, respectively. `map_genes.sh` generates a `maps/` directory with the transcriptome mapping BLAST results deposited. The input species IDs and path to the `maps/` directory should be input into SAMap.<jupyter_code>id1 = 'hu'
id2 = 'mo'
# passing in file names (SAMap will process the data with SAM and save the resulting objects to two `.h5ad` files.)
fn1 = '/path/to/file1.h5ad' #processed data will be automatically saved to `/path/to/file/file1_pr.h5ad`
fn2 = '/path/to/file2.h5ad' #processed data will be automatically saved to `/path/to/file/file2_pr.h5ad`
# runs SAMAP 
sm = SAMAP(fn1,fn2,id1,id2,f_maps = 'maps/')
samap = sm.run()

"""
# passing in already-processed SAM objects
sam1=SAM()
sam2=SAM()
sam1.load_data('/path/to/file1_pr.h5ad')
sam2.load_data('/path/to/file2_pr.h5ad')

sm = SAMAP(sam1,sam2,id1,id2,f_maps = 'maps/')
samap = sm.run()
"""<jupyter_output><empty_output><jupyter_text>To calculate alignment scores between cell types, we can use `get_mapping_scores`. This function will use the combined SAM object produced by SAMap to calculate alignment scores between cell types in the provided cell type annotation columns of `sam.adata.obs`. If no cell type annotations exist, the leiden clusters generated by SAM can be used (`k1=k2='leiden_clusters'`).The resulting tables show the highest-scoring alignment scores for each cell type in organism 1 (`D1`) and organism 2 (`D2`), respectively. The `n_top` parameter can be used to identify strongly-mapping subpopulations between clusters. If `n_top=0`, the alignment score will be averaged over all cells between two clusters. Otherwise, the alignment score will be averaged over the top `n_top` cells.For example, let's say there are 100 / 1000 cells in Cluster A that maps to a Cluster B in the other species. On average, the mapping score between cluster A and B will be small, because there are only a few cells in cluster A that map to cluster B. If we set `n_top` to 100, then the alignment score between A and B will be averaged over the 100 highest-scoring cells. As a result, the mapping score between A and B will be large.By default, `n_top=0`.<jupyter_code>k1 = 'cell_types' #cell types annotation key in `sam1.adata.obs`
k2 = 'cell_types' #cell types annotation key in `sam2.adata.obs`
D1,D2,MappingTable = get_mapping_scores(sm,k1,k2, n_top = 0)<jupyter_output><empty_output><jupyter_text>SAMap provides a class to find gene pairs enriched in different cell type pairs. The method entails finding gene pairs that contribute positively to the cross-species correlation between cell types and are differentially expressed in their respective mapped cell types.<jupyter_code>gpf = GenePairFinder(sm,k1=k1,k2=k2)<jupyter_output><empty_output><jupyter_text>`gpf.find_genes` can now be used to find gene pairs enriched in a cell type mapping.<jupyter_code>n1 = 'cell_type1' #cell type ID from organism 1 (must be present in `sam1.adata.obs[k1]`)
n2 = 'cell_type2' #cell type ID from organism 2 (must be present in `sam2.adata.obs[k1]`)
Gp,G1,G2 = gpf.find_genes(n1,n2)
#Gp are the gene pairs, G1 are the genes from organism 1, G2 are the genes from organism 2<jupyter_output><empty_output><jupyter_text>To get a table of enriched gene pairs from all cell type mappings that have an alignment score above some threshold (`thr`), you can use:<jupyter_code>gene_pairs = gpf.find_all(thr=0.1)<jupyter_output><empty_output><jupyter_text>Saving/Loading SAMap<jupyter_code>from samap.utils import save_samap, load_samap
# save to a .pkl file
save_samap(sm,'path/to/file') #including the file ending (.pkl) is optional
# load
sm = load_samap('path/to/file') #including the file ending (.pkl) is optional<jupyter_output><empty_output><jupyter_text>Visualizing SAMap resultsDisplaying a heatmap:<jupyter_code>sm.display_heatmap(key1='leiden_clusters',key2='leiden_clusters')<jupyter_output><empty_output><jupyter_text>Launching an interactive GUI (requires SAMGUI, see the README instructions in the [SAM](https://github.com/atarashansky/self-assembling-manifold) github repo):<jupyter_code>sm.gui()<jupyter_output><empty_output><jupyter_text>To launch a similar interactive GUI with both species overlaid in one tab, run:<jupyter_code>sm.samap.gui()<jupyter_output><empty_output><jupyter_text>To create a sankey plot, use `samap.analysis.sankey_plot`. This requires `holoviews` to be installed (`pip install holoviews`). It should already be installed if you're using the Docker image.<jupyter_code>k1 = 'cell_types' #cell types annotation key in `sam1.adata.obs`
k2 = 'cell_types' #cell types annotation key in `sam2.adata.obs`
D1,D2,MappingTable = get_mapping_scores(sm,k1,k2, n_top = 0)
sankey_plot(MappingTable)<jupyter_output><empty_output><jupyter_text>Functional Enrichment AnalysisWe can use the `samap.mapping.FunctionalEnrichment` class to perform functional enrichment analysis of gene pairs enriched in each cell type. The functional annotations can be GO terms, KOG terms, etc.We can get functional annotations for our data by mapping the transcriptome to the Eggnog database. Below, we clone the eggnog repository from github, download the eggnog database, and run the mapping.<jupyter_code># run the below line if you are not running the SAMap Docker image.
#!pip install git+https://github.com/eggnogdb/eggnog-mapper.git
# clone the repo (this only needs to be done once)
!git clone https://github.com/eggnogdb/eggnog-mapper.git # put in dockerfile
# download the database (this only needs to be done once)
!eggnog-mapper/download_eggnog_data.py -y<jupyter_output><empty_output><jupyter_text>`transcriptome1.fa` should be the path to your first transcriptome fasta file.\`transcriptome2.fa` should be the path to your second transcriptome fasta file, and so on.\`output_name` is the name of the output file.<jupyter_code># run EGGNOG
!eggnog-mapper/emapper.py -m diamond -i transcriptome1.fa --cpu 0 --itype CDS -o output_name1
!eggnog-mapper/emapper.py -m diamond -i transcriptome2.fa --cpu 0 --itype CDS -o output_name2<jupyter_output><empty_output><jupyter_text>The above command should eventually output an `output_name1.emapper.annotations` table.\This table has many columns, including a collection of GO terms for each transcript and KOG functional categories. In the below example, I perform functional enrichment analysis of the KOG categories (`best_OG_cat` column in the EGGNOG tables).`FunctionalEnrichment.calculate_enrichment()` returns three tables. The rows are cell types, and the columns are functional categories. `enrichment_scores` provides the `-log10` statistical significance p-values. High values indicate high enrichment of gene pairs with that functional category. `num_enriched_genes` is the number of genes with a particular functional annotation enriched in a particular cell type mapping. `enriched_genes` provides a `;`-separated list of the genes.<jupyter_code>from samap.analysis import FunctionalEnrichment
# load eggnog tables
A = pd.read_csv('output_name1.emapper.annotations',sep='\t',index_col=0)
B = pd.read_csv('output_name2.emapper.annotations',sep='\t',index_col=0)
# load SAMAP
sm = load_samap('samap_object.pkl')
# do functional enrichment
fe = FunctionalEnrichment(sm,[A,B],'best_OG_cat',['leiden_clusters','leiden_clusters'])
enrichment_scores,num_enriched_genes,enriched_genes = fe.calculate_enrichment(verbose=True)<jupyter_output><empty_output><jupyter_text>See below for a description of all inputs to `FunctionalEnrichment`:```FunctionalEnrichment(sms, DFS, col_key, keys, delimiter = '', align_thr = 0.1, limit_reference = False, n_top = 0)Parameters----------sms - list or tuple of SAMAP objectsDFS - list or tuple of pandas.DataFrame functional annotations (one per species present in the input SAMAP objects)col_key - str    The column name with functional annotations in the annotation DataFrames.keys - list or tuple of column keys from `.adata.obs` DataFrames (one per species present in the input SAMAP objects)    Cell type mappings will be computed between these annotation vectors.delimiter - str, optional, default ''    Some transcripts may have multiple functional annotations (e.g. GO terms or KOG terms) separated by    a delimiter. For KOG terms, this is typically no delimiter (''). For GO terms, this is usually a comma    (',').align_thr - float, optional, default 0.1    The alignment score below which to filter out cell type mappingslimit_reference - bool, optional, default False    If True, limits the background set of genes to include only those that are enriched in any cell type mappings    If False, the background set of genes will include all genes present in the input dataframes.n_top: int, optional, default 0    If `n_top` is 0, average the alignment scores for all cells in a pair of clusters.    Otherwise, average the alignment scores of the top `n_top` cells in a pair of clusters.    Set this to non-zero if you suspect there to be subpopulations of your cell types mapping    to distinct cell types in the other species.``` Enrichment scores can be plotted using `fe.plot_enrichment()`:```FunctionalEnrichment.plot_enrichment(self,cell_types = [], pval_thr=2.0,msize = 50)Create a plot summarizing the functional enrichment analysis.Parameters----------cell_types - list, default []    A list of cell types for which enrichment scores will be plotted. If empty (default),    all cell types will be plotted.pval_thr - float, default 2.0    -log10 p-values < 2.0 will be filtered from the plot.msize - float, default 50    The marker size in pixels for the dot plot.Returns-------fig - matplotlib.pyplot.Figureax - matplotlib.pyplot.AxesSubplot```<jupyter_code>fig,ax = fe.plot_enrichment(cell_types = [], pval_thr=2.0,msize = 50)<jupyter_output><empty_output><jupyter_text>Identifying Homolog/Paralog SubstitutionsOne of the key features of SAMap is the ability to automatically identify all instances where paralogs have much more similar expression patterns than their respective orthologs. To do this, all we need are the SAMAP objects and a list of ortholog pairs and paralog pairs. The paralogs can be either cross-species, within-species, or a mix of both. Within-species paralogs are automatically converted to cross-species paralogs by using the orthologs.***IMPORTANT***: For all genes in `ortholog_pairs` and `paralog_pairs`, this function expects the genes tobe prepended with their corresponding species IDs (i.e. `sm.id1` or `sm.id2`). This is crucial because some pairs of transcripts between transcriptomes may have the same names, so we need to make them unique by prepending with the speciesIDs. If `paralog_pairs` is left to its default value (`None`), then all genes in the homology graph that aren't orthologs are treated as if they were paralogs. In this case, the analysis would give you "homolog" substitutions rather than paralog substitutions.The documentation for `ParalogSubstitutions` is:```ParalogSubstitutions(sm, ortholog_pairs, paralog_pairs=None, psub_thr = 0.3):Identify paralog substitutions. For all genes in `ortholog_pairs` and `paralog_pairs`, this function expects the genes tobe prepended with their corresponding species IDs (i.e. `sm.id1` or `sm.id2`).Parameters----------sm - SAMAP objectortholog_pairs - n x 2 numpy array of ortholog pairsparalog_pairs - n x 2 numpy array of paralog pairs, optional, default None    If None, assumes every pair in the homology graph that is not an ortholog is a paralog.    Note that this would essentially result in the more generic 'homolog substitutions' rather    than paralog substitutions.    The paralogs can be either cross-species, within-species, or a mix of both. psub_thr - float, optional, default 0.3    Threshold for correlation difference between paralog pairs and ortholog pairs.    Paralog pairs that do not have greater than `psub_thr` correlation than their     corresponding ortholog pairs are filtered out.Returns-------RES - pandas.DataFrame    A table of paralog substitutions.```<jupyter_code>TABLE = ParalogSubstitutions(sm, ortholog_pairs, paralog_pairs=None, psub_thr = 0.3)
TABLE.head()<jupyter_output><empty_output><jupyter_text>Identifying orthologs  using EGGNOG.We can identify orthologs at different ancestral levels using the mapping tables output by EGGNOG.In the `eggNOG_OGs` column of the output mapping table are a comma-separated list of orthology groups at different ancestral levels that each transcript maps to. The levels are typically Eukaryota, Opisthokonta, Metazoa, Bilateria, Chordata, Vertebrata, and so on. We can call genes that map to the most recent common ancestor level (e.g. Vertebrata) as orthologs, and genes that map to a more ancestral level (e.g. Metazoa) as paralogs.For example, below we identify genes as orthologs if they map to the same Vertebrata orthology group and paralogs if they map to the same Metazoa orthology group.The documentation for `convert_eggnog_to_homologs` is:```   convert_eggnog_to_homologs(sm, A, B, og_key = 'eggNOG_OGs', taxon=2759)Gets an n x 2 array of homologs at some taxonomic level based on Eggnog results.    Parameters----------smp: SAMAP objectA: pandas.DataFrame, Eggnog output tableB: pandas.DataFrame, Eggnog output tableog_key: str, optional, default 'eggNOG_OGs'    The column name of the orthology group mapping results in the Eggnog output table.    taxon: int, optional, default 2759    Taxonomic ID corresponding to the level at which genes with overlapping orthology groups    will be considered homologs. Defaults to the Eukaryotic level.Returns-------homolog_pairs: n x 2 numpy array of homolog pairs.```<jupyter_code>ortholog_pairs = convert_eggnog_to_homologs(sm, A, B, og_key = 'eggNOG_OGs', taxon=7742)
paralog_pairs = convert_eggnog_to_homologs(sm, A, B, og_key = 'eggNOG_OGs', taxon=33208)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08298867719845414, 'fraction_numerical': 0.012746626889958642, 'mean_word_length': 5.36873920552677, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 73, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '398470', 'license': 'MIT', 'n_tokens_mistral': 4675, 'n_tokens_neox': 4352, 'n_words': 1747, 'path': 'SAMap_quickstart.ipynb', 'repo_name': 'Bassi-git/SAMap', 'size': '14749'}