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

获取脚本

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)