刚发的MIKE能直接从测序数据快速建树,第一次见,感觉挺神奇的,就来试试。
参考:
Wang, F., Wang, Y., Zeng, X., Zhang, S., Yu, J., Li, D., Zhang, X., & Schwartz, R. (2024). MIKE: an ultrafast, assembly- and alignment-free approach for phylogenetic tree construction. Bioinformatics. https://doi.org/10.1093/bioinformatics/btae154
个人理解就是相当于把整个基因组都切成kmer,然后看kmer种类的统计差异,可能类似于SNPcall之后用位点建树。(感觉理解上还是有出入,还望有仔细研读原文的大佬能不吝指摘)
24.4.8 加笔
差不多看明白了,统计的是kmer的种类。统计Jaccard系数,算的是两个物种kmer种类的交集与并集的比值。
举个例子:原文推荐用k21去切,那么得到的就是21个碱基长度的片段,这个片段的种类就有4的21次方个,但是一个物种理论上不太可能有所有的这些kmer,而且基因组越小的物种越不可能。如果两个物种的亲缘关系很近,那么它们共有的kmer种类就更多(在理想情况下)。
那么基因组大小变化幅度较大的内类群是否适合使用这种方法就是一个问题。
使用的kmer越长,物种自有的kmer种类占比就越大;kmer越短,共有的占比就越大。用理论物理常用的极限法设想一下的话,kmer长度为基因组长度或1个碱基的话,都没有分辨率,那么对于不同研究类群而言21恒定为最佳的值吗?不同的kmer取值是仅仅影响分辨率,还是说会影响拓扑?
先说体验结论:
很有趣,很快;
系统树的主要框架跟核基因树基本一致;
还有就是如原文提到的那样- more suitable for constructing phylogenetic trees at the genus level
流程搬运:
Kmer文件的准备
先下载KMC,KMC用于计数kmer
#先下载最新版本的kmc并上传至服务器
#kmc下载:https://github.com/refresh-bio/KMC/releases
#解压
tar -zxvf KMC3.2.4.linux.x64.tar.gz
#配置环境变量,将解压后出现的bin文件夹添加到.bashrc里,这里不详细赘述
#查看说明
kmc
#主要用法
#kmc [options] <input_file_name> <output_file_name> <working_directory>
#kmc [options] <@input_file_names> <output_file_name> <working_directory> 示例
#单端
kmc -k21 -t10 -fq 1_1.clean.fq k21_1 .
#双端
kmc -k21 -t10 -fq @1.lst k21_1 . lst文件内容如下
/1_1.clean.fq
/1_2.clean.fq
如果fasta或fastq文件过大,可能会报错:
Error: Wrong input file (kmc_core/fastq_reader.cpp: 850)
可以尝试临时调高系统的最大打开文件数
ulimit -n 10000
#查看
ulimit -a 如果依旧不行的话,可以尝试先转化成bam文件,因为我跑HybPiper的时候有生成bam文件,就直接拿来用了(注:这里是为了示范,所以将就着用了一下,但是hybpiper生成的bam应该是根据reference生成的,所以结果估计也是对应挑取单拷贝的reference,并不能代表整个基因组数据情况。真正使用的时候可以用bowtie2转,这里不赘述)。
mkdir kmc_tmp
#中间过程临时文件放置在临时文件夹中
kmc -k21 -t10 -fbam hybpiper/1/1.bam k21_1 kmc_tmp 生成文件 k21_1.kmc_pre k21_1.kmc_suf
接下来用 kmc_tools将二进制文件转换存储到txt文本中
kmc_tools transform k21_1 sort . dump -s k21_1.txt 至此就已经准备好一个样的kmer文件了,所有样本的kmer文件都准备好后就可以开始下一步
建树
安装Mike
git clone https://github.com/Argonum-Clever2/mike.git
cd src
make
# 需要提前安装R
Rscript install.r
#同理,将src文件夹加入环境变量
准备filelist,包含上一步生成kmer文件的绝对路径
绝对路径/k21_1.txt
.../k21_2.txt
.../k21_3.txt
...
生成sketched文件
./mike sketch -t 10 -l filelist -d . 生成的文件以.minhash.jac结尾
将sketched文件保存在一个sklist里面,参考上面的
统计Jaccard系数
mike compute -l sklist -L sklist -d . 生成 jaccard.txt
计算遗传距离
mike dist -l sklist -L sklist -d . 生成 dist.txt
生成树文件
#draw.r要加绝对路径才能用
Rscript ABSOLUTE_PATH/draw.r -f dist.txt -o dist.nwk
如果报错也可以在本地的Rstudio里运行
#install.packages("ape")
library(ape)
tree <- read.csv("absolute_path/dist.txt", sep='\t', header = TRUE, row.names = 1)
treedist <- as.dist(tree)
# bionj
tree <- bionj(treedist)
# nj
tree <- nj(treedist)
# output
write.tree(tree, "tree.nwk")
plot(tree)