import math
import os
import random
import subprocess
import colorcet as cc
import numpy
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import skbio
import skbio.stats.ordination
from sklearn.decomposition import PCA


#Path to current working directory.
PATH = os.getcwd()

#creates file name and title used when generating graphs
def make_filename(regex):
    filename = ""
    if regex.find('18') != -1:
        filename = filename + "2018 "
    else:
        filename = filename + "2019 "

    if regex.find('dna') != -1:
        filename = filename + "DNA "
    else:
        filename = filename + "RNA "

    if regex.find('16S') != -1:
        filename = filename + "16S "
    else:
        filename = filename + "28S "
    if regex.find('ctrl') != -1:
        filename = filename + "ctrl štartér"
    elif regex.find('lasa') != -1:
        filename = filename + "lasa štartér"
    else:
        filename = filename + "mesa štartér"
    return filename

#groups ASV rows based on same taxon
def group_ASVs_by_taxon(markerGene):

    with open(PATH + '/winezymes_'+ markerGene +'_18-19_dna-rna.tsv',) as file:
        reader = pd.read_csv(file, sep='\t')
        reader.drop(index=0, inplace=True)
        for i in reader.columns[4:]:
            reader[i] = pd.to_numeric(reader[i])
        filteredTable = pd.concat([reader.iloc[:, 0:3], reader.filter(regex='^18|^19')], axis=1, join='inner')
        filteredTable['Taxon'] = filteredTable['Taxon'].astype(str) + ';'
        filteredTable = filteredTable.drop(['Sequence'], axis=1)
        groupedDF = pd.DataFrame()
        usedList = list()
        indexList = list()

        for x in range(0, len(filteredTable)):
            string = str(filteredTable.iloc[x]['Taxon'])
            count = 0
            regex = ''
            flag = False
            for i in range(0, len(string)):
                if string[i] == ';':
                    count += 1
                if count == 9:
                    regex = string[:i]
                    count = 0
                    flag = True
                    if regex not in usedList:
                        usedList.append(regex)
                        tempdf = filteredTable[filteredTable.Taxon.str.match(regex)].sum(
                            numeric_only=True).to_frame().transpose()
                        tempdf.iloc[0]['Taxon'] = regex
                        tempdf.insert(loc=0, column='Taxon', value=regex)
                        indexList.append(filteredTable[filteredTable.Taxon.str.match(regex)].index)
                        if (groupedDF.empty):
                            groupedDF = tempdf
                        else:
                            groupedDF = pd.concat([groupedDF, tempdf], ignore_index=True).reset_index(drop=True)

        for i in range(0, len(usedList)):
            filteredTable.drop(index=filteredTable[filteredTable.Taxon.str.match(usedList[i])].index, inplace=True, axis=0)

        filteredTable.reset_index(drop=True, inplace=True)
        filteredTable['Taxon'] = filteredTable['Taxon'].str.rstrip(';')
        groupedDF = pd.concat([groupedDF, filteredTable], ignore_index=True).reset_index(drop=True)
        groupedDF.groupby('Taxon').sum(numeric_only=True).to_csv(
            PATH +'/winezymes_' + markerGene + '_18-19_dna-rna-grouped.tsv', sep='\t')

#visualise taxonomy of given samples
def taxonomy_analysis(markerGene):
    with open(PATH + '/winezymes_' + markerGene + '_18-19_dna-rna-grouped.tsv',) as file:
        reader = pd.read_csv(file, sep='\t')
        reader.drop(index=0, inplace=True)

        col = sns.color_palette(cc.glasbey_hv, len(reader['Taxon']))
        random.shuffle(col)
        regexList = list()
        regexList.append('18-rb-ctrl-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('18-rb-ctrl-rna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('18-rb-lasa-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('18-rb-lasa-rna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('18-rb-mesa-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('18-rb-mesa-rna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-ctrl-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-ctrl-rna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-lasa-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-lasa-rna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-mesa-dna-' + markerGene + '-m\d-rep[ABC]-\d')
        regexList.append('19-rb-mesa-rna-' + markerGene + '-m\d-rep[ABC]-\d')

        for i in reader.columns[1:]:
            reader[i] = pd.to_numeric(reader[i])

        for regex in regexList:

            x = regex.find('-m\d-')
            initialDF = reader.filter(regex=regex)
            initialDF = pd.DataFrame(initialDF.div(initialDF.sum(axis=0), axis=1))
            df = pd.DataFrame()
            if regex[1] == str(7):
                for i in range(1, 6):
                    tmpregex = regex[:x] + '-m' + str(i) + regex[x + 4:]
                    tmpDf = initialDF.filter(regex=tmpregex)
                    meanSeries = tmpDf.mean(axis=1)
                    meanSeries.name = 'm' + str(i)
                    if not math.isnan(meanSeries[1]):
                        df['m' + str(i)] = meanSeries

            else:
                for i in range(1, 4):
                    tmpregex = regex[:x] + '-m' + str(i) + regex[x + 4:]
                    tmpDf = initialDF.filter(regex=tmpregex)
                    meanSeries = tmpDf.mean(axis=1)
                    meanSeries.name = 'm' + str(i)
                    if not math.isnan(meanSeries[1]):
                        df['m' + str(i)] = meanSeries
            df = df.transpose()
            labels = list()
            if not df.empty:
                for i in range(1, len(df.iloc[0]) + 1):
                    isAbove5Percent = False
                    for j in range(0, len(df)):
                        if df.iloc[j][i] >= 0.05:
                            isAbove5Percent = True

                    if isAbove5Percent:
                        prevL = len(reader.iloc[i - 1][0])
                        l = reader.iloc[i - 1][0].rfind(';')
                        uncultured = True
                        while uncultured:
                            label = reader.iloc[i - 1][0][l + 1:prevL]
                            if label.find('uncultured') == -1:
                                uncultured = False
                            else:
                                prevL = l
                                l = reader.iloc[i - 1][0].rfind(';', 0, l)
                        if label.find('_') != -1:
                            if label.find('f__') != -1:
                                labels.append(label[label.find('_') + 2:])
                            else:
                                labels.append("$\it{" + label[label.find('_') + 2:] + "}$")
                        else:
                            labels.append("$\it{" + label + "}$")
                    else:

                        labels.append('_nolabel_', )

                df.plot(kind='bar', stacked=True, color=col, width=0.8, figsize=(5, 10)).spines[
                    ['right', 'top']].set_visible(False)

                plt.legend(labels=labels, loc='upper center', bbox_to_anchor=(0.5, -0.07), ncol=3, fontsize=15)
                plt.xticks(rotation=360, fontsize=15)
                plt.yticks(fontsize=15)
                plt.xlabel("Fázy", fontsize=15)
                plt.ylabel("Proporčné zastúpenie", fontsize=15)
                filename = make_filename(regex)
                plt.title(filename, fontsize=15)
                figure = plt.gcf()

                if not os.path.exists(PATH + "/stackplots"):
                    os.mkdir(PATH + "/stackplots")

                figure.savefig(
                   PATH + "/stackplots/" + filename + ".png", bbox_inches='tight')
                plt.clf()

#generate linechart  to see if yeast from starter survives
def starter_survival():
    with open(PATH + '/winezymes_28S_18-19_dna-rna-grouped.tsv') as file:
        reader = pd.read_csv(file, sep='\t')
        reader.drop(index=0, inplace=True)

        lasaregexList = list()
        mesaregexList = list()
        lasaregexList.append('18-rb-lasa-dna-28S-m\d-rep[ABC]-\d')
        lasaregexList.append('18-rb-lasa-rna-28S-m\d-rep[ABC]-\d')
        mesaregexList.append('18-rb-mesa-dna-28S-m\d-rep[ABC]-\d')
        mesaregexList.append('18-rb-mesa-rna-28S-m\d-rep[ABC]-\d')
        lasaregexList.append('19-rb-lasa-dna-28S-m\d-rep[ABC]-\d')
        lasaregexList.append('19-rb-lasa-rna-28S-m\d-rep[ABC]-\d')
        mesaregexList.append('19-rb-mesa-dna-28S-m\d-rep[ABC]-\d')
        mesaregexList.append('19-rb-mesa-rna-28S-m\d-rep[ABC]-\d')

        for i in reader.columns[4:]:
            reader[i] = pd.to_numeric(reader[i])
        reader.iloc[:, 4:] = pd.DataFrame(reader.iloc[:, 4:].div(reader.iloc[:, 4:].sum(axis=0), axis=1))
        sacha_regex = "Eukaryota;Fungi;Dikarya;Ascomycota;Saccharomycotina;Saccharomycetes;Saccharomycetales;Saccharomycetaceae;Saccharomyces"
        hanes_regex = "Eukaryota;Fungi;Dikarya;Ascomycota;Saccharomycotina;Saccharomycetes;Saccharomycetales;Saccharomycodaceae;Hanseniaspora"
        lasa_regex = "Eukaryota;Fungi;Dikarya;Ascomycota;Saccharomycotina;Saccharomycetes;Saccharomycetales;Saccharomycetaceae;Lachancea"
        mesa_regex = "Eukaryota;Fungi;Dikarya;Ascomycota;Saccharomycotina;Saccharomycetes;Saccharomycetales;Metschnikowiaceae;Metschnikowia"

        tmpdf = reader[reader.Taxon.str.match(sacha_regex)].sum().to_frame().transpose()
        basedf = tmpdf
        tmpdf = reader[reader.Taxon.str.match(hanes_regex)].sum().to_frame().transpose()
        basedf = basedf.append(tmpdf)
        tmpdf = reader[reader.Taxon.str.match(lasa_regex)].sum().to_frame().transpose()
        lasadf = basedf.append(tmpdf)
        tmpdf = reader[reader.Taxon.str.match(mesa_regex)].sum().to_frame().transpose()
        mesadf = basedf.append(tmpdf)
        mesadf = mesadf.reset_index().drop('index', axis=1)
        lasadf = lasadf.reset_index().drop('index', axis=1)

        for regex in lasaregexList:
            lasadffiltered = lasadf.filter(regex=regex)
            y = regex.find('-m\d-')
            anovadf = pd.DataFrame(columns=['id', 'phase', 'yeast', 'proportion'])
            for j in range(0, len(lasadffiltered.columns)):
                for i in range(0, 3):
                    x = len(lasadffiltered.columns[j])
                    id = lasadffiltered.columns[j][x - 3] + lasadffiltered.columns[j][x - 1]
                    proportion = lasadffiltered.iloc[i][j]
                    z = lasadffiltered.columns[j].find('-m')
                    phase = 'm' + lasadffiltered.columns[j][z + 2]
                    yeast = ''
                    if i == 0:
                        yeast = 'Saccharomyces'
                    if i == 1:
                        yeast = 'Hanseniaspora'
                    if i == 2:
                        yeast = 'Lachancea'
                    tmpdf = pd.DataFrame(columns=['id', 'phase', 'yeast', 'proportion'])
                    tmpdf['id'] = [id]
                    tmpdf['phase'] = [phase]
                    tmpdf['yeast'] = [yeast]
                    tmpdf['proportion'] = [proportion]
                    anovadf = anovadf.append(tmpdf, ignore_index=True)
            anovadf = anovadf.sort_values('phase')
            palette = {'Saccharomyces': 'C0', 'Hanseniaspora': 'C1', 'Lachancea': 'C2'}
            filename = make_filename(regex)
            sns.set_theme(style="whitegrid")
            g = sns.catplot(
                data=anovadf, x="phase", y="proportion", hue="yeast",
                capsize=.1, errorbar=("pi", 100), palette=palette,
                kind="point", estimator='median'
            ).set(title=filename, xlabel='Fázy', ylabel='Proporčné zastúpenie')

            g.despine(left=True)
            g._legend.set_title('Kvasinky')
            g._legend.set_title('Kvasinky')

            g._legend.texts[0].set_text('$\it{Saccharomyces}$')
            g._legend.texts[1].set_text('$\it{Hanseniaspora}$')
            g._legend.texts[2].set_text('$\it{Lachancea}$')
            figure = plt.gcf()

            if not os.path.exists(PATH + "/linechart"):
                os.mkdir(PATH + "/linechart")

            figure.savefig(
                PATH + "/linechart/" + filename + ".png", bbox_inches='tight')
            plt.clf()

        for regex in mesaregexList:
            mesadffiltered = mesadf.filter(regex=regex)
            y = regex.find('-m\d-')
            anovadf = pd.DataFrame(columns=['id', 'phase', 'yeast', 'proportion'])
            for j in range(0, len(mesadffiltered.columns)):
                for i in range(0, 3):
                    x = len(mesadffiltered.columns[j])
                    id = mesadffiltered.columns[j][x - 3] + mesadffiltered.columns[j][x - 1]
                    proportion = mesadffiltered.iloc[i][j]
                    z = mesadffiltered.columns[j].rfind('-m')
                    phase = 'm' + mesadffiltered.columns[j][z + 2]
                    yeast = ''
                    if i == 0:
                        yeast = 'Saccharomyces'
                    if i == 1:
                        yeast = 'Hanseniaspora'
                    if i == 2:
                        yeast = 'Metschnikowia'
                    tmpdf = pd.DataFrame(columns=['id', 'phase', 'yeast', 'proportion'])
                    tmpdf['id'] = [id]
                    tmpdf['phase'] = [phase]
                    tmpdf['yeast'] = [yeast]
                    tmpdf['proportion'] = [proportion]
                    anovadf = anovadf.append(tmpdf, ignore_index=True)
            anovadf = anovadf.sort_values('phase')
            palette = {'Saccharomyces': 'C0', 'Hanseniaspora': 'C1', 'Metschnikowia': 'C2'}
            filename = make_filename(regex)
            anovadf = anovadf.sort_values('phase')
            sns.set_theme(style="whitegrid")
            g = sns.catplot(
                data=anovadf, x="phase", y="proportion", hue="yeast",
                capsize=.1, errorbar=("pi", 100), palette=palette,
                kind="point", estimator='median'
            ).set(title=filename, xlabel='Fázy', ylabel='Proporčné zastúpenie')
            g._legend.set_title('Kvasinky')

            g._legend.texts[0].set_text('$\it{Saccharomyces}$')
            g._legend.texts[1].set_text('$\it{Hanseniaspora}$')
            g._legend.texts[2].set_text('$\it{Metschnikowia}$')
            g.despine(left=True)
            figure = plt.gcf()

            figure.savefig(
                PATH + "/linechart/" + filename + ".png", bbox_inches='tight')
            plt.clf()


#differential analysis using deseq2
def differential_analysis(markerGene):
    with open(PATH + '/winezymes_'+ markerGene +'_18-19_dna-rna-grouped.tsv') as file:
        reader = pd.read_csv(file, sep='\t')
        for i in reader.columns[4:]:
            reader[i] = pd.to_numeric(reader[i])
        regexList = list()
        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m1-repA-[123]', '18-rb-lasa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-lasa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m1-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-dna-' + markerGene + '-m1-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-lasa-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])

        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m1-repA-[123]', '18-rb-lasa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-lasa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m1-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-rna-' + markerGene + '-m1-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['18-rb-lasa-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])

        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m2-repA-[123]', '18-rb-lasa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-lasa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m2-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-dna-' + markerGene + '-m2-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-lasa-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])

        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m2-repA-[123]', '18-rb-lasa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-lasa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m2-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-rna-' + markerGene + '-m2-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['18-rb-lasa-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])

        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m3-repA-[123]', '18-rb-lasa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-lasa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-dna-' + markerGene + '-m3-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-ctrl-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-dna-' + markerGene + '-m3-repA-[123]', '18-rb-mesa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-lasa-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-mesa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])

        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m3-repA-[123]', '18-rb-lasa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-lasa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-ctrl-rna-' + markerGene + '-m3-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-ctrl-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['18-rb-lasa-rna-' + markerGene + '-m3-repA-[123]', '18-rb-mesa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['18-rb-lasa-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '18-rb-mesa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m1-repA-[123]', '19-rb-lasa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-lasa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m1-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-dna-' + markerGene + '-m1-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-lasa-dna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m1-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m1-repA-[123]', '19-rb-lasa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-lasa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m1-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-rna-' + markerGene + '-m1-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m1-repA-[123]'])
        regexList.append(['19-rb-lasa-rna-' + markerGene + '-m1-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m1-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m2-repA-[123]', '19-rb-lasa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-lasa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m2-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-dna-' + markerGene + '-m2-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-lasa-dna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m2-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m2-repA-[123]', '19-rb-lasa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-lasa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m2-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-rna-' + markerGene + '-m2-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m2-repA-[123]'])
        regexList.append(['19-rb-lasa-rna-' + markerGene + '-m2-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m2-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m3-repA-[123]', '19-rb-lasa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-lasa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-dna-' + markerGene + '-m3-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-ctrl-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-dna-' + markerGene + '-m3-repA-[123]', '19-rb-mesa-dna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-lasa-dna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-mesa-dna-' + markerGene + '-m3-rep[ABC]-[123]'])

        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m3-repA-[123]', '19-rb-lasa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-lasa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-ctrl-rna-' + markerGene + '-m3-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-ctrl-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])
        regexList.append(
            ['19-rb-lasa-rna-' + markerGene + '-m3-repA-[123]', '19-rb-mesa-rna-' + markerGene + '-m3-repA-[123]'])
        regexList.append(['19-rb-lasa-rna-' + markerGene + '-m3-rep[ABC]-[123]',
                          '19-rb-mesa-rna-' + markerGene + '-m3-rep[ABC]-[123]'])
        padjdf = pd.DataFrame()
        log2folddf = pd.DataFrame()
        significantdf = pd.DataFrame()
        for x in range(0, len(regexList)):
            ctrlgroup = reader.filter(regex=regexList[x][0])
            expgroup = reader.filter(regex=regexList[x][1])

            if len(ctrlgroup.columns) >= 3 and len(expgroup.columns) >= 3:
                joinedTable = pd.concat([ctrlgroup, expgroup], axis=1)
                joinedTable.insert(loc=0, column='Taxon', value=reader['Taxon'])
                joinedTable.to_csv(PATH +'/deseq2data.csv', sep=',', index=False)
                coldata = pd.DataFrame(data={'rep': [], 'additive': []})
                for i in range(1, len(joinedTable.columns)):
                    if 'ctrl' in joinedTable.columns[i]:
                        coldata = pd.concat(
                            [coldata, pd.DataFrame(data={'rep': [joinedTable.columns[i]], 'additive': ['ctrl']})])
                    if 'lasa' in joinedTable.columns[i]:
                        coldata = pd.concat(
                            [coldata, pd.DataFrame(data={'rep': [joinedTable.columns[i]], 'additive': ['lasa']})])
                    if 'mesa' in joinedTable.columns[i]:
                        coldata = pd.concat(
                            [coldata, pd.DataFrame(data={'rep': [joinedTable.columns[i]], 'additive': ['mesa']})])
                coldata.reset_index(drop=True, inplace=True)
                coldata.to_csv(PATH + '/deseq2coldata.csv', sep=',', index=False)

                pathToR = "C:/Program Files/R/R-4.2.2/bin/x64/Rscript.exe"
                pathToRScript = PATH + "/deseq2.R"
                subprocess.call([pathToR, pathToRScript, PATH], shell=True)
                with open(PATH + '/deseq2results.csv') as results:
                    resultsReader = pd.read_csv(results, sep=',')
                    year = str(regexList[x][0][0:2])
                    stage = ''
                    if 'm1' in regexList[x][0]:
                        stage = 'm1'
                    if 'm2' in regexList[x][0]:
                        stage = 'm2'
                    if 'm3' in regexList[x][0]:
                        stage = 'm3'
                    if 'dna' in regexList[x][0]:
                        material = 'dna'
                    else:
                        material = 'rna'
                    if 'repA' in regexList[x][0]:
                        replicates = 'only A'
                    else:
                        replicates = 'all'
                    if 'ctrl' in regexList[x][0]:
                        ctrl = 'ctrl'
                    else:
                        ctrl = 'lasa'
                    if 'lasa' in regexList[x][1]:
                        exp = 'lasa'
                    else:
                        exp = 'mesa'

                    baseSeries = (pd.Series(data=[year, stage, material, replicates, ctrl, exp],
                                            index=['year', 'stage', 'material', 'replicates', 'control_group',
                                                   'experiment_group']))
                    padjseries = (pd.Series(data=resultsReader['padj'].tolist(),
                                            index=resultsReader['Unnamed: 0']))
                    log2foldseries = (pd.Series(data=resultsReader['log2FoldChange'].tolist(),
                                                index=resultsReader['Unnamed: 0']))
                    significantseries = padjseries.apply(lambda x: 'x' if x < 0.01 else numpy.NaN)
                    padjseries = pd.concat([baseSeries, padjseries]).to_frame().transpose()
                    log2foldseries = pd.concat([baseSeries, log2foldseries]).to_frame().transpose()
                    significantseries = pd.concat([baseSeries, significantseries]).to_frame().transpose()
                    if padjdf.empty:
                        padjdf = padjseries
                    else:
                        padjdf = pd.concat([padjdf, padjseries])
                    if log2folddf.empty:
                        log2folddf = log2foldseries
                    else:
                        log2folddf = pd.concat([log2folddf, log2foldseries])
                    if significantdf.empty:
                        significantdf = significantseries
                    else:
                        significantdf = pd.concat([significantdf, significantseries])
                    results.close()
        log2folddf.reset_index(drop=True, inplace=True)
        padjdf.reset_index(drop=True, inplace=True)
        significantdf.reset_index(drop=True, inplace=True)
        sortcols = significantdf.iloc[:, 6:]
        sortcols.fillna(0)
        sortcols = sortcols.mask(sortcols == 'x', 1)
        for i in sortcols.columns:
            sortcols[:][i] = pd.to_numeric(sortcols[:][i])
        sortedcols = sorted(sortcols.columns, key=lambda x: sortcols[x].sum(), reverse=True)
        sortedsignificantdf = significantdf[sortedcols]
        sortedlog2folddf = log2folddf[sortedcols]
        sortedpadjdf = padjdf[sortedcols]
        significantdf = pd.concat([significantdf.iloc[:, :6], sortedsignificantdf], axis=1)
        log2folddf = pd.concat([significantdf.iloc[:, :6], sortedlog2folddf], axis=1)
        padjdf = pd.concat([significantdf.iloc[:, :6], sortedpadjdf], axis=1)
        for i in range(0, len(significantdf.columns)):
            pos = significantdf.columns[i].rfind(';')
            significantdf = significantdf.rename(columns={significantdf.columns[i]: significantdf.columns[i][pos + 1:]})
            log2folddf = log2folddf.rename(columns={log2folddf.columns[i]: log2folddf.columns[i][pos + 1:]})
            padjdf = padjdf.rename(columns={padjdf.columns[i]: padjdf.columns[i][pos + 1:]})
        log2folddf.to_csv(PATH + '/deseq2-log2fold-' + markerGene +'.csv')
        padjdf.to_csv(PATH + '/deseq2-padj-' + markerGene +'.csv')
        significantdf.to_csv(PATH + '/deseq2-significant-' + markerGene +'.csv')

#rarefaction of alpha diversity
def rarefaction(markerGene):
    with open(PATH + '/winezymes_'+ markerGene +'_18-19_dna-rna-grouped.tsv') as file:
        reader = pd.read_csv(file, sep='\t')
        reader.reset_index(drop=True, inplace=True)
        regexList = ['18-rb-ctrl-dna-16S-m1-rep[ABC]-[123]', '18-rb-lasa-dna-16S-m1-rep[ABC]-[123]',
                     '18-rb-mesa-dna-16S-m1-rep[ABC]-[123]',
                     '18-rb-ctrl-rna-16S-m1-rep[ABC]-[123]', '18-rb-lasa-rna-16S-m1-rep[ABC]-[123]',
                     '18-rb-mesa-rna-16S-m1-rep[ABC]-[123]',
                     '19-rb-ctrl-dna-16S-m1-rep[ABC]-[123]', '19-rb-lasa-dna-16S-m1-rep[ABC]-[123]',
                     '19-rb-mesa-dna-16S-m1-rep[ABC]-[123]',
                     '19-rb-ctrl-rna-16S-m1-rep[ABC]-[123]', '19-rb-lasa-rna-16S-m1-rep[ABC]-[123]',
                     '19-rb-mesa-rna-16S-m1-rep[ABC]-[123]',
                     '18-rb-ctrl-dna-16S-m2-rep[ABC]-[123]', '18-rb-lasa-dna-16S-m2-rep[ABC]-[123]',
                     '18-rb-mesa-dna-16S-m2-rep[ABC]-[123]',
                     '18-rb-ctrl-rna-16S-m2-rep[ABC]-[123]', '18-rb-lasa-rna-16S-m2-rep[ABC]-[123]',
                     '18-rb-mesa-rna-16S-m2-rep[ABC]-[123]',
                     '19-rb-ctrl-dna-16S-m2-rep[ABC]-[123]', '19-rb-lasa-dna-16S-m2-rep[ABC]-[123]',
                     '19-rb-mesa-dna-16S-m2-rep[ABC]-[123]',
                     '19-rb-ctrl-rna-16S-m2-rep[ABC]-[123]', '19-rb-lasa-rna-16S-m2-rep[ABC]-[123]',
                     '19-rb-mesa-rna-16S-m2-rep[ABC]-[123]',
                     '18-rb-ctrl-dna-16S-m3-rep[ABC]-[123]', '18-rb-lasa-dna-16S-m3-rep[ABC]-[123]',
                     '18-rb-mesa-dna-16S-m3-rep[ABC]-[123]',
                     '18-rb-ctrl-rna-16S-m3-rep[ABC]-[123]', '18-rb-lasa-rna-16S-m3-rep[ABC]-[123]',
                     '18-rb-mesa-rna-16S-m3-rep[ABC]-[123]',
                     '19-rb-ctrl-dna-16S-m3-rep[ABC]-[123]', '19-rb-lasa-dna-16S-m3-rep[ABC]-[123]',
                     '19-rb-mesa-dna-16S-m3-rep[ABC]-[123]',
                     '19-rb-ctrl-rna-16S-m3-rep[ABC]-[123]', '19-rb-lasa-rna-16S-m3-rep[ABC]-[123]',
                     '19-rb-mesa-rna-16S-m3-rep[ABC]-[123]']
        for i in reader.columns[1:]:
            reader[i] = pd.to_numeric(reader[i])
        finalfinalDataframe = pd.DataFrame()
        count = 0
        for regex in regexList:
            count += 1
            n = 10
            listOfSamples = list()
            filteredDF = reader.filter(regex=regex)
            for i in range(0, len(filteredDF.columns)):
                sampleList = list()
                dataList = list()
                for j in range(len(filteredDF)):
                    sampleList.extend([j] * int(filteredDF.iloc[j][i]))
                for x in range(n):
                    random.shuffle(sampleList)
                    listOfSamples.append(sampleList)
                xSeries = list()
                listOfTempLists = [[0] * len(filteredDF)] * n
                for g in range(0, len(sampleList), 100):
                    dataFrame = pd.DataFrame()

                    for z in range(n):
                        if i != 0:
                            randomList = listOfSamples[z][g - 100:g]
                        else:
                            randomList = listOfSamples[z][:g]
                        for a in randomList:
                            listOfTempLists[z][a] += 1
                        observed = skbio.diversity.alpha_diversity(metric="observed_otus", counts=listOfTempLists)
                        tempSeries = pd.Series(observed)
                        if dataFrame.empty:
                            dataFrame = pd.DataFrame({g: tempSeries})
                        else:
                            dataFrame = pd.concat([dataFrame, tempSeries], axis=1)

                    y = dataFrame.mean(axis=1).tolist()
                    dataList.append(y[0])
                    xSeries.extend([g])

                xSeries = pd.Series(xSeries)
                finalSeries = pd.Series(dataList)
                finalDataFrame = pd.DataFrame({'number of samples': xSeries, filteredDF.columns[i]: finalSeries})
                if finalfinalDataframe.empty:
                    finalfinalDataframe = finalDataFrame
                else:
                    finalfinalDataframe = finalfinalDataframe.merge(finalDataFrame, on='number of samples', how='outer')

            if regex.find('mesa') != -1:
                if not finalfinalDataframe.empty:
                    finalfinalDataframe = finalfinalDataframe.set_index('number of samples')
                    if regex.find('m1') != -1:
                        phase = 'm1'
                    if regex.find('m2') != -1:
                        phase = 'm2'
                    if regex.find('m3') != -1:
                        phase = 'm3'
                    title = make_filename(regex)
                    title = title[:title.rfind('mesa')] + " " + phase
                    c = {}
                    for col in finalfinalDataframe.columns:

                        if col.find('ctrl') != -1:
                            c[col] = 'r'
                        if col.find('mesa') != -1:
                            c[col] = 'b'
                        if col.find('lasa') != -1:
                            c[col] = 'g'

                    figure, ax = plt.subplots()
                    sns.lineplot(data=finalfinalDataframe, palette=c, ax=ax)
                    ax.set_title(title, fontsize=13)
                    ax.set_xlabel('Počet čítaní', fontsize=13)
                    ax.set_ylabel('Bohatosť', fontsize=13)
                    ax.set_xticklabels(ax.get_xticklabels(), fontsize=12)
                    ax.set_yticklabels(ax.get_yticklabels(), fontsize=12)
                    ax.spines[
                        ['right', 'top']].set_visible(False)
                    ax.grid(False)
                    ax.legend(title='Vzorky', fontsize=12)
                    sns.move_legend(ax, "upper left", bbox_to_anchor=(1, 1))
                    plt.setp(ax.get_legend().get_title(), fontsize='13')

                    if not os.path.exists(PATH + "/rarefaction"):
                        os.mkdir(PATH + "/rarefaction")

                    figure.savefig(
                        PATH + "/rarefaction/" + title + ".png", bbox_inches='tight')
                    finalfinalDataframe = pd.DataFrame()

#comparison of alpha diversity between groups of samples
def alpha_diversity_comparison(markerGene):
    with open(PATH + '/winezymes_' + markerGene + '_18-19_dna-rna-grouped.tsv') as file:
        reader = pd.read_csv(file, sep='\t')
        reader.reset_index(drop=True, inplace=True)
        tips = sns.load_dataset("tips")
        df = pd.DataFrame(columns=['year', 'phase', 'material', 'additive', 'observed_asvs', 'shannon'])

        for col in reader.loc[:, reader.columns != 'Taxon']:
            observed = pd.Series(skbio.diversity.alpha_diversity(metric="observed_otus", counts=reader[col]))[0]
            shannon = pd.Series(skbio.diversity.alpha_diversity(metric="shannon", counts=reader[col]))[0]
            year = col[:2]
            if col.find('dna') != -1:
                material = 'dna'
            else:
                material = 'rna'
            if col.find('m1') != -1:
                phase = 'm1'
            if col.find('m2') != -1:
                phase = 'm2'
            if col.find('m3') != -1:
                phase = 'm3'

            if col.find('ctrl') != -1:
                additive = 'ctrl'
            if col.find('lasa') != -1:
                additive = 'lasa'
            if col.find('mesa') != -1:
                additive = 'mesa'

            df = pd.concat([df, pd.DataFrame([{'year': year, 'phase': phase, 'material': material, 'additive': additive,
                                               'observed_asvs': observed, 'shannon': shannon}])])
            df.reset_index(drop=True, inplace=True)

        figure, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(15, 8))
        sns.set_context('paper')
        sns.boxplot(x="additive", y="observed_asvs",
                    hue="phase",
                    data=df.query("material == 'dna' and year == '18'"), hue_order=['m1', 'm2', 'm3'], ax=axes[0][0])
        axes[0][0].set_title('2018 DNA ' + markerGene + ' bohatosť', fontsize=14)
        axes[0][0].set_xlabel('')
        axes[0][0].set_ylabel('Bohatosť', fontsize=13)
        axes[0][0].spines[['right', 'top']].set_visible(False)
        axes[0][0].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][0].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][0].legend([], [], frameon=False)

        sns.boxplot(x="additive", y="observed_asvs",
                    hue="phase",
                    data=df.query("material == 'dna' and year == '19'"), hue_order=['m1', 'm2', 'm3'], ax=axes[1][0])
        axes[1][0].set_title('2019 DNA ' + markerGene + ' bohatosť', fontsize=14)
        axes[1][0].set_xlabel('Štartér', fontsize=13)
        axes[1][0].set_ylabel('Bohatosť', fontsize=13)
        axes[1][0].spines[['right', 'top']].set_visible(False)
        axes[1][0].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][0].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][0].legend([], [], frameon=False)
        sns.boxplot(x="additive", y="observed_asvs",
                    hue="phase",
                    data=df.query("material == 'rna' and year == '18'"), hue_order=['m1', 'm2', 'm3'], ax=axes[0][1])
        axes[0][1].set_title('2018 RNA ' + markerGene + ' bohatosť', fontsize=14)
        axes[0][1].set_xlabel('')
        axes[0][1].set_ylabel('')
        axes[0][1].spines[['right', 'top']].set_visible(False)
        axes[0][1].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][1].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][1].legend([], [], frameon=False)
        axes[0][1].legend(bbox_to_anchor=(1.02, 1), loc='upper left', borderaxespad=0, title='Fáza', fontsize=12)
        sns.boxplot(x="additive", y="observed_asvs",
                    hue="phase",
                    data=df.query("material == 'rna' and year == '19'"), hue_order=['m1', 'm2', 'm3'], ax=axes[1][1])
        axes[1][1].set_title('2019 RNA ' + markerGene + ' bohatosť', fontsize=14)
        axes[1][1].set_xlabel('Štartér', fontsize=13)
        axes[1][1].set_ylabel('')
        axes[1][1].spines[['right', 'top']].set_visible(False)
        axes[1][1].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][1].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][1].legend([], [], frameon=False)
        mng = plt.get_current_fig_manager()
        mng.resize(*mng.window.maxsize())
        figure = plt.gcf()
        plt.setp(axes[0][1].get_legend().get_title(), fontsize='14')

        if not os.path.exists(PATH + "/barplots"):
            os.mkdir(PATH + "/barplots")

        figure.savefig(
            PATH + '/barplots/' + markerGene + '-richness-phases.png',
            bbox_inches='tight')

        figure, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(15, 8))
        sns.set_context('paper')
        sns.boxplot(x="additive", y="shannon",
                    hue="phase",
                    data=df.query("material == 'dna' and year == '18'"), hue_order=['m1', 'm2', 'm3'], ax=axes[0][0])
        axes[0][0].set_title('2018 DNA ' + markerGene + ' rovnováha', fontsize=14)
        axes[0][0].set_xlabel('')
        axes[0][0].set_ylabel('Rovnováha', fontsize=13)
        axes[0][0].spines[['right', 'top']].set_visible(False)
        axes[0][0].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][0].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][0].legend([], [], frameon=False)
        sns.boxplot(x="additive", y="shannon",
                    hue="phase",
                    data=df.query("material == 'dna' and year == '19'"), hue_order=['m1', 'm2', 'm3'], ax=axes[1][0])
        axes[1][0].set_title('2019 DNA ' + markerGene + ' rovnováha', fontsize=14)
        axes[1][0].set_xlabel('Štartér', fontsize=13)
        axes[1][0].set_ylabel('Rovnováha', fontsize=13)
        axes[1][0].spines[['right', 'top']].set_visible(False)
        axes[1][0].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][0].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][0].legend([], [], frameon=False)
        sns.boxplot(x="additive", y="shannon",
                    hue="phase",
                    data=df.query("material == 'rna' and year == '18'"), hue_order=['m1', 'm2', 'm3'], ax=axes[0][1])
        axes[0][1].set_title('2018 RNA ' + markerGene + ' rovnováha', fontsize=14)
        axes[0][1].set_xlabel('')
        axes[0][1].set_ylabel('')
        axes[0][1].spines[['right', 'top']].set_visible(False)
        axes[0][1].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][1].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[0][1].legend([], [], frameon=False)
        axes[0][1].legend(bbox_to_anchor=(1.02, 1), loc='upper left', borderaxespad=0, title='Fáza', fontsize=12)
        sns.boxplot(x="additive", y="shannon",
                    hue="phase",
                    data=df.query("material == 'rna' and year == '19'"), hue_order=['m1', 'm2', 'm3'], ax=axes[1][1])
        axes[1][1].set_title('2019 RNA ' + markerGene + ' rovnováha', fontsize=14)
        axes[1][1].set_xlabel('Štartér', fontsize=13)
        axes[1][1].set_ylabel('')
        axes[1][1].spines[['right', 'top']].set_visible(False)
        axes[1][1].xaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][1].yaxis.set_tick_params(labelbottom=True, labelsize=13)
        axes[1][1].legend([], [], frameon=False)

        mng = plt.get_current_fig_manager()
        mng.resize(*mng.window.maxsize())
        figure = plt.gcf()
        plt.setp(axes[0][1].get_legend().get_title(), fontsize='14')

        if not os.path.exists(PATH + "/barplots"):
            os.mkdir(PATH + "/barplots")

        figure.savefig(
            PATH + '/barplots/' + markerGene + '-evennes-phases.png',
            bbox_inches='tight')
        plt.clf()
#PCoA of beta diversities
def PCoA_beta_diversity(markerGene):
    with open(PATH + '/winezymes_' + markerGene + '_18-19_dna-rna.tsv') as file:
        reader = pd.read_csv(file, sep='\t')
        reader.drop(index=0, inplace=True)
        for i in reader.columns[3:]:
            reader[i] = pd.to_numeric(reader[i])
        filterReader = reader.filter(
            regex='1[89]-rb-[cml][eat][sr][al]-[rd]na-' + markerGene + '-m[123]-rep[ABC]-[123]')
        filterReader = filterReader / filterReader.sum()
        transposedDna = filterReader.filter(
            regex='1[89]-rb-[cml][eat][sr][al]-dna-' + markerGene + '-m[123]-rep[ABC]-[123]').transpose()
        transposedRna = filterReader.filter(
            regex='1[89]-rb-[cml][eat][sr][al]-rna-' + markerGene + '-m[123]-rep[ABC]-[123]').transpose()

        for i in transposedDna.columns:
            transposedDna[i] = pd.to_numeric(transposedDna[i])
        for i in transposedRna.columns:
            transposedRna[i] = pd.to_numeric(transposedRna[i])

        pca = PCA(n_components=2).set_output(transform='pandas')
        dnaPca = pca.fit_transform(transposedDna)
        rnaPca = pca.fit_transform(transposedRna)
        dnaPca.reset_index(inplace=True)
        rnaPca.reset_index(inplace=True)
        dnaPca[['year', 'wine', 'additive', 'material', 'marker', 'phase', 'biorep', 'techrep']] = dnaPca[
            'index'].str.split('-', expand=True)
        rnaPca[['year', 'wine', 'additive', 'material', 'marker', 'phase', 'biorep', 'techrep']] = rnaPca[
            'index'].str.split('-', expand=True)
        dnaPca.set_index('index', inplace=True)
        rnaPca.set_index('index', inplace=True)

        beta_diversity_dna = skbio.diversity.beta_diversity(counts=transposedDna, ids=transposedDna.index,
                                                            metric="braycurtis")
        beta_diversity_rna = skbio.diversity.beta_diversity(counts=transposedRna, ids=transposedRna.index,
                                                            metric="braycurtis")

        dnaPca = skbio.stats.ordination.pcoa(distance_matrix=beta_diversity_dna, number_of_dimensions=3).samples
        rnaPca = skbio.stats.ordination.pcoa(distance_matrix=beta_diversity_rna, number_of_dimensions=3).samples
        tmp_dna = pd.DataFrame()
        for i in dnaPca.index:
            tmp_dna = pd.concat([tmp_dna, pd.Series(i.split('-'), index=['Rok', 'wine', 'Štartér', 'material',
                                                                         'marker', 'Fáza', 'biorep', 'techrep'])],
                                axis=1)
        tmp_dna = tmp_dna.transpose().set_index(dnaPca.index)
        dnaPca = pd.concat([dnaPca, tmp_dna], axis=1)

        tmp_rna = pd.DataFrame()
        for i in rnaPca.index:
            tmp_rna = pd.concat([tmp_rna, pd.Series(i.split('-'), index=['Rok', 'wine', 'Štartér', 'material',
                                                                         'marker', 'Fáza', 'biorep', 'techrep'])],
                                axis=1)
        tmp_rna = tmp_rna.transpose().set_index(rnaPca.index)
        rnaPca = pd.concat([rnaPca, tmp_rna], axis=1)
        cmap = {'18': 'b', '19': 'r'}
        size = {'m1': 15, 'm2': 45, 'm3': 75}
        marker = {'ctrl': 'o', 'lasa': 'X', 'mesa': '^'}

        figure, axes = plt.subplots(nrows=2, ncols=2, figsize=(11, 11))
        figure.suptitle('PCoA pre ' + markerGene + ' DNA')
        sns.scatterplot(data=dnaPca, x='PC1', y='PC2', hue='Rok', size='Fáza', style='Štartér', ax=axes[0][0]
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[0][0].set_box_aspect(1)
        axes[0][0].set(title='PC1 a PC2', xticklabels=[], yticklabels=[])
        sns.scatterplot(data=dnaPca, x='PC2', y='PC3', hue='Rok', size='Fáza', style='Štartér', ax=axes[1][1],
                        legend=False
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[1][1].set_box_aspect(1)
        axes[1][1].set(title='PC2 a PC3', xticklabels=[], yticklabels=[])
        sns.scatterplot(data=dnaPca, x='PC1', y='PC3', hue='Rok', size='Fáza', style='Štartér', ax=axes[1][0],
                        legend=False
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[1][0].set_box_aspect(1)
        axes[1][0].set(title='PC1 a PC3', xticklabels=[], yticklabels=[])

        axes[0][0].legend(bbox_to_anchor=(1.5, 0.9), loc='upper left', borderaxespad=0, fontsize=13)
        label_params = axes[0][0].get_legend_handles_labels()
        axes[0][1].clear()
        axes[0][1].set_axis_off()

        figure = plt.gcf()
        if not os.path.exists(PATH + "/pca"):
            os.mkdir(PATH + "/pca")
        figure.savefig(
            PATH + '/pca/2D-pcoa-' + markerGene + '-dna.png', bbox_inches='tight')

        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')

        for i in dnaPca.index:
            ax.scatter(dnaPca.loc[i]['PC1'], dnaPca.loc[i]['PC2'], dnaPca.loc[i]['PC3'],
                       c=cmap.get(dnaPca.loc[i]['Rok']),
                       marker=marker.get(dnaPca.loc[i]['Štartér']), s=size.get(dnaPca.loc[i]['Fáza']), alpha=0.5)

        ax.set(title='3D projekcia PCoA pre ' + markerGene + ' DNA', xticklabels=[], yticklabels=[], zticklabels=[],
               xlabel='PC1',
               ylabel='PC2', zlabel='PC3')
        ax.legend(*label_params, bbox_to_anchor=(1.02, 1), loc='upper left', borderaxespad=0)
        figure = plt.gcf()
        figure.savefig(
            PATH + '/pca/3D-pcoa-' + markerGene + '-dna.png', bbox_inches='tight')

        figure, axes = plt.subplots(nrows=2, ncols=2, figsize=(11, 11))
        figure.suptitle('PCoA pre ' + markerGene + ' RNA')
        sns.scatterplot(data=rnaPca, x='PC1', y='PC2', hue='Rok', size='Fáza', style='Štartér', ax=axes[0][0]
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[0][0].set_box_aspect(1)
        axes[0][0].set(title='PC1 a PC2', xticklabels=[], yticklabels=[])
        sns.scatterplot(data=rnaPca, x='PC2', y='PC3', hue='Rok', size='Fáza', style='Štartér', ax=axes[1][1],
                        legend=False
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[1][1].set_box_aspect(1)
        axes[1][1].set(title='PC2 a PC3', xticklabels=[], yticklabels=[])
        sns.scatterplot(data=rnaPca, x='PC1', y='PC3', hue='Rok', size='Fáza', style='Štartér', ax=axes[1][0],
                        legend=False
                        , palette=cmap, sizes=size, markers=marker, alpha=0.5)
        axes[1][0].set_box_aspect(1)
        axes[1][0].set(title='PC1 a PC3', xticklabels=[], yticklabels=[])

        axes[0][0].legend(bbox_to_anchor=(1.5, 0.9), loc='upper left', borderaxespad=0, fontsize=13)
        label_params = axes[0][0].get_legend_handles_labels()
        axes[0][1].clear()
        axes[0][1].set_axis_off()

        figure = plt.gcf()
        figure.savefig(
            PATH + '/pca/2D-pcoa-' + markerGene + '-rna.png', bbox_inches='tight')

        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')

        for i in rnaPca.index:
            ax.scatter(rnaPca.loc[i]['PC1'], rnaPca.loc[i]['PC2'], rnaPca.loc[i]['PC3'],
                       c=cmap.get(rnaPca.loc[i]['Rok']),
                       marker=marker.get(rnaPca.loc[i]['Štartér']), s=size.get(rnaPca.loc[i]['Fáza']), alpha=0.5)

        ax.set(title='3D projekcia PCoA pre ' + markerGene + ' RNA', xticklabels=[], yticklabels=[], zticklabels=[],
               xlabel='PC1',
               ylabel='PC2', zlabel='PC3')
        ax.legend(*label_params, bbox_to_anchor=(1.02, 1), loc='upper left', borderaxespad=0)
        figure = plt.gcf()
        figure.savefig(
            PATH + '/pca/3D-pcoa-' + markerGene + '-rna.png', bbox_inches='tight')
        plt.clf()


#CALLS TO THE FUCNTIONS

#group_ASVs_by_taxon(PATH + '/winezymes_16S_18-19_dna-rna.tsv', '16S')
#group_ASVs_by_taxon(PATH + '/winezymes_28S_18-19_dna-rna.tsv', '28S')

#taxonomy_analysis('16S')
#taxonomy_analysis('28S')

#starter_survival()

#differential_analysis('16S')
#differential_analysis('28S')

#rarefaction('16S')
#rarefaciton('28S')

#alpha_diversity_comparison('16S')
#alpha_diversity_comparison('28S')

#PCoA_beta_diversity('16S')
#PCoA_beta_diversity('28S')

