Github:
https://github.com/khigashi1987/Python_PCoA 3年前 python3
https://github.com/mwguthrie/python_PCoA 8年前 python
安装依赖
# 检测依赖
import numpy as np #ok
import itertools #ok
import math #ok
import argparse #ok
import pandas as pd #ok
import distance #
import matplotlib.pyplot as plt #
import scipy.stats #
from scipy.spatial.distance import squareform #
from sklearn import manifold #
# 缺少依赖
No module named 'sklearn'
No module named 'scipy'
No module named 'distance'
No module named 'matplotlib'
# 安装依赖
conda install scikit-learn # 同时安装了，scipy
conda install distance
conda install matplotlib
# JSD distance
conda install informationTheory
获取脚本
wget https://github.com/khigashi1987/Python_PCoA/archive/refs/heads/master.zip
tar -zxvf Python_PCoA-master.zip
cd Python_PCoA-master
python3 pcoa.py --help
########################################################
usage: pcoa.py [-h] [-f DATA_FILE] [-d {Jaccard,BrayCurtis,JSD}] [-b]
[-n N_ARROWS] [-g GROUP_FILE]
optional arguments:
-h, --help show this help message and exit
-f DATA_FILE, --file DATA_FILE
tab-separated text file. rows are variables, columns
are samples.
-d {Jaccard,BrayCurtis,JSD}, --distance_metric {Jaccard,BrayCurtis,JSD}
choose distance metric used for PCoA.
-b, --biplot output biplot (with calculating factor loadings).
-n N_ARROWS, --number_of_arrows N_ARROWS
how many top-contributing arrows should be drawed.
-g GROUP_FILE, --grouping_file GROUP_FILE
plot samples by same colors and markers when they
belong to the same group. Please indicate Tab-
separated 'Samples vs. Group file' ( first columns are
sample names, second columns are group names ).
计算bray-pcoa
source /public/home/zzumgg03/huty/softwares/miniconda3/etc/profile.d/conda.sh
conda activate py37
# success
python3 ./Python_PCoA-master/pcoa.py -f tpm.txt \
--distance_metric BrayCurtis
参数：
-f 输入，行-变量，列-样本
-d 距离算法{Jaccard,BrayCurtis,JSD}
结果，只有一张图，
结果只给了一张图，但我需要距离矩阵和坐标信息。
改写脚本，只计算bray-curtis距离
#!/usr/bin/env python3
# vim: set fileencoding=utf-8 :
import numpy as np
from scipy.spatial.distance import squareform
import itertools
import argparse
import pandas as pd
import distance
from sklearn import manifold
import scipy.stats
import matplotlib.pyplot as plt
# 函数
datafile = "./tpm.test" # 行-基因，列-样本，800,000,000 genes x 150 samples
data = pd.read_table(datafile, index_col=0)
matrix = data.values
from scipy.spatial.distance import braycurtis
X = matrix.T # 旋转90度
X = np.array(X) # 转矩阵
n_samples = X.shape[0] # shape 几行 [0] 几列 [1]
n_distance = int(n_samples * (n_samples - 1) / 2) # 比较的次数
d_array = np.zeros((n_distance))
for i, (idx1, idx2) in enumerate(itertools.combinations(range(n_samples),2)):
d_array[i] = braycurtis(X[idx1], X[idx2])
tmp = squareform(d_array)
df = pd.DataFrame(data = tmp[0:,0:],
columns = data.columns,
index = data.columns)
df.to_csv('./out_distance.txt', sep='\t')
得到距离矩阵
R中进行PCOA
# 读取距离矩阵
dist = read.table("out_distance.txt", header=T, sep="\t", row.names=1)
# 分组
meta = read.table("../metadata/merge_information_abundance.txt", sep="\t", header=T, row.names=1)
meta = meta[rownames(dist),]
group = data.frame(Layer = meta$Layer)
rownames(group) = rownames(meta)
# pcoa 总体
library(vegan)
library(ggplot2)
library(ape)
library("ggrepel")
options(scipen = 3)
pcoa = pcoa(dist)
pcoa_point = data.frame(pcoa$vectors[, c(1, 2)],
Layer = group$Layer)
## 关于分组需要自己定义
xylab = paste(c("PCoA1: ", "PCoA2: "), round(pcoa$values$Relative_eig[1:2]*100, 2), "%", sep = "")
p =
ggplot(pcoa_point, aes(x=Axis.1, y=Axis.2, color=Layer)) +
geom_point(size=3) +
geom_text_repel(aes(label = rownames(pcoa_point)),
size = 5, color = "black") +
stat_ellipse(level = 0.95, show.legend = F, size = 1) +
theme_classic() +
labs(x=xylab[1], y=xylab[2], title="Bray-Curtis PCoA") +
theme(legend.text=element_text(size=15)) +
theme(legend.title=element_text(face='bold', size=20)) +
theme(axis.title = element_text(size = 20)) +
theme(axis.text = element_text(size = 20),
axis.line = element_line(size = 1),
axis.ticks = element_line(size = 1))
ggsave(p, file="geneset_bray_pcoa.pdf", width=9)
