dbscan聚类算法的R实现

2019年3月21日 255次阅读 来源: 聚类算法

首先，先讲下需要解决的问题：

问题：挑选出了一条染色体上的一些gene位点，用dbscan算法检查下这些基因在位置上有没有聚集。

输入文件：（ID,start,end） 

gene0001        1       1323
gene0002        1483    2619
gene0003        2580    4889
gene0009        14089   14664

PS：基因已按照起始位置排序


DBSCAN核心思想：如果一个点，在距它Eps的范围内有不少于MinPts个点，则该点就是核心点。核心和它Eps范围内的邻居形成一个簇。在一个簇内如果出现多个点都是核心点，则以这些核心点为中心的簇要合并。


R代码实现：

<pre name="code" class="plain">args<-commandArgs(TRUE)
dat<-read.table(args[1],sep="\t",header=F,row.names=1)
outfile<-args[2]
dis<-function(start1,end1,start2,end2){ ###计算两点之间距离
        start_p=max(start1,start2);
        end_p=min(end1,end2);
        if (start_p<end_p){return(0);}
        else{
         dis=start_p-end_p;
         return(dis);
        }
}

center_cluster<-function(dis,center,n,e){  #####寻找核心点和范围内的邻居，e（Eps）,n(MinPts)
        cluster<-center;
        while(TRUE){
                if (min(cluster)!=1 && dis[center,min(cluster)-1]<=e){
                        left=dis[center,min(cluster)-1];
                }
                else{
                        left=10000;
                }
               if (max(cluster)!=nrow(dis) && dis[max(cluster)+1,center]<=e){
                        right= dis[max(cluster)+1,center];
                }
                else{
                        right=10000;
                }
                if (left >right){cluster<-c(cluster,max(cluster)+1);}
                else if (left <right){cluster<-c(min(cluster)-1,cluster);}
                else if (left==right &&left!=10000){cluster<-c(min(cluster)-1,cluster);}
                else{break;}
        }

        if (length(cluster)>=n){
                return (cluster);
        }
        else{
                return (NULL);
        }
}

region_cluster<-function(dis,n,e){  ###将有overlap的核心簇连接成region
        region<-rep(0,nrow(dis));
        r<-1;
        for (i in 1:nrow(dis)){
                cluster=center_cluster(dis,i,n,e);
                if (length(cluster)){
                        if (max(region[cluster]!=0)){
                                region[cluster]=max(region[cluster]);
                        }
                        else{
                                region[cluster]=r;
                                r<-r+1;
                        }
                }
        }
        return(region);
}

gene_dis<-matrix(0,nrow(dat),nrow(dat)); 
for (i in 1:nrow(dat)){
        for (j in 1:i){
                gene_dis[i,j]=dis(dat[i,1],dat[i,2],dat[j,1],dat[j,2]);
        }
}

result<-as.data.frame(region_cluster(gene_dis,6,3000));
rownames(result)<-rownames(dat);
colnames(result)<-"cluster";

if ( max(result$cluster)>0 ){
        out<- matrix(0,max(result$cluster),2)
        for (i in 1:max(result$cluster)){
                out[i,]=c( dat[min(which(result$cluster==i)),1],dat[max(which(result$cluster==i)),2])
        }
        out<- as.data.frame(out)
        colnames(out)<-c("start","end")
        rownames(out)<-paste("region",rownames(out),sep="")
        write.table(out,outfile,quote=F, col.names = NA,sep="\t")
}else{
        write.table(t(c("","start","end")),outfile,quote=F,sep="\t",col.names=F,row.names=F) 
}


输出格式：
        start   end
 region1  11111  12000
region2  12001 120002

dbscan已有R包，写的代码没有卵用，但还是祭奠下~~~ 《dbscan聚类算法的R实现》
     原文作者：聚类算法
     原文地址: https://blog.csdn.net/hy_3210/article/details/50595158
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。 
  点赞