如何制作私人网站,中国4a广告公司100强,公司部门解散 转岗不同意会怎么样,做网站网页需要什么技术文章目录 Comparem简介比较基因组统计基因组使用模式其他 安装使用基于基因组计算氨基酸一致性基于基因组蛋白计算氨基酸一致性 结果转变成矩阵参考 Comparem简介
CompareM 是一个支持进行大规模基因组比较分析的软件工具包。它提供跨基因组#xff08;如氨基酸一致性#x… 文章目录 Comparem简介比较基因组统计基因组使用模式其他 安装使用基于基因组计算氨基酸一致性基于基因组蛋白计算氨基酸一致性 结果转变成矩阵参考 Comparem简介
CompareM 是一个支持进行大规模基因组比较分析的软件工具包。它提供跨基因组如氨基酸一致性和单个基因组如密码子使用率的统计计算。 同时可以并行化以便能够扩展到数千个基因组。主要功能
比较基因组统计
基因组之间的平均氨基酸一致性AAI通过计算查询基因组与参考数据库之间的 AAI 进行物种分类
基因组使用模式
密码子使用氨基酸使用k 8 的 kmer 使用情况如四核苷酸终止密码子使用
其他
识别基因侧向转移LGT 的二核苷酸和密码子使用模式使用差异矩阵、分层聚类树和热图进行数据探索
安装
# 安装
(base)$ mamba install -c bioconda comparem使用
基于基因组计算氨基酸一致性
使用aai_wf 流程计算氨基酸一致性:
# 待比较基因组序列
(base) [yutaomyosin test]$ ls
GCA_001780165.1_genomic.fa HTR8_metabat2_bin.67.fa GCA_003235575.1_genomic.fa# 运行aai_wf
(base) [yutaomyosin test]$ time comparem aai_wf -c 30 -x fa . aaiwf_out aaiwf.log
real 0m23.745s
user 0m40.652s
sys 0m1.594s
# -c 线程
# -x 基因组后缀
# . 带比较基因组目录
# aaiwf_out 输出目录(base) [yutaomyosin test]$ cat aaiwf.log
[2022-01-24 11:25:04] INFO: CompareM v0.1.2
[2022-01-24 11:25:04] INFO: comparem aai_wf -c 30 -x fa . aaiwf_out
[2022-01-24 11:25:04] INFO: Identifying genes within genomes:Finished processing 3 of 3 (100.00%) genomes. # 共3个基因组总共3种AAI比较方式
[2022-01-24 11:25:19] INFO: Identified genes written to: aaiwf_out/genes # 对基因组预测蛋白
[2022-01-24 11:25:19] INFO: Appending genome identifiers to query genes.
[2022-01-24 11:25:19] INFO: Creating DIAMOND database (be patient!).
[2022-01-24 11:25:19] INFO: Performing self similarity sequence between genomes (be patient!). # 在基因组之间进行相似性比较
[2022-01-24 11:25:19] INFO: Sorting table with hits (be patient!).
[2022-01-24 11:25:20] INFO: Sequence similarity results written to: aaiwf_out/similarity
[2022-01-24 11:25:20] INFO: Calculating length of genes.
[2022-01-24 11:25:20] INFO: Indexing sorted hit table.
[2022-01-24 11:25:20] INFO: Calculating AAI between all 3 pairs of genomes: #计算3对基因组的AAIFinished processing 3 of 3 (100.00%) pairs.
[2022-01-24 11:25:20] INFO: Summarizing AAI results.
[2022-01-24 11:25:20] INFO: AAI between genomes written to: aaiwf_out/aai/aai_summary.tsv其中当前目录包含一组FASTA格式的基因组结果被写入一个名为aai_output的目录30个处理器应被用于计算结果。可以看到comparem最终是对所有基因组两两之间进行比较不考虑顺序所以是一个组合情况可以通过R的chose(n, k)得到最终的组合数例如3个基因组的最终需要进行3次比较16个基因组最终需要进行120次比较。还可以指定一些可选的参数。这包括用于定义基因组间互为最佳命中即同源物的序列相似性参数。默认情况下e值–evalue、序列同一性百分比–per_identity和比对长度百分比–per_aln_len参数被设置为1e-5、30%和70%。当指定要处理的基因组目录时CompareM只处理扩展名为fna的文件。这可以用-x(–file_ext)参数来改变。此外如果基因组已经由氨基酸蛋白质序列表示相对于基因组核苷酸序列这必须用–蛋白质标志来指定。否则将使用Prodigal从头识别基因。通过使用-cpus参数指定的多个处理器可以大大减少计算所有成对AAI值的时间。
基于基因组蛋白计算氨基酸一致性
对于基因组序列默认使用prodigal预测基因input_file参数表示要比较的基因组集合可以是i一个文本文件其中每一行表示一个基因组的位置或者ii一个包含所有要比较的基因组/蛋白氨基酸的目录。基因组/蛋白的序列必须是FASTA格式。output_dir表示所有输出文件目录。 对于蛋白序列直接使用faa氨基酸序列
nohup time comparem aai_wf --proteins -c 30 -x gz GTDBr214_479_B.anthracis_gene GTDBr214_479_B.anthracis_gene_aai aaiwf.log
# aai_wf AAI工作流程
# --proteins 指定输入文件是蛋白序列
# -c 线程
# -x 输入蛋白序列后缀名称结果
(base) [yutaomyosin Two_new_classes]$ ls aaiwf_out/
aai comparem.log genes similarity
(base) [yutaomyosin aaiwf_out]$ head aai/aai_summary.tsv
#Genome A Genes in A Genome B Genes in B # orthologous genes Mean AAI StdOrthologous fraction (OF)
HTR8_metabat2_bin.67 2502 GCA_001780165.1_genomic 3086 497 47.82 10.13 19.86
HTR8_metabat2_bin.67 2502 GCA_003235575.1_genomic 2464 430 47.97 10.59 17.45
GCA_001780165.1_genomic 3086 GCA_003235575.1_genomic 2464 965 52.99 11.96 39.16
(base) [yutaomyosin aaiwf_out]$成对的AAI统计数据在输出文件./output_dir/aai/aai_summary.tsv中提供。该文件由8列组成具体含义如下其中第6列即是AAI值。
1-第一个基因组的标识符
2-第一个基因组中的基因数
3-第二个基因组的标识符
4-第二个基因组中的基因数
5-两个基因组之间确定的直系同源基因的数量
6-直系同源基因的平均氨基酸一致性AAI。
7-直系同源基因的AAI的标准偏差
8-两个基因组之间的直系亲缘关系OF定义为直系亲缘关系的基因数除以其中一个基因组的最小基因数。转变成矩阵
上述长列表数据可以通过如下代码转换为矩阵随后可按照如下方式可视化20个基因组左右时适合。
if(!requireNamespace(pacman)){install.packages(pacman)}
pacman::p_load(igraph, corrplot, ggsci, ggplot2 )
# 用于将长表变成宽表
setwd(/Users/yut/Documents/12个海洋元基因组/MAGs_gene_functions/Two_novel_class/)
f Eisenbacteria_aai_summary(1).tsv # 读入compareM aai_wf输出长格式表
f Krumholzibacteriota_aai_summary(1).tsv
d - read.table(f, header F)[c(1, 3, 6)] # 取出第一个基因组第二个基因组及平均AAI
g graph.data.frame(d, directed FALSE) #使用igraph读入数据框
mat - get.adjacency(g, attr V6 # 显示的AAI值, spars F # 非稀疏没有的填充0)
mat[mat 0] - 100 # 自身AAI 100
#mat - round(mat / 100, 2) # 将百分比转化成小数便于作图
mat# 设置颜色映射
my.col - colorRampPalette(c(#4e9cb8 #最小值的颜色,#f2f1f1 # 中间颜色, #de6589 # 最大值颜色))(10) # 取10个连续色
# corrplot
pdf(Eisenbacteria_aai_summary.pdf, height 15, width 15)
pdf(Krumholzibacteriota_aai_summary.pdf, height 15, width 15)
p - corrplot(mat, method circle # cell shape, type lower # lower triangle, order hclust # 排序方式为层级聚类original, AOE, FPC, hclust, alphabet, hclust.method complete # 聚类方法 ward, ward.D, ward.D2, single, complete, average, mcquitty, median or centroid, col.lim c(0, 100) # 设置颜色mapping值的范围, is.corr F # 非相关系数的值,col my.col # 颜色, tl.col black # x/y坐标字体颜色, addCoef.col T # 显示值, title AAI (%))
dev.off()
# 不支持ggsave(plot p, Eisenbacteria_aai_summary.pdf, height 12, width 12, device pdf) 参考
comparem github
