function [bestclust, GeneGroup, GeneGroupList]=SelectCell(matrix, genelist, clsn, reps, conc, genelist_TFhistone, bestclust_TF, bestclust_Histone)
% matrix: gene expression profile (genes x conditions)
% genelist: genelist of matrix
% clsn: the number of cluster to be obtained
% reps: replication time (eg. 1000)
% conc: concentration parameter of vMF distribution (eg. 10)
% genelist_TFhistone: genelist of TF-HM
% bestclust_TF: cluster ID on TF-binding (genes x cluster ID)
% bestclust_Histone: cluster ID on histone modification (genes x cluster ID)


GeneGroup=[];
GeneGroupList={};

ID_mvsKB=[];
ID_KBvsm=[];
n=length(genelist_TFhistone);
m=length(genelist);
for i=1:n
    for j=1:m
        if strcmpi(genelist_TFhistone(i), genelist(j))==1
            ID_mvsKB=[ID_mvsKB;j];
            ID_KBvsm=[ID_KBvsm;i];
        end
    end
end

genelist_part=genelist_TFhistone(ID_KBvsm);
bestclust_TF_part=bestclust_TF(ID_KBvsm);
bestclust_Histone_part=bestclust_Histone(ID_KBvsm);

zeromatrix=zeros(n,10);
for i=1:n
    zeromatrix(i,bestclust_TF(i))=1;
end
table_H=zeromatrix;

zeromatrix=zeros(n,10);
for i=1:n
    zeromatrix(i,bestclust_Histone(i))=1;
end
table_KB=zeromatrix;

table_H_part=table_H(ID_KBvsm,:);
table_KB_part=table_KB(ID_KBvsm,:);

M=matrix(ID_mvsKB,:);
W=M*M';

norm=normvector(W,clsn);

[bestclust] = movmf_m(W,norm,clsn,reps,conc);

for i=1:clsn
    
    disp('----------------------------------------------------------------------')
    disp('cluster ID is')
    disp(i)
    table_clsn=table_H_part(bestclust==i,:)'*table_KB_part(bestclust==i,:);
    %image(table_clsn*4)
    ave_clsn=mean(reshape(table_clsn,100,1));
    std_clsn=std(reshape(table_clsn,100,1));
    nM_clsn=(table_clsn-repmat(ave_clsn,10,10))./repmat(std_clsn,10,10);
    tcdf_clsn=tcdf(nM_clsn,99);
    
    selected_clsn=find(tcdf_clsn>0.99);
    TFID_selected_clsn=rem(selected_clsn,10);
    just=find(TFID_selected_clsn==0);
    TFID_selected_clsn(just)=10; 
    
    hisID_selected_clsn=fix(selected_clsn./10)+1;
    if isempty(just)==0
        hisID_selected_clsn(just)=fix(selected_clsn(just)./10);
    end
    
    disp('IDs of selected cells are')
    display('    TF  Histone  Expression')
    i_clsn=repmat(i,size(TFID_selected_clsn),1);
    obtained=[TFID_selected_clsn hisID_selected_clsn i_clsn];
    
    disp(obtained)
    GeneGroup=[GeneGroup;obtained];
    
    [a,b]=size(obtained);
    
    table=[bestclust_TF_part bestclust_Histone_part bestclust];

    for I=1:a
        id_GeneGroup=find(table(:,1)==obtained(I,1) & table(:,2)==obtained(I,2) & table(:,3)==obtained(I,3));
        genes={genelist_part(id_GeneGroup)};
        GeneGroupList=[GeneGroupList genes];
    end
      
    %pause
end



