function [h,p] = plotFamilyBarGraphComparisonByDrug(fam2plot, legend2plot, SIC1_idx, SIC2_idx, SIC3_idx, SIC4_idx, L_taxonomy, L_5_relative_abundance_cleaned_forexp, L_5_abs_abundance_forexp) 

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

% 14 mb plus BtWT
for i=1:3
    % create vectors for relative abundance
    % without drug, +saxa, +sita
    SIC1_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp([SIC1_idx(i),SIC1_idx(i)+3,SIC1_idx(i)+6], idx_fam)';
    SIC2_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp([SIC2_idx(i),SIC2_idx(i)+3,SIC2_idx(i)+6], idx_fam)';
    SIC3_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp([SIC3_idx(i),SIC3_idx(i)+3,SIC3_idx(i)+6], idx_fam)';
    SIC4_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp([SIC4_idx(i),SIC4_idx(i)+3,SIC4_idx(i)+6], idx_fam)';
    
    % create vectors for absolute abundance
    SIC1_abs_reps(i,:) = L_5_abs_abundance_forexp([SIC1_idx(i),SIC1_idx(i)+3,SIC1_idx(i)+6], idx_fam)';
    SIC2_abs_reps(i,:) = L_5_abs_abundance_forexp([SIC2_idx(i),SIC2_idx(i)+3,SIC2_idx(i)+6], idx_fam)';
    SIC3_abs_reps(i,:) = L_5_abs_abundance_forexp([SIC3_idx(i),SIC3_idx(i)+3,SIC3_idx(i)+6], idx_fam)';
    SIC4_abs_reps(i,:) = L_5_abs_abundance_forexp([SIC4_idx(i),SIC4_idx(i)+3,SIC4_idx(i)+6], idx_fam)';
end

% mean relative abundance 
SIC1_rel_mean = mean(SIC1_rel_reps);
SIC2_rel_mean = mean(SIC2_rel_reps);
SIC3_rel_mean = mean(SIC3_rel_reps);
SIC4_rel_mean = mean(SIC4_rel_reps);

% std relative abundance
SIC1_rel_std = std(SIC1_rel_reps);
SIC2_rel_std = std(SIC2_rel_reps);
SIC3_rel_std = std(SIC3_rel_reps);
SIC4_rel_std = std(SIC4_rel_reps);

% mean relative abundance 
SIC1_abs_mean = mean(SIC1_abs_reps);
SIC2_abs_mean = mean(SIC2_abs_reps);
SIC3_abs_mean = mean(SIC3_abs_reps);
SIC4_abs_mean = mean(SIC4_abs_reps);

% std relative abundance
SIC1_abs_std = std(SIC1_abs_reps);
SIC2_abs_std = std(SIC2_abs_reps);
SIC3_abs_std = std(SIC3_abs_reps);
SIC4_abs_std = std(SIC4_abs_reps);
    
% make data structures for bar plot
d_rel = []; s_rel = []; d_abs = []; s_abs = [];
for i = 1:3
    % relative abundance
    d_rel = [d_rel; [SIC1_rel_mean(i) SIC2_rel_mean(i) SIC3_rel_mean(i) SIC4_rel_mean(i)]];
    s_rel = [s_rel; [SIC1_rel_std(i) SIC2_rel_std(i) SIC3_rel_std(i) SIC4_rel_std(i)]];
    
    % absolute abundance
    d_abs = [d_abs; [SIC1_abs_mean(i) SIC2_abs_mean(i) SIC3_abs_mean(i) SIC4_abs_mean(i)]];
    s_abs = [s_abs; [SIC1_abs_std(i) SIC2_abs_std(i) SIC3_abs_std(i) SIC4_abs_std(i)]];
end

% plot relative abundance
figure
b = bar(d_rel, 'grouped'); hold all;

% Calculate the number of groups and number of bars in each group
[ngroups,nbars] = size(d_rel);
% Get the x coordinate of the bars
x = nan(nbars, ngroups);
for i = 1:nbars
    x(i,:) = b(i).XEndPoints;
end
% Plot the errorbars
errorbar(x',d_rel,s_rel,'k','linestyle','none'); hold all;
xticklabels({'control','+saxa','+sita'})
%xlabel('passages');
ylabel(strcat(fam2plot, ' rel. abundance'));
set(gca,'FontSize',14)
ylim([0 0.3])
xlim([0.3 3.5])
legend(legend2plot, 'location', 'northeast');

% plot absolute abundance
figure
b = bar(d_abs, 'grouped'); hold all;

[ngroups,nbars] = size(d_abs); % Calculate the number of groups and number of bars in each group
x = nan(nbars, ngroups); % Get the x coordinate of the bars

for i = 1:nbars
    x(i,:) = b(i).XEndPoints;
end

% Plot the errorbars
errorbar(x',d_abs,s_abs,'k','linestyle','none'); hold all;
xticklabels({'control','+saxa','+sita'})
%xlabel('passages');
ylabel(strcat(fam2plot, ' abs. abundance'));
set(gca,'FontSize',14)
ylim([0 10])
xlim([0.3 3.5])
legend(legend2plot, 'location', 'northeast');



h = [];
conds = ["control","saxa","sita"];
for i = 1:3
    varname = strcat(conds(i),'rel');
    [h(1),p.(varname)(1,2)] = ttest(SIC1_rel_reps(1:3,i), SIC2_rel_reps(1:3,i));
    [h(2),p.(varname)(1,3)] = ttest(SIC1_rel_reps(1:3,i), SIC3_rel_reps(1:3,i));
    [h(3),p.(varname)(1,4)] = ttest(SIC1_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    [h(7),p.(varname)(2,3)] = ttest(SIC2_rel_reps(1:3,i), SIC3_rel_reps(1:3,i));
    [h(8),p.(varname)(2,4)] = ttest(SIC2_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    [h(9),p.(varname)(3,4)] = ttest(SIC3_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    
    varname = strcat(conds(i), 'abs');
    [h(4),p.(varname)(1,2)] = ttest(SIC1_abs_reps(1:3,i), SIC2_abs_reps(1:3,i));
    [h(5),p.(varname)(1,3)] = ttest(SIC1_abs_reps(1:3,i), SIC3_abs_reps(1:3,i));
    [h(6),p.(varname)(1,4)] = ttest(SIC1_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
    [h(10),p.(varname)(2,3)] = ttest(SIC2_abs_reps(1:3,i), SIC3_abs_reps(1:3,i));
    [h(11),p.(varname)(2,4)] = ttest(SIC2_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
    [h(12),p.(varname)(3,4)] = ttest(SIC3_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
end