<jupyter_start><jupyter_text>GSD: Assessing ambiguous nts in PB 'complete' genomesChecks a collection of PacBio sequenced yeast 'complete' genomes from [Yue et al 2017](https://www.ncbi.nlm.nih.gov/pubmed/28416820) ambiguous/gap-representing residues. Comparing to the SGD reference sequence from [here](http://sgd-archive.yeastgenome.org/sequence/S288C_reference/chromosomes/fasta/) as well. 'Complete' here is meant to refer to the analyses examing the mitochondrial and nuclear combined as one unit, unlike `GSD Assessing_ambiguous_nts_in_nuclear_PB_genomes.ipynb` where was just nuclear genome. And although Yue et al 2017 kept these two genomes separate, this is moving towards where will also analyze the 1011 genome collection where mitochondrial and nuclear genome mixed in one file. As will be seen here, the SGD reference sequence is easily available as one file with mitochondrial genome included; however, `GSD Assessing_ambiguous_nts_in_nuclear_PB_genomes.ipynb` shows it is fairly easy to put together one without the mitochondrial using the FASTA files of the individual chromosomes that are available from a different yeastgenome.org page.References for sequence data:  - [Contrasting evolutionary genome dynamics between domesticated and wild yeasts.<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>. Nat Genet. 2017 Jun;49(6):913-924. doi: 10.1038/ng.3847. Epub 2017 Apr 17. PMID: 28416820](https://www.ncbi.nlm.nih.gov/pubmed/28416820)- [Life with 6000 genes. <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>. Science. 1996 Oct 25;274(5287):546, 563-7. PMID: 8849441](https://www.ncbi.nlm.nih.gov/pubmed/8849441)-----  PreparationGet scripts and sequence data necessary.<jupyter_code>!pip install pyfaidx<jupyter_output>Collecting pyfaidx
  Downloading https://files.pythonhosted.org/packages/75/a5/7e2569527b3849ea28d79b4f70d7cf46a47d36459bc59e0efa4e10e8c8b2/pyfaidx-0.5.5.2.tar.gz
Requirement already satisfied: six in /srv/conda/lib/python3.6/site-packages (from pyfaidx) (1.12.0)
Requirement already satisfied: setuptools>=0.7 in /srv/conda/lib/python3.6/site-packages (from pyfaidx) (40.6.3)
Building wheels for collected packages: pyfaidx
  Running setup.py bdist_wheel for pyfaidx ... [?25ldone
[?25h  Stored in directory: /home/jovyan/.cache/pip/wheels/54/a2/b4/e242e58d23b2808e191b214067880faa46cd2341f363886e0b
Successfully built pyfaidx
Installing collected packages: pyfaidx
Successfully installed pyfaidx-0.5.5.2<jupyter_text>Get the genomes data by running these commands.<jupyter_code>import pandas as pd
# Prepare for getting PacBio (Yue et al 2017 sequences)
#make a list of the strain designations
yue_et_al_strains = ["S288C","DBVPG6044","DBVPG6765","SK1","Y12",
                     "YPS128","UWOPS034614","CBS432","N44","YPS138",
                     "UFRJ50816","UWOPS919171"]
# Get & unpack the genome sequences from strains 
for s in yue_et_al_strains:
    !curl -LO http://yjx1217.github.io/Yeast_PacBio_2016/data/Nuclear_Genome/{s}.genome.fa.gz
    !curl -OL http://yjx1217.github.io/Yeast_PacBio_2016/data/Mitochondrial_Genome/{s}.mt.genome.fa.gz
    !gunzip -f {s}.genome.fa.gz
    !gunzip -f {s}.mt.genome.fa.gz
    # add the mitochondrial genome content onto the nuclear
    !cat {s}.genome.fa {s}.mt.genome.fa > temp.genome.fa
    !mv temp.genome.fa {s}.genome.fa
!rm *.mt.genome.fa
# add identifiers to each `chr` so results for each strain clear later
chromosome_id_prefix = "chr"
def add_strain_id_to_description_line(file,strain_id):
    '''
    Takes a file and edits every description line to add 
    strain_id after the caret.
    
    Saves the fixed file
    '''
    import sys
    output_file_name = "temp.txt"
    # prepare output file for saving so it will be open and ready
    with open(output_file_name, 'w') as output_file:

        # read in the input file
        with open(file, 'r') as input_handler:
            # prepare to give feeback later or allow skipping to certain start
            lines_processed = 0

            for line in input_handler:
                lines_processed += 1
                if line.startswith(">"):
                    rest_o_line = line.split(">")
                    new_line = ">"+strain_id + rest_o_line[1]
                else:
                    new_line = line
                
                # Send text to output
                output_file.write(new_line)

    
    # replace the original file with edited
    !mv temp.txt {file}
    # Feedback
    sys.stderr.write("\n{} chromosome identifiers tagged.".format(file))

for s in yue_et_al_strains:
    add_strain_id_to_description_line(s+".genome.fa",s)
# Get SGD reference sequence that includes nuclear and mitochondrial sequence as one file,
# among others. I'll use file name for the reference genome worked out 
# in `GSD Assessing_ambiguous_nts_in_nuclear_PB_genomes.ipynb`, so more of the
# previously worked out code will work.
!curl -OL http://sgd-archive.yeastgenome.org/sequence/S288C_reference/genome_releases/S288C_reference_genome_Current_Release.tgz
!tar -xzf S288C_reference_genome_Current_Release.tgz
!rm S288C_reference_genome_Current_Release.tgz
!mv S288C_reference_genome_R64-2-1_20150113/S288C_reference_sequence_R64-2-1_20150113.fsa ./SGD_REF.genome.fa
!rm -rf S288C_reference_genome_R64-2-1_20150113
# Make a list of all `genome.fa` files, excluding `genome.fa.nhr` and `genome.fa.nin` and `genome.fansq`
# The excluding was only necessary because I had run some BLAST queries preiminarily in development. Normally, 
# it would just be the `.re.fa` at the outset. (But keeping because removal ability could be useful.)
fn_to_check = "genome.fa" 
genomes = []
import os
import fnmatch
for file in os.listdir('.'):
    if fnmatch.fnmatch(file, '*'+fn_to_check):
        if not file.endswith(".nhr") and not file.endswith(".nin") and not file.endswith(".nsq") :
            genomes.append(file)
genomes<jupyter_output><empty_output><jupyter_text>Now you are prepared to analyze each PacBio-sequenced genome.  Assessing the genomes in regards to stretches of N base calls Total number of Ns present<jupyter_code>from pyfaidx import Fasta
query = "N"
#query = "X"
#query = "A" # CONTROL. should occur!
for g in genomes:
    chrs = Fasta(g)
    ambiguous = []
    for x in chrs:
        #print(x.name)
        if query in str(x):
            ambiguous.append(x.name)
    #print(ambiguous)
    if not ambiguous:
        print ("No ambiguous nucleotides or gaps in {}.".format(g))
    else:
        print("There are occurences of {} in {}.".format(query,g))<jupyter_output>There are occurences of N in DBVPG6044.genome.fa.
No ambiguous nucleotides or gaps in SGD_REF.genome.fa.
There are occurences of N in DBVPG6765.genome.fa.
There are occurences of N in CBS432.genome.fa.
There are occurences of N in YPS128.genome.fa.
There are occurences of N in UFRJ50816.genome.fa.
There are occurences of N in S288C.genome.fa.
There are occurences of N in SK1.genome.fa.
There are occurences of N in UWOPS034614.genome.fa.
There are occurences of N in UWOPS919171.genome.fa.
There are occurences of N in YPS138.genome.fa.
There are occurences of N in N44.genome.fa.
There are occurences of N in Y12.genome.fa.<jupyter_text>Another way to assess, just count all the letters present:<jupyter_code>%%time
from pyfaidx import Fasta
import pandas as pd
import collections
nt_counts = {}
for g in genomes:
    strain_id = g.split(".genome.fa")[0]
    concatenated_seqs = ""
    chrs = Fasta(g)
    for x in chrs:
        #print(x.name)
        concatenated_seqs += str(x)
    nt_counts[strain_id] = collections.Counter(concatenated_seqs)
nt_count_df = pd.DataFrame.from_dict(nt_counts, orient='index').fillna(0)
nt_count_df["Total_nts"] = nt_count_df.sum(1)
def percent_calc(items):
    '''
    takes a list of two items and calculates percentage of first item
    within total (second item)
    '''
    return items[0]/items[1]
nt_count_df['% N'] = nt_count_df[['N','Total_nts']].apply(percent_calc, axis=1)
nt_count_df = nt_count_df.style.format({'Total_nts':'{:.2E}','% N':'{:.2%}'})
nt_count_df<jupyter_output><empty_output><jupyter_text>Examining the number and size of stretches of Ns present Make a dataframe of the data and also make a summary one with the say, top five(?) number instances, for each.<jupyter_code>%%time
# count frequency of blocks of Ns in the genomes
import re
from collections import defaultdict
from pyfaidx import Fasta
import pandas as pd
import collections


min_number_Ns_in_row_to_collect = 3
pattern_obj = re.compile("N{{{},}}".format(min_number_Ns_in_row_to_collect), re.I)  # adpated from
# code worked out in `collapse_large_unknown_blocks_in_DNA_sequence.py`, which relied heavily on
# https://stackoverflow.com/a/250306/8508004

len_match_dict_by_strain = {}

#genomes = ["N44.genome.fa"]

for g in genomes:
    len_match_dict = defaultdict(int)
    strain_id = g.split(".genome.fa")[0]
    records = Fasta(g)
    for record in records:
        for m in pattern_obj.finditer(str(record)):
            len_match_dict[len(m.group())] += 1
    #before passing final dictionary for strain to
    #collection, add an entry of zero for size of 1 so the
    #strain with no stretches will be in final dataframe
    if not len_match_dict:
        len_match_dict[1]=0
    len_match_dict_by_strain[strain_id] = len_match_dict

#stretches_size_freq_df = pd.DataFrame.from_dict(len_match_dict_by_strain, orient='index').fillna(0)
#stretches_size_freq_df = pd.DataFrame.from_dict(len_match_dict_by_strain).fillna(0).stack().reset_index() # if wanted all stretches of size to show count even if none present
stretches_size_freq_df = pd.DataFrame.from_dict(len_match_dict_by_strain).stack().reset_index() 
stretches_size_freq_df.columns = ['stretch_size','strain','#_instances'] 
stretches_size_freq_df = stretches_size_freq_df[['strain','stretch_size','#_instances']]
stretches_size_freq_df = stretches_size_freq_df.sort_values(['stretch_size','#_instances'],ascending=[0,0]).reset_index(drop=True)<jupyter_output>CPU times: user 3.16 s, sys: 0 ns, total: 3.16 s
Wall time: 3.43 s<jupyter_text>Because of use of `%%time` in above cell, easiest to display resulting dataframe using next cell.<jupyter_code>stretches_size_freq_df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10394990366088631, 'fraction_numerical': 0.04046242774566474, 'mean_word_length': 4.232358870967742, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 56, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 2, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '100212', 'license': 'MIT', 'n_tokens_mistral': 3911, 'n_tokens_neox': 3439, 'n_words': 988, 'path': 'notebooks/GSD/GSD Assessing_ambiguous_nts_in_complete_PB_genomes.ipynb', 'repo_name': 'fomightez/cl_sq_demo-binder', 'size': '10465'}