<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Get information from GFF file1.1&nbsp;&nbsp;Convert GFF to Pandas DataFrame1.2&nbsp;&nbsp;Check genes1.3&nbsp;&nbsp;Generate nucleotide fasta files for CDS2&nbsp;&nbsp;Get COGs2.1&nbsp;&nbsp;Run EggNOG Mapper2.2&nbsp;&nbsp;Add COG class and KEGG orthology2.3&nbsp;&nbsp;Uniprot ID mapping3&nbsp;&nbsp;Clean up and save annotation3.1&nbsp;&nbsp;Final statistics3.2&nbsp;&nbsp;Fill missing values<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import sys,os
from Bio import SeqIO
sys.path.append('/home/tahani/Documents/github/modulome/modulome/')
from modulome.config import *
from modulome.util import *
organism = 'Synechococcus_elongatus'
kegg_organism_code = 'syf'
org_dir = os.path.join(DATA_DIR,'organisms',organism)
seq_dir = os.path.join(org_dir,'sequence_files')<jupyter_output><empty_output><jupyter_text>Get information from GFF file  Convert GFF to Pandas DataFrameCode is borrowed from modulome-nextflow to make sure dataframes will work with log-tpm files<jupyter_code>annot_list = []
for filename in os.listdir(seq_dir):
    if filename.endswith('.gff3'):
        gff = os.path.join(seq_dir,filename)
        annot_list.append(gff2pandas(gff))
keep_cols = ['refseq','start','end','strand','gene_name','locus_tag','old_locus_tag','gene_product','ncbi_protein']
DF_annot = pd.concat(annot_list)[keep_cols]
DF_annot.set_index('locus_tag',drop=True,inplace=True)

# Some gene products have commas in the name, which are represented by %2C. 
# Should we change this or keep this as is? (i.e. 2,3-butanediol)
DF_annot.head()<jupyter_output><empty_output><jupyter_text>Check genes<jupyter_code>tpm_file = os.path.join(org_dir,'1_log_tpm_final.csv')
DF_log_tpm = pd.read_csv(tpm_file,index_col=0)
# Check that the genes are the same in the expression dataset as in the annotation dataframe.
# Mismatched genes are listed below
test = DF_annot.sort_index().index == DF_log_tpm.sort_index().index
DF_annot[~test]<jupyter_output><empty_output><jupyter_text>Generate nucleotide fasta files for CDS<jupyter_code>cds_list = []
for filename in os.listdir(seq_dir):
    if filename.endswith('.fasta'):
        fasta = os.path.join(seq_dir,filename)
        seq = SeqIO.read(fasta,'fasta')
        
        # Get gene information for genes in this fasta file
        df_genes = DF_annot[DF_annot.refseq == seq.id]
        for i,row in df_genes.iterrows():
            cds = seq[row.start-1:row.end]
            if row.strand == '-':
                cds = seq[row.start-1:row.end].reverse_complement()
            cds.id = row.name
            cds.description = row.gene_name if pd.notnull(row.gene_name) else row.name
            cds_list.append(cds)
SeqIO.write(cds_list,os.path.join(seq_dir,'CDS.fna'),'fasta')<jupyter_output><empty_output><jupyter_text>Get COGs  Run EggNOG Mapper1. Go to [http://eggnog-mapper.embl.de/].1. Upload the CDS.fna file from your organism directory (within the sequence_files folder)1. Make sure to limit the taxonomy to the correct level1. After the job is submitted, you must follow the link in your email to run the job.1. Once the job completes (after ~4 hrs), download the annotations file.1. Save the annotation file to `data/eggNOG/.annotations`  Add COG class and KEGG orthology<jupyter_code>DF_eggnog = pd.read_csv(os.path.join(DATA_DIR,'eggNOG',organism+'.annotations'),sep='\t',skiprows=4,header=None)
eggnog_cols = ['query_name','seed eggNOG ortholog','seed ortholog evalue','seed ortholog score',
               'Predicted taxonomic group','Predicted protein name','Gene Ontology terms',
               'EC number','KEGG_ko','KEGG_Pathway','KEGG_Module','KEGG_Reaction',
               'KEGG_rclass','BRITE','KEGG_TC','CAZy','BiGG Reaction','tax_scope',
               'eggNOG OGs','bestOG_deprecated','COG','eggNOG free text description']
DF_eggnog.columns = eggnog_cols
# Strip last three rows as they are comments
DF_eggnog = DF_eggnog.iloc[:-3]

# Set locus tag as index
DF_eggnog = DF_eggnog.drop_duplicates('query_name').set_index('query_name')
DF_eggnog.index.name = 'locus_tag'

# Only keep KEGG orthology and COG
DF_eggnog = DF_eggnog[['KEGG_ko','COG']]
DF_eggnog.head()
DF_annot_cog = pd.concat([DF_annot,DF_eggnog],axis=1,sort=False)
DF_annot_cog.index.name = 'locus_tag'
assert(len(DF_annot_cog) == len(DF_annot))<jupyter_output><empty_output><jupyter_text>Uniprot ID mapping<jupyter_code># Try the uniprot ID mapping tool - Use EMBL for Genbank file and P_REFSEQ_AC for Refseq file
mapping_uniprot = uniprot_id_mapping(DF_annot.ncbi_protein.fillna(''),input_id='P_REFSEQ_AC',output_id='ACC',
                             input_name='ncbi_protein',output_name='uniprot')

# Merge with current annotation
DF_annot_cog_uniprot = pd.merge(DF_annot_cog.reset_index(),mapping_uniprot,how='left',on='ncbi_protein')
DF_annot_cog_uniprot.set_index('locus_tag',inplace=True)
assert(len(DF_annot) == len(DF_annot_cog_uniprot))
DF_annot_cog_uniprot.head()<jupyter_output><empty_output><jupyter_text>Clean up and save annotation<jupyter_code># Temporarily remove warning
pd.set_option('mode.chained_assignment', None)
order = ['gene_name','start','end','strand','gene_product','uniprot','COG','KEGG_ko']
DF_annot_final = DF_annot_cog_uniprot[order]
# Strip "ko:" from KEGG orthologies
DF_annot_final['KEGG_ko'] = [item[3:] if isinstance(item,str) else item for item in DF_annot_final['KEGG_ko']]
# Only include one COG group
DF_annot_final['COG'] = [item[0] if isinstance(item,str) else item for item in DF_annot_final['COG']]
DF_annot_final.head()<jupyter_output><empty_output><jupyter_text>Final statistics<jupyter_code>sns.set_style('ticks')
fig,ax = plt.subplots()
DF_annot_final.count().plot(kind='bar',ax=ax)
ax.set_ylabel('# of Values',fontsize=18)
ax.tick_params(labelsize=16)
counts = DF_annot_final.COG.value_counts()
plt.pie(counts.values,labels=counts.index);<jupyter_output><empty_output><jupyter_text>Fill missing values<jupyter_code># Fill in missing gene names with locus tag names
DF_annot_final['tmp_name'] = DF_annot_final.copy().index.tolist()
DF_annot_final.gene_name.fillna(DF_annot_final.tmp_name,inplace=True)
DF_annot_final.drop('tmp_name',axis=1,inplace=True)

# Fill missing COGs with X
DF_annot_final['COG'].fillna('X',inplace=True)
DF_annot_final.to_csv(os.path.join(DATA_DIR,'organisms',organism,'1_gene_info_final.csv'))
#Tahani 8/9/2020 11:31pm<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12612472851380702, 'fraction_numerical': 0.008687558175612783, 'mean_word_length': 6.171301446051168, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '247696', 'license': 'MIT', 'n_tokens_mistral': 2422, 'n_tokens_neox': 2315, 'n_words': 444, 'path': 'Old/Notebooks/compile_gene_annotation.ipynb', 'repo_name': 'talbulus/Synechococcus_elongatus', 'size': '6446'}