Энтеротипирование Driver Genera PCA сюжеты - PullRequest
Новая
       37

Энтеротипирование Driver Genera PCA сюжеты

0 голосов
/ 08 апреля 2020

Я в основном следовал учебному пособию при построении сюжета для энтеротипирования, как описано:

https://enterotype.embl.de/

Но после анализа между классами, как я могу определить роды драйверов для каждого кластера? У меня есть таблица, отформатированная с помощью Sample ID в столбцах, OTU в строках вместе с относительным изобилием.
Вот мой скрипт: 

data=read.table("d:/enterotype_aspiratedata.txt", header=T, row.names=1, dec=".", sep="\t")

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=NULL

for (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)

noise.removal <- function(dataframe, percent=0.01, top=NULL){
  dataframe->Matrix
  bigones <- rowSums(Matrix)*100/(sum(rowSums(Matrix))) > percent 
  Matrix_1 <- Matrix[bigones,]
  print(percent)
  return(Matrix_1)
}
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) 

cat("\nplot", k, "clusters between-class analysis (BCA).\n") 

bca.df <- as.data.frame(obs.bet$ls)
colnames(bca.df)[1:2] <- c("PC1", "PC2")
bca.df$cluster <- data.cluster

labs <- rownames(obs.bet$ls)

s.class(obs.bet$ls,fac=as.factor(data.cluster),
        cstar = 1,axesell = TRUE, cellipse = 1.5, grid=F,col=rainbow(3),
        cpoint = 0)
points(obs.bet$ls[grep("N\\.", rownames(obs.bet$ls)),], pch = 15, col = "#FFE4C4")
points(obs.bet$ls[grep("C\\.", rownames(obs.bet$ls)),], pch = 16, col = "#008B8B")
points(obs.bet$ls[grep("A\\.", rownames(obs.bet$ls)),], pch = 19, col = "#006400")

cbind(data.cluster,obs.bet$ls)
obs.pca$c1

#plot 2
obs.pcoa=dudi.pco(data.dist, scannf=F, nf=3)
s.class(obs.pcoa$li, fac=as.factor(data.cluster), sub="Principal coordiante analysis",
        cstar = 1,axesell = TRUE, cellipse = 1.5, grid=F,col=rainbow(3),
        cpoint = 0)

points(obs.bet$ls[grep("N\\.", rownames(obs.bet$ls)),], pch = 0, col = "#FFE4C4")
points(obs.bet$ls[grep("C\\.", rownames(obs.bet$ls)),], pch = 1, col = "#008B8B")
points(obs.bet$ls[grep("A\\.", rownames(obs.bet$ls)),], pch = 19, col = "#006400")

введите описание изображения здесь

...