<jupyter_start><jupyter_text>Protein Feature Figures  10K Contigs<jupyter_code>bac = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Bac10k_AA_Features.csv')
viral = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Phage10k_AA_Features.csv')
viral["Type"] = "Viral"
bac["Type"] = "Bacterial"
df_combined = pd.concat([viral, bac])

fig = plt.figure(figsize = (20, 12))

ax1 = fig.add_subplot(221) #top left
ax2 = fig.add_subplot(222) #top right
ax3 = fig.add_subplot(223) #bottom left
ax4 = fig.add_subplot(224) #bottom right

ax1.annotate("A", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax2.annotate("B", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax3.annotate("C", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax4.annotate("D", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")



ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Gene_Density'], palette="Set3", showfliers=False, ax=ax1)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Gene Density (Genes per 1K Bases)", fontsize=16)

ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Median Operon Length'], palette="Set3", showfliers=False, ax=ax2)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Median Operon Length", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Percent_Overlap'], palette="Set3", showfliers=False, ax=ax3)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)

ax.set_ylabel("Percentage of Genes Overlapping", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['AA_Length'], palette="Set3", showfliers=False, ax=ax4)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)


ax.set_ylabel("Median Amino Acid Length", fontsize=16)

fig.savefig("../figures/Figure_1_VirBrant_Protein_Features/Figure_1_10k.jpg")
print("Gene Density")
print(kruskal(viral['Gene_Density'], bac['Gene_Density']))
print("--------")
print("Median Operon Length")
print(kruskal(viral['Median Operon Length'], bac['Median Operon Length']))
print("--------")
print("Percentage of Genes Overlapping")
print(kruskal(viral['Percent_Overlap'], bac['Percent_Overlap']))
print("--------")
print("Median Amino Acid Length")
print(kruskal(viral['AA_Length'], bac['AA_Length']))<jupyter_output>Gene Density
KruskalResult(statistic=18799.040839515485, pvalue=0.0)
--------
Median Operon Length
KruskalResult(statistic=10443.506185165259, pvalue=0.0)
--------
Percentage of Genes Overlapping
KruskalResult(statistic=9105.193653235168, pvalue=0.0)
--------
Median Amino Acid Length
KruskalResult(statistic=15255.432738095467, pvalue=0.0)<jupyter_text>5K Contigs<jupyter_code>bac = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Bac5k_AA_Features.csv')
viral = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Phage5k_AA_Features.csv')
viral["Type"] = "Viral"
bac["Type"] = "Bacterial"
df_combined = pd.concat([viral, bac])

fig = plt.figure(figsize = (20, 12))

ax1 = fig.add_subplot(221) #top left
ax2 = fig.add_subplot(222) #top right
ax3 = fig.add_subplot(223) #bottom left
ax4 = fig.add_subplot(224) #bottom right

ax1.annotate("A", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax2.annotate("B", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax3.annotate("C", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax4.annotate("D", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")



ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Gene_Density'], palette="Set3", showfliers=False, ax=ax1)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Gene Density (Genes per 1K Bases)", fontsize=16)

ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Median Operon Length'], palette="Set3", showfliers=False, ax=ax2)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Median Operon Length", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Percent_Overlap'], palette="Set3", showfliers=False, ax=ax3)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)

ax.set_ylabel("Percentage of Genes Overlapping", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['AA_Length'], palette="Set3", showfliers=False, ax=ax4)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)


ax.set_ylabel("Median Amino Acid Length", fontsize=16)

fig.savefig("../figures/Figure_1_VirBrant_Protein_Features/Figure_Supplemental_5k.jpg")
print("Gene Density")
print(kruskal(viral['Gene_Density'], bac['Gene_Density']))
print("--------")
print("Median Operon Length")
print(kruskal(viral['Median Operon Length'], bac['Median Operon Length']))
print("--------")
print("Percentage of Genes Overlapping")
print(kruskal(viral['Percent_Overlap'], bac['Percent_Overlap']))
print("--------")
print("Median Amino Acid Length")
print(kruskal(viral['AA_Length'], bac['AA_Length']))<jupyter_output>Gene Density
KruskalResult(statistic=29679.77969458356, pvalue=0.0)
--------
Median Operon Length
KruskalResult(statistic=14026.884842893342, pvalue=0.0)
--------
Percentage of Genes Overlapping
KruskalResult(statistic=12530.12264168525, pvalue=0.0)
--------
Median Amino Acid Length
KruskalResult(statistic=18611.872533397534, pvalue=0.0)<jupyter_text>3K Contigs<jupyter_code>bac = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Bac3k_AA_Features.csv')
viral = pd.read_csv('../data/VirBrant_Training_Data/Protein_Features/Phage3k_AA_Features.csv')
viral["Type"] = "Viral"
bac["Type"] = "Bacterial"
df_combined = pd.concat([viral, bac])

fig = plt.figure(figsize = (20, 12))

ax1 = fig.add_subplot(221) #top left
ax2 = fig.add_subplot(222) #top right
ax3 = fig.add_subplot(223) #bottom left
ax4 = fig.add_subplot(224) #bottom right

ax1.annotate("A", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax2.annotate("B", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax3.annotate("C", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")
ax4.annotate("D", xy=(-0.1, 1.05), xycoords="axes fraction", fontsize=20, weight="bold")



ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Gene_Density'], palette="Set3", showfliers=False, ax=ax1)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Gene Density (Genes per 1K Bases)", fontsize=16)

ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Median Operon Length'], palette="Set3", showfliers=False, ax=ax2)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)
ax.set_ylabel("Median Operon Length", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['Percent_Overlap'], palette="Set3", showfliers=False, ax=ax3)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)

ax.set_ylabel("Percentage of Genes Overlapping", fontsize=16)


ax = sns.boxplot(x=df_combined["Type"], y=df_combined['AA_Length'], palette="Set3", showfliers=False, ax=ax4)
ax.set_title("")
ax.set_xlabel("")
ax.xaxis.set_tick_params(labelsize=20)
ax.yaxis.set_tick_params(labelsize=14)


ax.set_ylabel("Median Amino Acid Length", fontsize=16)

fig.savefig("../figures/Figure_1_VirBrant_Protein_Features/Figure_Supplemental_3k.jpg")
print("Gene Density")
print(kruskal(viral['Gene_Density'], bac['Gene_Density']))
print("--------")
print("Median Operon Length")
print(kruskal(viral['Median Operon Length'], bac['Median Operon Length']))
print("--------")
print("Percentage of Genes Overlapping")
print(kruskal(viral['Percent_Overlap'], bac['Percent_Overlap']))
print("--------")
print("Median Amino Acid Length")
print(kruskal(viral['AA_Length'], bac['AA_Length']))<jupyter_output>Gene Density
KruskalResult(statistic=43416.64351025498, pvalue=0.0)
--------
Median Operon Length
KruskalResult(statistic=31201.90289236799, pvalue=0.0)
--------
Percentage of Genes Overlapping
KruskalResult(statistic=23789.984206461988, pvalue=0.0)
--------
Median Amino Acid Length
KruskalResult(statistic=20591.060407121047, pvalue=0.0)<jupyter_text>Methods to Aquire Protein Features (Built into VirBrant)** Methods detailed below does not include files (due to size)   Get ORFs using Prodigal<jupyter_code>%%bash
prodigal -i [input_fasta_file] -o genes_output.txt -a Amino_Acids.faa -p meta<jupyter_output><empty_output><jupyter_text>Get Headers and Sizes of Amino Acids and Contigs<jupyter_code>%%bash
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' [input_fasta_file] > Sizes_of_input_fasta_file.txt


awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' Amino_Acids.faa > Sizes_of_Amino_Acids.txt<jupyter_output><empty_output><jupyter_text>Get Operon Features<jupyter_code>f = open("Sizes_of_input_fasta_file.txt", "r")
lengths = f.read().split(">")[1:]
df = pd.DataFrame(lengths)
df['ID'] = [x.split("\n")[0] for x in df[0]]
df['Size'] = [x.split("\n")[1] for x in df[0]]
df.head()
df_lengths = df[["ID", "Size"]]

import pandas as pd
f = open("Sizes_of_Amino_Acids.txt", "r")
lengths = f.read().split(">")[1:]
df = pd.DataFrame(lengths)
df['ID'] = [x.split("\n")[0] for x in df[0]]
df['Size'] = [x.split("\n")[1] for x in df[0]]
df.head()
df = df[["ID", "Size"]]

Gene_Count = len(df)
df['Genome'] = ["_".join(x.split('_')[0:3]).lstrip() for x in df["ID"]]
df['Strand'] = [pd.to_numeric(x.split('#')[3].rstrip().lstrip()) for x in df["ID"]]
df['GC'] = [pd.to_numeric(x.split('#')[4].split(';')[5].split("=")[1]) for x in df["ID"]]
df['Start'] = [pd.to_numeric(x.split('#')[1].rstrip().lstrip()) for x in df["ID"]]
df['Stop'] = [pd.to_numeric(x.split('#')[2].rstrip().lstrip()) for x in df["ID"]]
df.drop("ID", axis=1, inplace=True)

Groups_df = [v for k, v in df.groupby('Genome', as_index=False)]
def Operon_Distribution(list_of_strains):
    """
    Input: Pandas Strand Column 
    Output: Number of Operons, mean size of operons, stdev of operons, median of operons
    Usage: 
    
    """
    from statistics import median, mean, stdev
    gene_strand_list = []
    strands = list(list_of_strains)
    counts = 1
    if len(strands)==1:
        return(len(strands), len(strands), len(strands), len(strands))
    for i in range(0, len(strands)-1):
        if strands[i]==strands[i+1]:
            counts+=1
        
        ## Deal with End Case
            if (i+2)==len(strands):
                gene_strand_list.append(counts)
        else:
            gene_strand_list.append(counts)
            counts = 1
    return(len(gene_strand_list), max(gene_strand_list), mean(gene_strand_list), median(gene_strand_list))
## This chunk creates a dictionary containing genome strand absolute values
genome_strands = {}
for file in range(len(Groups_df)):
    df2 = pd.DataFrame(Groups_df[file])
    Gene_Count = len(df2)

    ## Get GC Percentage
    GC = df2["GC"].median()
    
    AA_Length = df2["Size"].median()
    
    ## Get Overlap Percentage
    Overlap = 0
    if len(df2)==1:
        Percent_Overlap=0
    else:
        for i in range(len(df2)-1):
            if df2.iloc[i]['Strand'] != df2.iloc[i+1]['Strand']:
                next
            else:
                if df2.iloc[i]['Stop'] > df2.iloc[i+1]['Start']:
                    Overlap += 1
        Percent_Overlap = Overlap/len(df2)  
    
    
    ids = list(df2["Genome"])[0]
    length = df_lengths[df_lengths['ID'].str.match(ids)]

    if len(length)==0:
        print("Empty: ")
        continue
    lengths = int(length.iloc[0][1])
    Gene_Density = float(Gene_Count)/(lengths/1000)
    num_operons, max_operon, mean_operon_length, median_operon_length = Operon_Distribution(df2["Strand"])
    
    genome_strands[ids]=[Gene_Density, median_operon_length, Percent_Overlap, AA_Length]
    
df_features = pd.DataFrame(genome_strands)
df_features = df_features.transpose()
df_features.columns = ["Gene_Density", "Median Operon Length", "Percent_Overlap", "AA_Length"]
Groups_df=1
df_features.to_csv('Input_Fasta_File_AA_Features.csv')
### End of Script<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1825309706758664, 'fraction_numerical': 0.043750980084679314, 'mean_word_length': 5.950953678474114, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '473054', 'license': 'MIT', 'n_tokens_mistral': 5400, 'n_tokens_neox': 4950, 'n_words': 820, 'path': 'scripts/Figure_1_VirBrant_Protein_Features.ipynb', 'repo_name': 'Strong-Lab/Viral_Classification_in_Metagenomics', 'size': '12754'}