diff --git a/.gitignore b/.gitignore
index 538c8c5..a25bbf7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 .DS_Store
 *~
+bin
+data
diff --git a/README.md b/README.md
index e2ef7d6..704f8c9 100644
--- a/README.md
+++ b/README.md
@@ -34,8 +34,9 @@ Directory structure:
   - bams: contains the aligned reads as BAM Files
   - logs: contains the log files of each alignment run
 - gene_counts: contains the read counts for each gene in the reference.
+- analysis: contains this repository (all the scripts used for analysis) as well as generated data, external programs, etc.
 
-## Steps
+## Preprocessing
 
 ### 1. MD5 checksum
 
@@ -152,3 +153,53 @@ Then we add a header to that file with labels and sample names:
 ```bash
 paste -d '\t' <(echo "GeneId") <(echo "Length") <(head -1 gene_counts_180117.txt.summary | cut -f2- | tr -s '\t' '\n' | awk 'OFS="\t"{split($1, a, "_"); print a[1]}' | paste -s -d '\t' -) | cat - gene_counts_180117_only_mapped.txt > gene_counts_180117_only_mapped_with_header.txt
 ```
+
+## Analysis
+
+### 1. Normalization
+
+Removed entries with RPKM values lower than 2 and applied deSeq normalisation:
+```
+analysis_rpkm_deSeq_normalization.m
+```
+
+### 2. Fold change
+
+Visualize fold change for each experiment time-range pair.
+```
+analysis_foldchange.m
+analysis_clustergram_foldchange.m
+```
+
+### 3. Binomial fit test
+Visualize binomial fit test for each experiment time-range pair.
+```
+analysis_binomial_test.m
+```
+
+### 4. Development of gene expression over time
+Analyze the changes in gene expression over time (for each experiment) for each gene. This classifies each gene into one of 27 possible combinations/conditions: ddu, uud, ndd, ... etc, where:
+- "d" = downregulated
+- "n" = no regulation (normal)
+- "u" = upregulated
+```
+analysis_expression_development.m
+```
+
+### 5. Gene ontology
+
+Run gene ontology annotation for each condition/combination file generated in the previous step. We use the perl script [GOAnalysis](https://software.scic.corp.brain.mpg.de/projects/MPIBR-Bioinformatics/GOAnalysis) by Georgi Tushev.
+
+```bash
+for file in target_genes_lists/*.txt; do file_name=$(basename "$file"); file_name="${file_name%.*}"; perl ../../bin/GOAnalysis/GOAnalysis.pl -obo annotations/go-basic.obo -ann annotations/gene_association.rgd.gz -background background_list.txt -target $file > "go_"$file_name".txt"; done
+```
+
+Get all unique GO terms available in the annotated files:
+```bash
+cat go_* | cut -f1 | sort -u > all_go_terms.txt
+```
+
+Create table with the P-value of each condition for each GO Term available:
+```
+analysis_gene_ontology.m
+```
diff --git a/analysis_binomial_test.m b/analysis_binomial_test.m
new file mode 100644
index 0000000..c66e548
--- /dev/null
+++ b/analysis_binomial_test.m
@@ -0,0 +1,27 @@
+clc
+clear variables;
+close all;
+
+% get records
+fRead = fopen('../gene_counts/gene_counts_180117_only_mapped_with_header.txt','r');
+% read first line (header)
+header = fgetl(fRead);
+header = strsplit(header);
+header(2) = []; % remove Length from header, not necessary for output file
+
+% read counts
+txt = textscan(fRead,'%s %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','delimiter','\t');
+fclose(fRead);
+symbols = txt{1}; % gene names
+length = txt{2}; % gene lengths
+counts = [txt{3:end}] + 1; % counts for each sample
+
+% remove entries with rpkm values lower than 2
+idx_threshold = calculate_rpkm(counts, length);
+counts(idx_threshold,:) = [];
+symbols(idx_threshold) = [];
+
+% create figures for each experiment time range (5min, 30min, 60min)
+stats05 = NBinFitTest(counts(:, [2,5]), counts(:, [1,4]), 'Constant', true);
+stats30 = NBinFitTest(counts(:, [7,10]), counts(:, [6,9]), 'Constant', true);
+stats60 = NBinFitTest(counts(:, [12,15]), counts(:, [11,14,16]), 'Constant', true);
diff --git a/analysis_clustergram_foldchange.m b/analysis_clustergram_foldchange.m
new file mode 100644
index 0000000..f48523a
--- /dev/null
+++ b/analysis_clustergram_foldchange.m
@@ -0,0 +1,49 @@
+clc
+clear variables;
+close all;
+
+% get records
+fRead = fopen('../gene_counts/gene_counts_180117_only_mapped_with_header.txt','r');
+% read first line (header)
+header = fgetl(fRead);
+header = strsplit(header);
+header(1:2) = []; % remove Length from header, not necessary for output file
+
+% read counts
+txt = textscan(fRead,'%s %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','delimiter','\t');
+fclose(fRead);
+symbols = txt{1}; % gene names
+length = txt{2}; % gene lengths
+counts = [txt{3:end}] + 1; % counts for each sample
+
+% remove entries with rpkm values lower than 2
+idx_threshold = calculate_rpkm(counts, length);
+counts(idx_threshold,:) = [];
+symbols(idx_threshold) = [];
+
+% create figures for each experiment time range (5min, 30min, 60min)
+stats05 = BinFoldChangeTest(counts(:, [2,5]), counts(:, [1,4]), 100, false);
+stats30 = BinFoldChangeTest(counts(:, [7,10]), counts(:, [6,9]), 100, false);
+stats60 = BinFoldChangeTest(counts(:, [12,15]), counts(:, [11,14,16]), 100, false);
+
+% idx_x contains genes meaninfully underexpressed
+% idx_y contains genes meaninfully overexpressed
+idx_common = stats05.idx_x | stats05.idx_y | stats30.idx_x | stats30.idx_y | stats60.idx_x | stats60.idx_y;
+M = [stats05.M, stats30.M, stats60.M];
+M = M(idx_common,:);
+
+R = tiedrank(M)./size(M,1);
+
+clustergram(R,'Standardize', 2,'cluster',1, 'ColumnLabels', {'05 min', '30 min', '60 min'});
+
+%{
+mtx = calculate_deSeq(counts);
+mtx = mtx(idx_common,:);
+idx_remove = [3,8,13];
+mtx(:,idx_remove) = [];
+header(idx_remove) = [];
+
+% disp(symbols(stats05.idx_y));
+R = tiedrank(mtx)./size(mtx,1);
+clustergram(R,'Standardize',2,'cluster',1)
+%}
\ No newline at end of file
diff --git a/analysis_expression_development.m b/analysis_expression_development.m
new file mode 100644
index 0000000..b798b71
--- /dev/null
+++ b/analysis_expression_development.m
@@ -0,0 +1,78 @@
+clc
+clear variables;
+close all;
+
+% get records
+fRead = fopen('../gene_counts/gene_counts_180117_only_mapped_with_header.txt','r');
+% read first line (header)
+header = fgetl(fRead);
+header = strsplit(header);
+header(1:2) = []; % remove Length from header, not necessary for output file
+
+% read counts
+txt = textscan(fRead,'%s %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','delimiter','\t');
+fclose(fRead);
+symbols = txt{1}; % gene names
+length = txt{2}; % gene lengths
+counts = [txt{3:end}] + 1; % counts for each sample
+
+% remove entries with rpkm values lower than 2
+idx_threshold = calculate_rpkm(counts, length);
+counts(idx_threshold,:) = [];
+symbols(idx_threshold) = [];
+
+% create genes background list
+fileID = fopen('data/gene_ontology/background_list.txt','w');
+formatSpec = '%s\n';
+[nrows,ncols] = size(symbols);
+for row = 1:nrows
+    fprintf(fileID,formatSpec,symbols{row,:});
+end
+fclose(fileID);
+
+% create figures for each experiment time range (5min, 30min, 60min)
+stats05 = BinFoldChangeTest(counts(:, [2,5]), counts(:, [1,4]), 100, false);
+stats30 = BinFoldChangeTest(counts(:, [7,10]), counts(:, [6,9]), 100, false);
+stats60 = BinFoldChangeTest(counts(:, [12,15]), counts(:, [11,14,16]), 100, false);
+
+M = [stats05.M, stats30.M, stats60.M];
+
+% get all gene expression combinations over time
+% 'd' means significantly downregulated
+% 'n' means no significant change in expression (normal)
+% 'u' means significantly upregulated
+% i.e.: u,n,n means gene upregulated in sample01, normal in sample02, normal in sample03
+
+% http://stackoverflow.com/questions/4165859/generate-all-possible-combinations-of-the-elements-of-some-vectors-cartesian-pr#4169488
+x = ['d', 'n', 'u']; % Set of possible gene expression states
+K = 3; % Length of each permutation
+s = struct('d', 'idx_x', 'n', 'idx_none', 'u', 'idx_y'); % reference of each state to corresponding index from fold change stats
+
+% Create all possible permutations (with repetition) of states stored in x
+C = cell(K, 1);  % Preallocate a cell array
+[C{:}] = ndgrid(x); % Create K grids of values
+y = cellfun(@(x){x(:)}, C); % Convert grids to column vectors
+y = [y{:}]; % Obtain all permutations
+
+[nrows,ncols] = size(y);
+exp = [];
+exp_header = {};
+for row = 1:nrows
+    exp_idx = (stats05.(s.(y(row,1))) & stats30.(s.(y(row,2))) & stats60.(s.(y(row,3))));
+    exp_geneIds = symbols(exp_idx, :);
+    combination_label = y(row,:);
+    % write gene IDs to file
+    fileName = sprintf('data/gene_ontology/%s.txt', combination_label);
+    fileHandler = fopen(fileName,'w');
+    formatSpec = '%s\n';
+    fprintf(fileHandler,formatSpec,exp_geneIds{:});
+    fclose(fileHandler);
+    % append expression combination label to header
+    exp_header{end+1} = combination_label;
+    % write index for all genes
+    exp = [exp exp_idx];
+end
+% create table with row and column labels
+exp_table = array2table(exp, 'RowNames', symbols, 'VariableNames', exp_header);
+
+
diff --git a/analysis_foldchange.m b/analysis_foldchange.m
new file mode 100644
index 0000000..b748054
--- /dev/null
+++ b/analysis_foldchange.m
@@ -0,0 +1,27 @@
+clc
+clear variables;
+close all;
+
+% get records
+fRead = fopen('../gene_counts/gene_counts_180117_only_mapped_with_header.txt','r');
+% read first line (header)
+header = fgetl(fRead);
+header = strsplit(header);
+header(1:2) = []; % remove Length from header, not necessary for output file
+
+% read counts
+txt = textscan(fRead,'%s %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','delimiter','\t');
+fclose(fRead);
+symbols = txt{1}; % gene names
+length = txt{2}; % gene lengths
+counts = [txt{3:end}] + 1; % counts for each sample
+
+% remove entries with rpkm values lower than 2
+idx_threshold = calculate_rpkm(counts, length);
+counts(idx_threshold,:) = [];
+symbols(idx_threshold) = [];
+
+% create figures for each experiment time range (5min, 30min, 60min)
+stats05 = BinFoldChangeTest(counts(:, [2,5]), counts(:, [1,4]), 100, true);
+stats30 = BinFoldChangeTest(counts(:, [7,10]), counts(:, [6,9]), 100, true);
+stats60 = BinFoldChangeTest(counts(:, [12,15]), counts(:, [11,14,16]), 100, true);
\ No newline at end of file
diff --git a/analysis_gene_ontology.m b/analysis_gene_ontology.m
new file mode 100644
index 0000000..5ffd00d
--- /dev/null
+++ b/analysis_gene_ontology.m
@@ -0,0 +1,31 @@
+clc
+clear variables;
+close all;
+
+% get all unique GO terms
+fRead = fopen('data/gene_ontology/all_go_terms.txt','r');
+ref_go_terms = textscan(fRead,'%s','delimiter','\t');
+ref_go_terms = ref_go_terms{1};
+fclose(fRead);
+
+% get go annotation files
+annot_files = dir('data/gene_ontology/go_*.txt');
+
+go_matrix = zeros(size(ref_go_terms, 1), size(annot_files, 1));
+for i = 1:size(annot_files, 1)
+    file = annot_files(i);
+    fRead = fopen(sprintf('data/gene_ontology/%s',file.name),'r');
+    content = textscan(fRead,'%s %s %s %n %n %n %n %n','delimiter','\t');
+    qry_go_terms = content{1};
+    p_values = content{end};
+    fclose(fRead);
+    
+    [~, idx_ref, idx_qry] = intersect(ref_go_terms, qry_go_terms);
+    
+    go_matrix(idx_ref, i) = p_values(idx_qry);
+
+end
+
+go_table = array2table(go_matrix, 'RowNames', ref_go_terms, 'VariableNames', cellfun(@(s) strtok(s, '.') ,{annot_files.name}, 'UniformOutput', false));
+
+clustergram(go_matrix,'Standardize', 2,'cluster',1, 'ColumnLabels', cellfun(@(s) strtok(s, '.') ,{annot_files.name}, 'UniformOutput', false));
diff --git a/analysis_rpkm_deSeq_normalization.m b/analysis_rpkm_deSeq_normalization.m
new file mode 100644
index 0000000..24c9b5a
--- /dev/null
+++ b/analysis_rpkm_deSeq_normalization.m
@@ -0,0 +1,37 @@
+clc
+clear variables;
+close all;
+
+% get records
+fRead = fopen('../gene_counts/gene_counts_180117_only_mapped_with_header.txt','r');
+% read first line (header)
+header = fgetl(fRead);
+header = strsplit(header);
+header(2) = []; % remove Length from header, not necessary for output file
+
+% read counts
+txt = textscan(fRead,'%s %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n %n','delimiter','\t');
+fclose(fRead);
+symbols = txt{1}; % gene names
+length = txt{2}; % gene lengths
+counts = [txt{3:end}]; % counts for each sample
+
+% remove entries with rpkm values lower than 2
+idx_threshold = calculate_rpkm(counts, length);
+counts(idx_threshold,:) = [];
+symbols(idx_threshold) = [];
+
+% apply deSeq normalization to counts
+normalized_counts = calculate_deSeq(counts);
+
+% get normalized list of gene counts
+normalized_cell = [symbols num2cell(normalized_counts)];
+% print list to text file
+normalized_table = cell2table(normalized_cell, 'VariableNames', header);
+% create file
+writetable(normalized_table, 'rpkm_and_deSeq_normalized_data.txt', 'Delimiter','\t');
+
+% plot correlation matrix
+header(1) = []; % remove gene ids column from header
+plotCorrelationMatrix(normalized_counts, header, [16,16], 'correlation_matrix_rpkmBigger2_and_deSeq');
+
diff --git a/figures/binomial_fit_test_s12-s15_s11-s14-s16.fig b/figures/binomial_fit_test_s12-s15_s11-s14-s16.fig
new file mode 100644
index 0000000..5071ba1
Binary files /dev/null and b/figures/binomial_fit_test_s12-s15_s11-s14-s16.fig differ
diff --git a/figures/binomial_fit_test_s2-s5_s1-s4.fig b/figures/binomial_fit_test_s2-s5_s1-s4.fig
new file mode 100644
index 0000000..0bb86c1
Binary files /dev/null and b/figures/binomial_fit_test_s2-s5_s1-s4.fig differ
diff --git a/figures/binomial_fit_test_s7-s10_s6-s9.fig b/figures/binomial_fit_test_s7-s10_s6-s9.fig
new file mode 100644
index 0000000..558b437
Binary files /dev/null and b/figures/binomial_fit_test_s7-s10_s6-s9.fig differ
diff --git a/figures/clustergram_fold_change.fig b/figures/clustergram_fold_change.fig
new file mode 100644
index 0000000..9f5afab
Binary files /dev/null and b/figures/clustergram_fold_change.fig differ
diff --git a/figures/correlation_matrix_rpkmBigger2.png b/figures/correlation_matrix_rpkmBigger2.png
new file mode 100644
index 0000000..fed19d5
Binary files /dev/null and b/figures/correlation_matrix_rpkmBigger2.png differ
diff --git a/figures/correlation_matrix_rpkmBigger2_and_deSeq.png b/figures/correlation_matrix_rpkmBigger2_and_deSeq.png
new file mode 100644
index 0000000..335d163
Binary files /dev/null and b/figures/correlation_matrix_rpkmBigger2_and_deSeq.png differ
diff --git a/figures/fold_change_s12-s15_s11-s14-s16.fig b/figures/fold_change_s12-s15_s11-s14-s16.fig
new file mode 100644
index 0000000..3bd71cb
Binary files /dev/null and b/figures/fold_change_s12-s15_s11-s14-s16.fig differ
diff --git a/figures/fold_change_s2-s5_s1-s4.fig b/figures/fold_change_s2-s5_s1-s4.fig
new file mode 100644
index 0000000..8d00188
Binary files /dev/null and b/figures/fold_change_s2-s5_s1-s4.fig differ
diff --git a/figures/fold_change_s7-s10_s6-s9.fig b/figures/fold_change_s7-s10_s6-s9.fig
new file mode 100644
index 0000000..a189d92
Binary files /dev/null and b/figures/fold_change_s7-s10_s6-s9.fig differ
diff --git a/figures/gene_ontology_clustergram.fig b/figures/gene_ontology_clustergram.fig
new file mode 100644
index 0000000..02a011c
Binary files /dev/null and b/figures/gene_ontology_clustergram.fig differ
diff --git a/helpers/BinFoldChangeTest.m b/helpers/BinFoldChangeTest.m
new file mode 100644
index 0000000..eca4b94
--- /dev/null
+++ b/helpers/BinFoldChangeTest.m
@@ -0,0 +1,57 @@
+function [stats] = BinFoldChangeTest(X,Y,nbins,draw)
+
+    nrm = calculate_deSeq([X,Y]);
+    Xnrm = nrm(:,1:size(X,2));
+    Ynrm = nrm(:,(size(X,2)+1):end);
+    A = log10(mean(nrm, 2));
+    M = log2(mean(Ynrm,2)) - log2(mean(Xnrm,2));
+
+    % create bins
+    idx_bin = zeros(size(A,1),1);
+    idx_bin(1:nbins:end) = 1;
+    idx_bin = cumsum(idx_bin);
+
+    % resort
+    [~, idx_sort] = sort(A,'descend');
+    idx_rev(idx_sort) = (1:size(idx_sort,1))';
+    idx_res = false(size(A,1),1);
+    Z = zeros(size(A,1),1);
+    for k = 1 : max(idx_bin)
+        idx_now = idx_bin == k;
+        M_now =  M(idx_now(idx_rev));
+        Z_now = zscore(M_now);
+        idx_good = abs(Z_now) >= 1.65;
+        idx_res(idx_now(idx_rev)) = idx_good;
+        Z(idx_now(idx_rev)) = Z_now;
+    end
+    
+    idx_x = idx_res & (M < median(M)) & (A >= 1);
+    idx_y = idx_res & (M >= median(M)) & (A >= 1);
+    idx_none = ~(idx_x | idx_y);
+
+    if draw
+        
+        figure('Color','w');
+        hold(gca,'on');
+        plot(A(idx_none),M(idx_none), '.', 'color', [.65,.65,.65]);
+        plot(A(idx_x), M(idx_x), '.', 'color', [30,144,255]./255);
+        plot(A(idx_y), M(idx_y), '.', 'color', [148,0,211]./255);
+        plot([min(A),max(A)],[median(M),median(M)],'k');
+        hold(gca,'off');
+        
+        set(gca, 'box','off',...
+                 'xlim',[min(A), max(A)],...
+                 'ylim',[-max(abs(M)),max(abs(M))]);
+        
+    end
+    
+    % pack in structure
+    stats.nrm = nrm;
+    stats.A = A;
+    stats.M = M;
+    stats.Z = Z;
+    stats.idx_x = idx_x;
+    stats.idx_y = idx_y;
+    stats.idx_none = idx_none;
+    
+end
\ No newline at end of file
diff --git a/helpers/NBinFitTest.m b/helpers/NBinFitTest.m
new file mode 100644
index 0000000..fd5242d
--- /dev/null
+++ b/helpers/NBinFitTest.m
@@ -0,0 +1,40 @@
+function [stats] = NBinFitTest(X, Y, vartest, plot_tf)
+
+    nrm = calculate_deSeq([X,Y]);
+    Xnrm = nrm(:,1:size(X,2));
+    Ynrm = nrm(:,(size(X,2)+1):end);
+    A = log10(mean(nrm, 2));
+    M = log2(mean(Ynrm,2)) - log2(mean(Xnrm,2));
+    
+    t = nbintest(X, Y, 'VarianceLink', vartest);
+    idx_x = (t.pValue <= 0.05) & (M < median(M));
+    idx_y = (t.pValue <= 0.05) & (M >= median(M));
+    idx_none = ~(idx_x | idx_y);
+
+    if plot_tf
+        
+        figure('Color','w');
+        hold(gca,'on');
+        plot(A(idx_none),M(idx_none), '.', 'color', [.65,.65,.65]);
+        plot(A(idx_x), M(idx_x), '.', 'color', [30,144,255]./255);
+        plot(A(idx_y), M(idx_y), '.', 'color', [148,0,211]./255);
+        plot([min(A),max(A)],[median(M),median(M)],'k');
+        hold(gca,'off');
+        
+        set(gca, 'box','off',...
+                 'xlim',[min(A), max(A)],...
+                 'ylim',[-max(abs(M)),max(abs(M))]);
+        
+    end
+    
+    % pack in structure
+    stats.nrm = nrm;
+    stats.A = A;
+    stats.M = M;
+    stats.p = t.pValue;
+    stats.idx_x = idx_x;
+    stats.idx_y = idx_y;
+    stats.idx_none = idx_none;
+    
+end
+
diff --git a/helpers/calculate_deSeq.m b/helpers/calculate_deSeq.m
new file mode 100644
index 0000000..3880f86
--- /dev/null
+++ b/helpers/calculate_deSeq.m
@@ -0,0 +1,18 @@
+function [norm_counts,size_factors] = calculate_deSeq( raw_counts )
+    % applying deSeq normalization 
+    % assuming equal median between columns
+
+    % calculate geometrical means
+    geo_means = exp(mean(log(raw_counts),2));
+
+    % calculate ratios
+    ratios = bsxfun(@rdivide, raw_counts(geo_means > 0, :), geo_means(geo_means > 0));
+
+    % calculate size_factors
+    size_factors = median(ratios, 1);
+
+    % calculate normalized counts
+    norm_counts = bsxfun(@rdivide, raw_counts, size_factors);
+
+
+end
\ No newline at end of file
diff --git a/helpers/calculate_rpkm.m b/helpers/calculate_rpkm.m
new file mode 100644
index 0000000..f9814f2
--- /dev/null
+++ b/helpers/calculate_rpkm.m
@@ -0,0 +1,7 @@
+function idx_threshold = calculate_rpkm( raw_counts, length )
+
+    rpkm = bsxfun(@rdivide, 1e6.*raw_counts(:,2:end), sum(raw_counts(:,2:end),1));
+    rpkm = bsxfun(@rdivide, 1e3.*rpkm, length);
+
+    idx_threshold = all(rpkm < 2, 2);
+end
\ No newline at end of file
diff --git a/helpers/plotCorrelationMatrix.m b/helpers/plotCorrelationMatrix.m
new file mode 100644
index 0000000..7073ad5
--- /dev/null
+++ b/helpers/plotCorrelationMatrix.m
@@ -0,0 +1,99 @@
+function [] = plotCorrelationMatrix(cnts,label,fonts, file_name)
+% plots correlation matrix
+
+    % figure size
+    fig_h = 800;
+    fig_w = 800;
+
+    % get dimensions
+    dim_w = size(cnts,2);
+    
+    % get borders
+    k = 1;
+    left = 0.04;
+    bottom = 0.04;
+    margin = 0.01;
+    width = (1 - left - dim_w*margin)/dim_w;
+    height = width;
+    
+    % create figure
+    hFig = figure('Color','w',...
+       'Tag','hMainWindow',...
+       'Position',[200,100,fig_w,fig_h],...
+       'MenuBar','none',...
+       'ToolBar','none',...
+       'Name','SeqStatAnalysis',...
+       'NumberTitle','off');
+    hAx = zeros(dim_w*dim_w,1);
+
+    % create idx matrix
+    idx_all = (1:dim_w*dim_w);
+    idx_mtx = reshape(idx_all, dim_w, dim_w);
+    idx_plot = tril(idx_mtx,-1);
+    idx_plot(idx_plot == 0) = [];
+    idx_rsq = triu(idx_mtx,1);
+    idx_rsq(idx_rsq == 0) = [];
+
+    for i = 1 : dim_w
+        for j = 1 : dim_w
+            hAx(k) = axes('Parent',hFig,'Box','on',...
+                      'XTick',[],'YTick',[],...
+                      'Color','w','Units','normalized',...
+                      'Position',[left,bottom,width,height]);
+        
+            y = log2(cnts(:,j)+1);
+            x = log2(cnts(:,i)+1);
+                  
+            % plot data
+            if any(idx_plot == k)
+                hold(hAx(k),'on');
+                plot(x,y,'.','Color',[.8,.8,.8]);
+                plot([0,20],[0,20],'k');
+                hold(hAx(k),'off');
+                set(hAx(k),'Box','off','Layer','top',...
+                       'XTick',[],'XLim',[-2,20],...
+                       'YTick',[],'YLim',[-2,20]);
+            end
+        
+            % plot histogram
+            if (i == j)
+                [f,x_bin] = hist(x,10);
+                bar(hAx(k),x_bin,f,0.2,'EdgeColor','none','FaceColor',[.4,.4,.4]);
+                set(hAx(k),'Box','off','Layer','top',...
+                       'XTick',[],'XLim',[0,20],...
+                       'YTick',[],'YLim',[0,10000]);
+            end
+            % plot rsq
+            if any(idx_rsq == k)
+                rsq = corr(x,y);
+                text(0.5,0.5,sprintf('%.4f',rsq),'Parent',hAx(k),...
+                'VerticalAlignment','middle','HorizontalAlignment','center',...
+                'FontSize',fonts(2));
+                set(hAx(k),'XTick',[],'XLim',[0,1],'YTick',[],'YLim',[0,1]);
+            end
+        
+        
+            % add X-Label
+            if(i == 1)
+                xlabel(hAx(k),label{j},'FontSize',fonts(1));
+            end
+             
+            % add Y-Label
+            if(j == 1)
+                ylabel(hAx(k),label{i},'FontSize',fonts(1));
+            end
+    
+            % update coordinates
+            left = left + margin + width;
+            k = k + 1;
+        
+        end
+
+        % update coordinates
+        left = 0.04;
+        bottom = bottom + margin + height;
+    end
+
+    print(hFig,'-dpng','-r300',[file_name '.png']);
+
+end
\ No newline at end of file