%load file
load('diff_genes_across_clusters.mat')

%write a for loop that for each cluster...
for cluster = 1:num_clusters
    %sort the cells into those that belong in the cluster
    cluster_cells = (cluster_labels==cluster);
    %sort the cells into those that do NOT belong in the cluster
    noncluster_cells = (cluster_labels~=cluster);
    %initialize/reset data structures to store results for each gene
    %for every cluster for each iteration of the for loop
    cluster_genes=[];
    cluster_pvalues_per_gene=[];
    gene_counter=0;

    %for every gene...
    for gene = 1:num_genes
    %get the gene expression values for cells in cluster 1
        expr_cluster=A_log(cluster_cells,gene);
    %get the gene expression values for cells NOT in cluster 1
        expr_noncluster=A_log(noncluster_cells,gene);

    %check for differential gene expression between cells in and out
    %of cluster 1
        [~,p] = ttest2(expr_cluster,expr_noncluster,'Vartype','unequal');

    %store results only if...

    %1. p-value from t-test is significant...
        if p<threshold
          %AND 2. mean gene expression is higher in cells in cluster 1 vs noncluster 1
              if mean(expr_cluster)>mean(expr_noncluster)
              %count the gene
                  gene_counter=gene_counter+1;
              %remember the gene by its index number
                  cluster1_genes(gene_counter)=gene;
              %remember p-value for the gene
                  cluster1_pvalues_per_gene(gene_counter)=p;
              end
        end
    end

    %store the information for each cluster in a data structure called ge_clusters

    %get the number of the cluster
    ge_clusters(cluster).name = cluster;
    %get the indices of the differentially highly expressed genes for the cluster
    ge_clusters(cluster).genes = cluster_genes;
    %get the p-values of the differentially highly expressed genes for the cluster
    ge_clusters(cluster).pvalues = cluster_pvalues_per_gene;
end