这是那篇原始的文章：Arumugam, M., Raes, J., et al. (2011) Enterotypes of the human gut microbiome, Nature,doi://10.1038/nature09944 在谷歌上一搜，作者竟然做了个分析肠型的教程在这，学习一下：http://enterotyping.embl.de/enterotypes.html 这是2018年大佬们的共识文章：这是国人翻译的这篇文章，http://blog.sciencenet.cn/blog-3334560-1096828.html 当然，如果你只需要获得自己的结果或者自己课题的结果，不需要跑代码的，有最新的网页版分型，更好用，网址也放在这，同样也是上面翻译的那篇文章里提到的网址：http://enterotypes.org/ 只需要把菌属的含量比例文件上就能很快得到结果。

下面我就边学习边做来尝试着来个分析，并把代码放在这里备忘。其实作者已经整理好了代码，我学习一下，争取实现对手上的数据进行分析。

首先下载测试数据，

wget http://enterotyping.embl.de/MetaHIT_SangerSamples.genus.txt
wget http://enterotyping.embl.de/enterotypes_tutorial.sanger.R

跑跑示例数据，排排错

我表示对R语言还只是一知半解的状态，所以，先跑下，然后能用上自己的数据， 当个工具用就暂知足啦。我是黑苹果10.11的系统，运行这个软件提示少了Xquartz，于是装了个，windows和linux应该不需要。原代码中还提示『没有”s.class”这个函数』，百度了一下发现有个老兄的新浪博客说了是这个包，于是加了句library(ade4)就ok了。 Xquartz的下载地址Mac 10.6+：https://dl.bintray.com/xquartz/downloads/XQuartz-2.7.11.dmg

#Uncomment next two lines if R packages are already installed
#install.packages("cluster")
#install.packages("clusterSim")
library(cluster)
library(clusterSim)
#BiocManager::install("genefilter")
library(ade4)#Download the example data and set the working directory
#setwd('<path_to_working_directory>')
data=read.table("../MetaHIT_SangerSamples.genus.txt", header=T, row.names=1, dec=".", sep="\t")
data=data[-1,]dist.JSD <- function(inMatrix, pseudocount=0.000001, ...) {
 KLD <- function(x,y) sum(x *log(x/y))
 JSD<- function(x,y) sqrt(0.5 * KLD(x, (x+y)/2) + 0.5 * KLD(y, (x+y)/2))
 matrixColSize <- length(colnames(inMatrix))
 matrixRowSize <- length(rownames(inMatrix))
 colnames <- colnames(inMatrix)
 resultsMatrix <- matrix(0, matrixColSize, matrixColSize) inMatrix = apply(inMatrix,1:2,function(x) ifelse (x==0,pseudocount,x)) for(i in 1:matrixColSize) {
   for(j in 1:matrixColSize) {
     resultsMatrix[i,j]=JSD(as.vector(inMatrix[,i]),
                            as.vector(inMatrix[,j]))
   }
 }
 colnames -> colnames(resultsMatrix) -> rownames(resultsMatrix)
 as.dist(resultsMatrix)->resultsMatrix
 attr(resultsMatrix, "method") <- "dist"
 return(resultsMatrix)
}data.dist=dist.JSD(data)pam.clustering=function(x,k) { # x is a distance matrix and k the number of clusters
 require(cluster)
 cluster = as.vector(pam(as.dist(x), k, diss=TRUE)$clustering)
 return(cluster)
}data.cluster=pam.clustering(data.dist, k=3)require(clusterSim)
nclusters = index.G1(t(data), data.cluster, d = data.dist, centrotypes = "medoids")nclusters=NULLfor (k in 1:20) {
 if (k==1) {
   nclusters[k]=NA
 } else {
   data.cluster_temp=pam.clustering(data.dist, k)
   nclusters[k]=index.G1(t(data),data.cluster_temp,  d = data.dist,
                         centrotypes = "medoids")
 }
}plot(nclusters, type="h", xlab="k clusters", ylab="CH index",main="Optimal number of clusters")obs.silhouette=mean(silhouette(data.cluster, data.dist)[,3])
cat(obs.silhouette) #0.1899451#data=noise.removal(data, percent=0.01)## plot 1
obs.pca=dudi.pca(data.frame(t(data)), scannf=F, nf=10)
obs.bet=bca(obs.pca, fac=as.factor(data.cluster), scannf=F, nf=k-1)
dev.new()
s.class(obs.bet$ls, fac=as.factor(data.cluster), grid=F,sub="Between-class analysis", col=c(3,2,4))#plot 2
obs.pcoa=dudi.pco(data.dist, scannf=F, nf=3)
dev.new()
s.class(obs.pcoa$li, fac=as.factor(data.cluster), grid=F,sub="Principal coordiante analysis", col=c(3,2,4))

上图，稍微调整下

, col=c(3,2,4)这个代码是给三个聚类上不同的颜色，还没搞清楚怎么给画的圈上色来实现理好看的效果，相信对于熟悉R语言的同学是小菜一碟。, cell=0, cstar=0是不显示圈和边线，只显示散点。 不加这两个参数，只用上面的代码，图如下：

加上两个参数的图片，就和教程里的最后面的两张图一样：

rsync -P -avz -e “ssh -p5566 -i /drives/C/Users/Admin/Desktop/cityu_jhli_1.id” /drives/d/Kraken_12.tar.gz cityu_jhli_1@172.16.22.11:/BIGDATA1/cityu_jhli_1/mhyleung/database/findfungi

[Admin.DESKTOP-7JT504C] ➤ scp -P 5566 -i /drives/C/Users/Admin/Desktop/cityu_jhli_1.id -r /drives/d/Kraken_12.tar.gz cityu_jhli_1@172.16.22.11:/BIGDATA1/cityu_jhli_1/mhyleung/database/findfungi

ssh -p5566 -i /drives/C/Users/Admin/Desktop/cityu_jhli_1.id   cityu_jhli_1@172.16.22.11

Docker运行

默认创建mkcp + 随机一种伪装头配置文件：

docker run -d --name v2ray --privileged --restart always --network host jrohy/v2ray

自定义v2ray配置文件:

docker run -d --name v2ray --privileged -v /path/config.json:/etc/v2ray/config.json --restart always --network host jrohy/v2ray

查看v2ray配置:

docker exec v2ray bash -c "v2ray info"

warning: 如果用centos，需要先关闭防火墙

systemctl stop firewalld.service
systemctl disable firewalld.service

https://github.com/MultiQC

https://github.com/ranger/ranger

https://zhuanlan.zhihu.com/p/34369349

the_plot <- function()
{
  x <- seq(0, 1, length.out = 100)
  y <- pbeta(x, 1, 10)
  plot(
    x,
    y,
    xlab = "False Positive Rate",
    ylab = "Average true positive rate",
    type = "l"
  )
}

png(
  "test.png",
  width     = 3.25,
  height    = 3.25,
  units     = "in",
  res       = 1200,
  pointsize = 4
)
par(
  mar      = c(5, 5, 2, 2),
  xaxs     = "i",
  yaxs     = "i",
  cex.axis = 2,
  cex.lab  = 2
)
the_plot()
dev.off()

当然，更好的解决方案是放弃这种基本的图形和使用一个系统，将处理你的分辨率缩放。例如，

library(ggplot2)

ggplot_alternative <- function()
{
  the_data <- data.frame(
    x <- seq(0, 1, length.out = 100),
    y = pbeta(x, 1, 10)
  )

ggplot(the_data, aes(x, y)) +
    geom_line() +
    xlab("False Positive Rate") +
    ylab("Average true positive rate") +
    coord_cartesian(0:1, 0:1)
}

ggsave(
  "ggtest.png",
  ggplot_alternative(),
  width = 3.25,
  height = 3.25,
  dpi = 1200
)

1199 git clone https://github.com/GPZ-Bioinfo/tmap.git
1200 cd tmap
1201 ll
1202 python setup.py install
1203 ll
1204 cd ../
1205 rm -rf tmap
1206 deactivate
1207 rmvirtualenv tmap_ENV
1208 mkvirtualenv -p /usr/bin/python3.4m tmap_ENV

pip3 install pypiwin32

conda install scipy

sudo pip3 install matplotlib

Workspace loaded from ~/.RData]

> setwd("/home/shenzy/work/beast/51samples")
> library(seqinr)
> data=read.fasta("51strain_core_gene_alignment.aln")
> library(ape)

Attaching package: ‘ape’

The following objects are masked from ‘package:seqinr’:

    as.alignment, consensus

> write.nexus.data(data,file="51strain_core_gene_alignment.aln.nexus", format="DNA")