dbscan聚类算法的R实现

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

问题：挑选出了一条染色体上的一些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包，写的代码没有卵用，但还是祭奠下~~~