当前位置： 首页 > wzjs >正文

做网站 江门域名ip查询

wzjs 2025/8/22 1:56:46
做网站 江门,域名ip查询,有什么简单的网站项目,设计公司品牌策划目录 BLAST比对 导入环境 构建SM对象 整合跨物种对象 保存结果 保存整合anndata及可视化 导入相关包 提取anndata对象 整合UMAP可视化 指定基因可视化 桑基图可视化 相似性热图可视化 之前的帖子详细介绍了如何从参考基因组和注释中提取蛋白序列，利用提…

目录

BLAST比对

导入环境

构建SM对象

整合跨物种对象

保存结果

保存整合anndata及可视化

导入相关包

提取anndata对象

整合UMAP可视化

指定基因可视化

桑基图可视化

相似性热图可视化

之前的帖子详细介绍了如何从参考基因组和注释中提取蛋白序列，利用提取的编码蛋白序列，以及经过预处理的单细胞/空间转录组的H5ad格式数据，我们就可以进行SAMap的跨物种数据整合

代码思路详见GitHub - atarashansky/SAMap: SAMap: Mapping single-cell RNA sequencing datasets from evolutionarily distant organisms.

BLAST比对

首先我们需要将提取的蛋白编码序列，也就是最长转录本的序列进行一一比对，有多少个物种就需要依次进行比对，

如有物种ABC，我们需要进行的比对则如下，每个物种需要用两个字母的缩写来代称，我这里对应的是aa/bb/cc，比对结果文件夹也会用上述缩写名称表示为aabb/aacc/bbcc，每个文件夹下是两个txt比对文件

bash map_genes.sh --tr1 A_protein_longest.fa --t1 prot --n1 aa --tr2 B_protein_longest.fa --t2 prot --n2 bb
bash map_genes.sh --tr1 A_protein_longest.fa --t1 prot --n1 aa --tr2 C_protein_longest.fa --t2 prot --n2 cc
bash map_genes.sh --tr1 B_protein_longest.fa --t1 prot --n1 bb --tr2 C_protein_longest.fa --t2 prot --n2 cc

官方提供的比对脚本如下所示（可以自己进行单独的BLAST，但是要保持同样的文件命名规则才可以继续进行整合）

#!/bin/bash#parse arguments
die() {printf '%s\n' "$1" >&2exit 1
}_help() {printf "Usage: map_genes.sh [--tr1] [--t1] [--n1] [--tr2] [--t2] [--n2]\n\t[--tr1]: path to transcriptome/proteome 1\n\t[--t1]: is 1 a transcriptome [nucl] or proteome [prot]\n\t[--n1]: two character identifier of 1\n\t[--tr2]: path to transcriptome/proteome 2\n\t[--t2]: is 2 a transcriptome [nucl] or proteome [prot]\n\t[--n2]: two character identifier of 2"
}if [ $# -eq 0 ]; then_helpexit 1
fiwhile :; docase $1 in-h|-\?|--help)_helpexit;;--tr1)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thentr1=$2shiftelsedie 'ERROR: "--tr1" requires a non-empty option argument.'fi;;--tr2)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thentr2=$2shiftelsedie 'ERROR: "--tr2" requires a non-empty option argument.'fi;;--n1)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thenn1=$2shiftelsedie 'ERROR: "--n1" requires a non-empty option argument.'fi;;--n2)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thenn2=$2shiftelsedie 'ERROR: "--n2" requires a non-empty option argument.'fi;;--t1)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thent1=$2shiftelsedie 'ERROR: "--t1" requires a non-empty option argument.'fi;;--t2)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thent2=$2shiftelsedie 'ERROR: "--t2" requires a non-empty option argument.'fi;;--threads)       # Takes an option argument; ensure it has been specified.if [ "$2" ]; thenn_threads=$2shiftfi;;--)              # End of all options.shiftbreak;;-?*)printf 'WARN: Unknown option (ignored): %s\n' "$1" >&2;;*)               # Default case: No more options, so break out of the loop.breakesacshift
doneif [ -z ${n_threads+x} ]; thenn_threads=8
fi#read -e -p 'Transcriptome 1: ' tr1
#read -e -p 'Transcriptome type: ' t1
#read -e -p 'Transcriptome ID: ' n1
#read -e -p 'Transcriptome 2: ' tr2
#read -e -p 'Transcriptome type: ' t2
#read -e -p 'Transcriptome ID: ' n2n1="${n1:0:2}"
n2="${n2:0:2}"mkdir -p "maps"
mkdir -p "maps/${n1}${n2}"if [ ! -f "${tr1}.nhr" ] && [ ! -f "${tr1}.phr" ]
then
makeblastdb -in $tr1 -dbtype $t1
fiif [ ! -f "${tr2}.nhr" ] && [ ! -f "${tr2}.phr" ]
then
makeblastdb -in $tr2 -dbtype $t2
fiif [[ "$t1" == "nucl" && "$t2" == "nucl" ]]
then
echo "Running tblastx in both directions"
tblastx -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
tblastx -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fiif [[ "$t1" == "nucl" && "$t2" == "prot" ]]
then
echo "Running blastx from 1 to 2 and tblastn from 2 to 1"
blastx -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
tblastn -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fiif [[ "$t1" == "prot" && "$t2" == "nucl" ]]
then
echo "Running tblastn from 1 to 2 and blastx from 2 to 1"
tblastn -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
blastx -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fiif [[ "$t1" == "prot" && "$t2" == "prot" ]]
then
echo "Running blastp in both directions"
blastp -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
blastp -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fi

导入环境

进行完物种比对之后就可以正式开始整合流程，首先导入所有需要的包

from samap.mapping import SAMAP
from samap.analysis import (get_mapping_scores, GenePairFinder,sankey_plot, chord_plot, CellTypeTriangles, ParalogSubstitutions, FunctionalEnrichment,convert_eggnog_to_homologs, GeneTriangles)
from samalg import SAM
import pandas as pd
import scanpy as sc
from anndata import AnnData

我们可以检查一下之前比对的结果文件是否有误

A=pd.read_csv('/data/work/BLAST/maps_out/maps/aabb/aa_to_bb.txt',sep='\t',index_col=0,header=None)
B=pd.read_csv('/data/work/BLAST/maps_out/maps/aabb/bb_to_aa.txt',sep='\t',index_col=0,header=None)
A.head()

比对文件的前两列应该是对应我们单细胞/空间转录组的基因名称

构建SM对象

接下来可以导入我们单细胞或者空间转录组的分析数据，构建SAMap分析对象sm

fn1 = 'A.h5ad'
fn2 = 'B.h5ad'
fn3 = 'C.h5ad'
filenames = {'aa':fn1,'bb':fn2,'cc':fn3}
#假设比对数据路径为/data/work/BLAST/maps_out/maps/
sm = SAMAP(filenames, f_maps = '/data/work/BLAST/maps_out/maps/', save_processed=True)

整合跨物种对象

紧接着进行整合过程，日志中会描述整合的细节

sm.run(pairwise=True)
samap = sm.samap # SAM object with three species stitched together
sm.scatter()#可视化整合对象

整合操作结束之后，这个sm对象里就包含了我们所需要的所有结果了

保存结果

我一般会将这里的sm对象先保存成pkl格式，下次读取使用即可

import pickle
with open('/data/work/sm_aabbcc.pkl', 'wb') as file:pickle.dump(sm, file)
print("对象已保存到本地文件。")#with open('/data/work/sm_aabbcc.pkl', 'rb') as file:
#    sm = pickle.load(file)
#print("对象已从本地文件加载。")
print(dir(sm))

这一步可以检查sm包含的所有数据结构，即

['ITER_DATA', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'gnnm', 'gnnm_refined', 'gns', 'gns_dict', 'gui', 'ids', 'pairwise', 'plot_expression_overlap', 'query_gene_pair', 'query_gene_pairs', 'refine_homology_graph', 'run', 'run_time', 'run_umap', 'samap', 'sams', 'scatter', 'smap']

保存整合anndata及可视化

为了后续的其他分析，我们一般会提取SAMap结果的anndata格式数据

导入相关包

import scanpy as sc
import matplotlib.pyplot as plt
import networkx as nx
import umap
from scipy.sparse import csr_matrix
import seaborn as sns
import holoviews as hv

提取anndata对象

samap = sm.samap
adata = samap.adata
adata

这里的adata就是我们需要的整合数据

整合UMAP可视化

sc.pl.umap(adata, color=['species'], title='UMAP projection colored by species', size=10,save='_species_projection.pdf')

指定基因可视化

这里要注意整合数据中不同物种的对应基因要加上之前定义的缩写的前缀，“XX_gene”这样

sc.pl.umap(adata, color='SB_BDNF', title='BDNF of sb',size=5,save='_SB_BDNF_projection.pdf')
sc.pl.umap(adata, color='NA_BDNF', title='BDNF of na',size=5,save='_NA_BDNF_projection.pdf')
sc.pl.umap(adata, color='SA_BDNF', title='BDNF of sa',size=5,save='_SA_BDNF_projection.pdf')

桑基图可视化

我们也可以用经典的sankeyplot表示不同物种分群之间的关系，这里的‘annotation’和‘layer’是转录组数据中的注释标签

keys = {'SB':'annotation','NA':'annotation','SA':'annotation','HM':'layer','MU':'layer'}
D,MappingTable = get_mapping_scores(sm,keys,n_top = 0)
#可视化比对表格检查
D.head()
MappingTable.head()
sankey_plot(MappingTable, align_thr=0.05, species_order=['SB', 'NA', 'SA', 'HM', 'MU'])

（但是这里自带的函数最多只能绘制三个物种的可视化，不知道为啥。。）

相似性热图可视化

跨物种比较最常用的就是热图可视化，这里我也利用seaborn进行了尝试

import seaborn as sns
import matplotlib.pyplot as plt# 假设 MappingTable 是一个 Pandas DataFrame，且已经包含了相似性分数
# 如果 MappingTable 是一个 NumPy 数组，也可以直接使用sns.set(style="white")# 绘制热图
plt.figure(figsize=(15, 12))  # 增大热图的尺寸
sns.heatmap(MappingTable,  # 数据源annot=True,  # 是否在热图上显示数值fmt=".2f",  # 数值的格式cmap="coolwarm",  # 色彩映射，coolwarm 是一个常见的选择linewidths=0.5,  # 单元格之间的分隔线宽度cbar_kws={"label": "Similarity Score"},  # 色标标签annot_kws={"size": 10},  # 设置注释的字体大小
)# 设置标题和标签
plt.title("Similarity Heatmap between Clusters", fontsize=18)  # 增大标题字体大小
plt.xlabel("Clusters", fontsize=16)  # 增大X轴标签字体大小
plt.ylabel("Clusters", fontsize=16)  # 增大Y轴标签字体大小# 设置坐标轴刻度标签的字体大小
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)# 保存可视化文件
plt.savefig("/data/work/figures/samap_heatmap.pdf", format="pdf", bbox_inches="tight")
plt.savefig("/data/work/figures/samap_heatmap.png", format="png", bbox_inches="tight")# 显示热图
plt.show()

http://www.dtcms.com/wzjs/436665.html

相关文章：

  • 企业网站建设方案书 范本百度指数人群画像哪里查询
  • 网站怎么做图片动态图片不显示了西安百度推广优化托管
  • pageadmin的最新版本seo接单一个月能赚多少钱
  • 佛山企业网站建设电话他达拉非的副作用和危害
  • 泰州做企业网站发布信息的免费平台
  • 网站策划书ppt网站建设及网站推广
  • 做网站开发学什么语言好seo单页快速排名
  • 深圳龙华建设发展公司石景山区百科seo
  • 网站建设江门免费推广网站入口
  • 企业网站推广工具网站免费下载安装
  • 上海网站建设哪家好网络服务平台
  • 竞拍网站模板关键词推广营销
  • 互动网站案例在线识别图片找原图
  • 网站运营分析竞争对手今日要闻新闻
  • 有哪些免费的服务器沈阳百度seo排名优化软件
  • 谷歌俄语网站网站seo外链平台
  • 网站开发h5技术sem是什么意思啊
  • 织梦仿商城网站精品成品网站源码
  • wordpress文章摘要插件seo网站技术培训
  • 湘潭做网站选择磐石网络浙江seo
  • b2b2c网站开发广东省广州市白云区
  • 班级网站源代码下载seo怎么赚钱
  • 厦门网站建设外贸方象科技专注于什么领域
  • 手机下载视频网站模板网络营销推广的方法有哪些
  • 蚌埠网站制作哪里有今日时政新闻
  • 网站建设内容的重点优化培训内容
  • 网站宣传单模板做整站优化
  • 网站开网站开发设计公司广州最新发布最新
  • 企业网站制作服务可以打广告的平台
  • 境外社交网站上做推广国际军事最新消息今天