function [h,p,SIC_rel_reps] = plotFamilyBarGraphComparison(fam2plot, legend2plot, SIC_idx, L_taxonomy, L_relative_abundance_forexp) 

idx_fam = find(~cellfun(@isempty,strfind(L_taxonomy{5},fam2plot),'un',1)); % find certain family
n = length(SIC_idx);

% rows: SICS
% columns: replicates
for i=1:n
    % create vectors for relative abundance
    % do 3 replicates for now
    SIC_rel_reps(i,:) = L_relative_abundance_forexp{5}(SIC_idx(i):SIC_idx(i)+3, idx_fam)';
    SIC_rel_mean(i) = mean(SIC_rel_reps(i,:));
    SIC_rel_std(i) = std(SIC_rel_reps(i,:));
end

d_rel = SIC_rel_mean;
s_rel = SIC_rel_std;
% make data structures for bar plot
% d_rel = []; s_rel = [];
% for i = 1:7
%     % relative abundance
%     d_rel = [d_rel; [SIC_rel_mean(i)]];
%     s_rel = [s_rel; [SIC_rel_std(i)]]; 
% end

% plot relative abundance
figure
b = bar(d_rel); hold all;
set(gca,'xticklabel',legend2plot)
xtickangle(45)


% Plot the errorbars
errorbar(d_rel,s_rel,'k','linestyle','none', 'HandleVisibility','off'); hold all;
for i = 1:n
    scatter(repelem(i,4),SIC_rel_reps(i,:), [],'black', 'HandleVisibility','off'); hold all;
end

ylabel(strcat(fam2plot, ' rel. abundance'));
%xlim([0.3 3.5])
    
h = zeros(n,n);
p = zeros(n,n);
for i = 1:n
    for j = i:n
        [h(i,j), p(i,j)] = ttest2(SIC_rel_reps(i,1:4), SIC_rel_reps(j,1:4));
    end
end