Average Nucleotide Identity (ANI) 计算

# install blast
mkdir blast
cd blast
wget http://mirrors.vbi.vt.edu/mirrors/ftp.ncbi.nih.gov/blast/executables/legacy/2.2.26/blast-2.2.26-universal-macosx.tar.gz
tar zxvf blast-2.2.26-universal-macosx.tar.gz -C /opt/biosoft
# blast path: /opt/biosoft/blast-2.2.26/bin/blastall
# formatdb path: /opt/biosoft/blast-2.2.26/bin/formatdb
# fastacmd path: /opt/biosoft/blast-2.2.26/bin/fastacmd
cd ..
# install MUMmer
mkdir MUMmer
cd MUMmer
wget https://sourceforge.net/projects/mummer/files/mummer/3.23/MUMmer3.23.tar.gz
tar zxvf MUMmer3.23.tar.gz -C /opt/biosoft
cd /opt/biosoft/MUMmer3.23
make check
make install
# nucmer path: /opt/biosoft/MUMmer3.23/nucmer
# install Jspecies
mkdir /opt/biosoft/jspecies
cd /opt/biosoft/jspecies
wget http://imedea.uib-csic.es/jspecies/jars/jspecies1.2.1.jar
java -jar -Xms1024m -Xmx1024m jspecies1.2.jar

 -X 空比对对应的下降值 [默认 150]
 -q 碱基错配罚分 [默认 -1]
 -F 是否过滤query序列的简单重复和低复杂度序列 [F]
 -e e-value [默认 1e-15]
 -a 运行线程数 [默认 2]
 Identity (%) >= 30
 Alignment (%) >= 70 #more than 30% overall sequence identity (recalculated to an identity along the entire sequence) over an alignable region of at least 70% of their length
 length 1020 #query序列被连续切割成1020bp的片段，然后将这些片段blastn比对到reference序列，比对的结果用于后续计算

       
 -b 低得分区大于此值时，停止比对 [默认 200]
 -c 一簇匹配的序列最少需要有大于此长度的匹配。此值越大，则小的匹配越少 [默认 65]
 -g 在一簇匹配的序列中，两个邻近的matches之间的gap的最大阈值。此值越大，则越容易找到更大匹配区 [默认 90]
 -l single match的最小长度 [默认 20]

# pip3 install
pip3 install pyani
echo "PATH=/home/linuxbrew/.linuxbrew/bin:$PATH" >> ~/.bashrc
source ~/.bashrc
# local install
git clone git@github.com:widdowquinn/pyani.git
cd pyani
python3 setup.py install

# example
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIm_output -m ANIm -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIb_output -m ANIb -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIblastall_output -m ANIblastall -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_TETRA_output -m TETRA -g

# required
  -i INDIRNAME, --indir INDIRNAME    输入文件所在目录
  -o OUTDIRNAME, --outdir OUTDIRNAME 输出结果目录
 
# optional
  -v, --verbose         输出详细结果
  -f, --force           强制文件的写入
  -s FRAGSIZE, --fragsize FRAGSIZE
                        ANIb比对时序列的片段大小 [默认 1020]
  -l LOGFILE, --logfile LOGFILE
                        输出日志文件
  --skip_nucmer         跳过运行NUCmer
  --skip_blastn         跳过运行BLASTN
  --noclobber           不重做已存在文件
  --nocompress          不压缩或删除比较结果
  -g, --graphics        生成ANI的热图
  --gformat GFORMAT     图形文件输出格式 [pdf|png|jpg|svg] [默认 pdf,png,eps]
  --gmethod {mpl,seaborn} 图形输出方法 [默认 mpl]
  --labels LABELS       序列标记文件的路径
  --classes CLASSES     序列分类文件的路径
  -m {ANIm,ANIb,ANIblastall,TETRA}, --method {ANIm,ANIb,ANIblastall,TETRA}
                        ANI比对方法 [默认 ANIm]
  --scheduler {multiprocessing,SGE}
                        Job scheduler [默认 multiprocessing]
  --workers WORKERS     运行线程数 [默认 全部线程]
  --SGEgroupsize SGEGROUPSIZE
                        SGE组中job数 [默认 10000]
  --SGEargs SGEARGS     Additional arguments for qsub
  --maxmatch            Override MUMmer to allow all NUCmer matches
  --nucmer_exe NUCMER_EXE
                        NUCmer 程序路径
  --filter_exe FILTER_EXE
                        delta-filter 程序路径
  --blastn_exe BLASTN_EXE
                        BLASTN+ 程序路径
  --makeblastdb_exe MAKEBLASTDB_EXE
                        BLAST+ makeblastdb 程序路径
  --blastall_exe BLASTALL_EXE
                        BLASTALL 程序路径
  --formatdb_exe FORMATDB_EXE
                        BLAST formatdb 程序路径
  --write_excel         输出Excel格式文件
  --rerender            Rerender graphics output without recalculation
  --subsample SUBSAMPLE
                        Subsample a percentage [0-1] or specific number (1-n)
                        of input sequences
  --seed SEED           Set random seed for reproducible subsampling.
  --jobprefix JOBPREFIX
                        Prefix for SGE jobs (default ANI).

Rscript heatmap_ANI.R $OUTPUT/ANIb_percentage_identity.tab

#!/usr/bin/env Rscript
library("ggplot2")
library("reshape")
 
userprefs <- commandArgs(trailingOnly = TRUE)
file_name <- userprefs[1]
 
# read file
data.ANI <- read.table(file_name)
 
# cluster and order matrix
row.order <- hclust(dist(data.ANI))$order # clustering
col.order <- hclust(dist(t(data.ANI)))$order
dat_new <- data.ANI[row.order, col.order] # re-order matrix accoring to clustering
 
# melt to dataframe
df_molten_dat <- melt(as.matrix(dat_new)) # reshape into dataframe
names(df_molten_dat)[c(1:3)] <- c("shewanella1", "shewanella2", "ANI")
 
# make plot
p1 <- ggplot(data = df_molten_dat,
aes(x = shewanella1, y = shewanella2, fill = ANI)) +
geom_raster() +
scale_fill_distiller(palette = "RdYlBu", limits=c(0.8, 1)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
ggtitle("Shewanella ANI heatmap")+
geom_text(aes(label = round(ANI, 2)))+
theme(axis.title=element_text(size=16), strip.text.x=element_text(size=16),
legend.title=element_text(size=15), legend.text=element_text(size=14),
axis.text = element_text(size=14), title=element_text(size=20),
strip.background=element_rect(fill=adjustcolor("lightgray",0.2))
)
 
# export plot
png("Heatmap_ANI_highres.png", width=15, height=15, units="in", pointsize=12, res=500)
p1
dev.off()