<jupyter_start><jupyter_text>This notebook was created to find all of the High Quality SNV calls for all of the isolates in a DataBase with ~33k VCF files from processed FASTQ files that were downloaded from NCBI<jupyter_code>import vcf

%matplotlib inline
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.ticker as ticker
from itertools import compress
from pylab import MaxNLocator
import seaborn as sns; sns.set()
from matplotlib.colors import LogNorm
from matplotlib import gridspec
import ast
import itertools
import seaborn as sns
from sklearn.preprocessing import StandardScaler

import fastcluster
from sklearn import cluster, datasets
import scipy.cluster.hierarchy as hier
from sklearn.cluster import KMeans
import time
import sys

import Bio
from Bio.Alphabet import IUPAC
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio.Blast import NCBIXML
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio import pairwise2
from Bio import SeqIO
from Bio.Graphics import GenomeDiagram
from Bio.SeqUtils import GC

from Bio.Align.Applications import MuscleCommandline
from StringIO import StringIO
from Bio import AlignIO
from Bio.Align import AlignInfo
from Bio.Seq import MutableSeq
import itertools
import gzip

import networkx as nx
import scipy
import pickle<jupyter_output><empty_output><jupyter_text>[1] This section of the notebook loads all of the Lineage Calls that were made on the isolates in rollingDB as well as the mean coverage / isolates for all of the isolates in rollingDB.    We runs through the pilon-generated VCF for each isolate (N = 33,875) and scrapes the High Quality Single Nucleotide Variant calls made against H37Rv. The DataFrame containing these high quality SNV calls (with the Reference Position, Reference Allele and Alterante Allele) for each isolate is then stored in a directory containing all of the DFs for all isolates.   Side-Note: Some isolate have ~20x coverage but we require 25 reads at a given position to make a SNP call  Scrape Luca's DB<jupyter_code>#genomic data directory
rolling_DB_dir = '/n/data1/hms/dbmi/farhat/rollingDB/genomic_data/'

#get all folders (each folder corresponds to a different sequenced isolate)
isolate_directories = os.listdir(rolling_DB_dir)

#iterate through each sequenced isolate
isolate_i = 0

for isolate_ID in isolate_directories:

    #directory that stores files for each sequenced isolate
    directory_for_sequenced_isolate = rolling_DB_dir + isolate_ID

    #check to see if megapipe successfully ran on sequenced isolate
    try:
        #existence of a PILON directory and corresponding VCF file
        if ( 'pilon' in os.listdir(directory_for_sequenced_isolate) ):
            
            #existence of a REDUCED VCF and existence of FULL VCF file in relevent directories 
            if ( 'vcf' in list( itertools.chain.from_iterable( [filename.split('.') for filename in os.listdir(directory_for_sequenced_isolate + '/pilon/')] ) ) ) and ( isolate_ID + '_full' in list( itertools.chain.from_iterable( [filename.split('.') for filename in os.listdir(directory_for_sequenced_isolate + '/pilon/')] ) ) ):
                
                #we have a valid VCF file so megapipe ran successfully, let's keep the variant call information for this sequenced isolate
                
                #PILON-generated VCF DATA
                ########################################################################################################################
                
                #dictionaries to construct DataFrame from later that hold information about each high-quality SNP found for a specific isolate from rollingDB
                ref_dict = {} #reference allele
                alt_dict = {} #alternate allele
                alt_AF_dict = {} #frequency of alternate allele in population (0.90 - 1.0)
                
                #load in VCF file if one was found
                vcf_reader = vcf.Reader( open(directory_for_sequenced_isolate + '/pilon/' + isolate_ID + '.vcf' , 'r') )

                #iterate through each Variant Call 
                for record in vcf_reader:
                    
                    #check to see if the call is a PASS by Pilon
                    if record.FILTER == []:

                        #check to see if variant is SNP (length of Reference Allele is 1 base ; there's only 1 alternate allele ; length of alternate allele is 1 base)
                        if (len(record.REF) == 1) and (len(record.ALT) == 1) and (len(str(record.ALT[0])) == 1):

                            ##### Retrieve Relevant information for filtering quality of Base Call #####
                            # Mean Base Quality @ locus
                            BQ = record.INFO['BQ']
                            # Mean Mapping Quality @ locus
                            MQ = record.INFO['MQ']
                            # Number of Reads w/ Deletion 
                            DC = record.INFO['DC']
                            # Number of Reads w/ Insertion
                            IC = record.INFO['IC']
                            # Depth of Valid Reads in Pileup
                            VD = record.INFO['DP']

                            ### Filtering Criteria
                            #---> Mean Base Quality > 20
                            #---> Mean Mapping Quality > 30
                            #---> No Reads Supporting Insertions
                            #---> No Reads Supporting Deletions
                            #---> Number of High Quality Reads >= 25
                            if (BQ > 20) and (MQ > 30) and (DC == 0) and (IC == 0) and (VD >= 25): #SNP passed filtering criteria!

                                # Filtering Criteria for mutant allele frequency calculation (range: 0.90 - 1.0)
                                ref_allele = str( record.REF )
                                alt_allele = str( record.ALT[0] )
                                ref_pos = int( record.POS )
                                
                                # Calculate Alternate Allele Frequency
                                ########################################################################
                                #position-base dictionaries (order of output in Pilon)
                                base_order = {0:'A' , 1:'C' , 2:'G', 3:'T'}
                                base_order_r = {'A':0 , 'C':1 , 'G':2 , 'T':3}

                                # Quality-Percentage - Percentage of As, Cs, Gs, Ts weighted by Q & MQ at locus
                                QP = record.INFO['QP']

                                if (sum(np.array(QP) == 0) == 3): #1 Allele supported by Reads (Alternate), store information for this Base Call (3 bases have NO reads supporting)

                                    #create a Boolean list for the Bases that have at least 1% of support from reads
                                    supported_bases_boolean = list( np.array( QP ) > 0 )

                                    #positions of Bases [A,C,G,T] with at least 1% of support from reads
                                    supported_bases_i = [i for i,x in enumerate(supported_bases_boolean) if x == True] #returns list of length 
                                    supported_base = base_order[supported_bases_i[0]]

                                    #if the Reference Allele is supported
                                    if supported_base == ref_allele:
                                        alt_allele_frequency = 0.0 #no alternate allele

                                    #if an Alternate Allele is supported:
                                    if supported_base != ref_allele:
                                        alt_allele_frequency = 1.0 #alternate allele

                                elif (sum(np.array(QP) == 0) == 2): #2 Alleles supported, make sure 1 is the Reference & confirm the Alternate (2 bases have NO reads supporting)

                                    #position of the Reference Allele in [A,C,G,T]
                                    ref_allele_base_i = base_order_r[ref_allele] 

                                    #create a Boolean list for the Bases that have at least 1% of support from reads
                                    supported_bases_boolean = list( np.array( QP ) > 0 )

                                    #positions of Bases [A,C,G,T] with at least 1% of support from reads
                                    supported_bases_i = [i for i,x in enumerate(supported_bases_boolean) if x == True] #returns list of length 2

                                    #check that Reference Allele is supported by at least 1% of the reads
                                    if ref_allele_base_i in supported_bases_i:

                                        #find the Alternate Allele that was supported by at least 1% of the reads
                                        for supported_base_i in supported_bases_i:

                                            if supported_base_i != ref_allele_base_i: #then must be alternate allele

                                                alt_allele = base_order[supported_base_i] #Alternate Allele
                                                alt_allele_frequency = float(QP[supported_base_i]) / 100.0 #weighted percentage of reads supporting Alternate Allele
                                                break
                                ########################################################################
                                
                                # use only consensus SNPs (where alternate allele frequency >= 90%)
                                if alt_allele_frequency >= 0.90:
                                
                                    # Store good-quality SNV in SNV dictionaries
                                    ref_dict[ref_pos] = ref_allele
                                    alt_dict[ref_pos] = alt_allele

                #construct DataFrame from dictionaries
                SNV_DF = pd.DataFrame( index = pd.Series( ref_dict.keys() ) )
                SNV_DF['ref'] = pd.Series(ref_dict)
                SNV_DF['alt'] = pd.Series(alt_dict)

                #make reference-position a column and re-index, sort according to reference position
                SNV_DF['pos'] = SNV_DF.index
                SNV_DF.sort_values(by = 'pos' , ascending = True , inplace = True)
                SNV_DF.reset_index(inplace = True , drop = True)

                #pickle DataFrame holding High Quality SNVs from isolate
                SNV_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB/' + isolate_ID + '_SNVs.pkl')
                ########################################################################################################################
                
                #keep track of progress   
                isolate_i += 1
                if isolate_i % np.ceil(0.05*len(isolate_directories)) == 0:
                    print float(isolate_i) / float(len(isolate_directories))
        
    except OSError: #hit some file that is not another directory with genomic data 
        continue<jupyter_output>0.0500090232801
0.10001804656
0.15002706984
0.20003609312
0.2500451164
0.30005413968
0.35006316296
0.400072186241
0.450081209521
0.500090232801
0.550099256081
0.600108279361
0.650117302641<jupyter_text>[2] This section of the notebook constructs a script that iterates over *subsets* of the HQ SNV DataFrames for each isolate and creates a 'catalogue' of all of the HQ SNVs found in rollingDB. All of the relevent information is stored in dictionaries before being transformed into a DataFrame with all of the unique SNVs (indexed by Reference Position) and the number of isolates that contain a SNV at that position.<jupyter_code>import vcf
import os
import pandas as pd
import numpy as np
from pylab import MaxNLocator
import ast
import itertools
import time
import sys
import itertools
import gzip
import scipy<jupyter_output><empty_output><jupyter_text>We are going to split this up into several jobs, each iterating through N isolates. Each job will produce a DF with all SNVs detected in any of the N isolates & corresponding alleles.  Once all jobs have finished, we will merge all DFs together to create one master DF of all SNVs detected in any of the 11,034 isolates in rollingDB<jupyter_code>#INPUTS
first_isolate_index = int(sys.argv[1]) #INPUT - first isolate index
last_isolate_index = int(sys.argv[2]) #INPUT - last isolate index
job_number = str(sys.argv[3]) #INPUT - job number
#list of pickles DataFrame files for SNVs for each isolate in rollingDB
SNV_DF_list = os.listdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB')

#subset to only the pickled DFs for the isolates that will be searched for SNVs in this job
SNV_DF_list = SNV_DF_list[first_isolate_index:last_isolate_index]
def store_SNV_info_from_SNV(SNV_DF_row):

    ref_allele , alt_allele , ref_pos = list(SNV_DF_row)

    #SNV is not documented at this Reference Position
    if ref_pos not in ref_dict.keys():

        #store a new row for SNV at this reference position; ref allele, alt allele 1, nan, nan, 1 isolate has polymorphism here (so far)
        ref_dict[ref_pos] = ref_allele
        alt1_dict[ref_pos] = alt_allele
        alt2_dict[ref_pos] = np.nan
        alt3_dict[ref_pos] = np.nan
        num_isolates_dict[ref_pos] = 1

    #a SNV is already documented at this Reference Position
    else:

        #check to see if specific allele is documented, if not add it
        if alt_allele not in [alt1_dict[ref_pos] , alt2_dict[ref_pos] , alt3_dict[ref_pos]]:

            #need to document presence of new allele in isolate pool (either going to be alt2 or alt3)
            if ( type(alt2_dict[ref_pos])==float ) and ( np.isnan(alt2_dict[ref_pos]) ): #alt2 has not been found yet

                alt2_dict[ref_pos] = alt_allele

            elif ( type(alt3_dict[ref_pos])==float ) and ( np.isnan(alt3_dict[ref_pos]) ): #alt3 has not been found yet

                alt3_dict[ref_pos] = alt_allele

        #increment # of isolates that have a polymorphism here
        num_isolates_dict[ref_pos] = num_isolates_dict[ref_pos] + 1
#create dictionaries to store SNV info for each isolate
ref_dict = {}
alt1_dict = {}
alt2_dict = {}
alt3_dict = {}
num_isolates_dict = {}

#iterate through all DataFrames for each isolate in rollingDB
for pickled_SNV_DF in SNV_DF_list:
    
    #load SNV DF for isolates from pickled object
    SNV_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB/' + pickled_SNV_DF)

    #extract important info from SNV df
    _ = SNV_DF.apply(lambda row: store_SNV_info_from_SNV(row) , axis = 1 )  
#construct DataFrame from dictionaries
SNV_DF = pd.DataFrame( index = pd.Series( ref_dict.keys() ) )
SNV_DF['ref'] = pd.Series(ref_dict)
SNV_DF['alt1'] = pd.Series(alt1_dict)
SNV_DF['alt2'] = pd.Series(alt2_dict)
SNV_DF['alt3'] = pd.Series(alt3_dict)
SNV_DF['num_isolates'] = pd.Series(num_isolates_dict)

#make reference-position a column and re-index
SNV_DF['pos'] = SNV_DF.index

#sort according to reference position
SNV_DF.sort_values(by = 'pos' , ascending = True , inplace = True)

SNV_DF.reset_index(inplace = True , drop = True)<jupyter_output><empty_output><jupyter_text>Output Pickle files for downstream analysis<jupyter_code>SNV_DF.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/seperate_HQ_SNV_DFs/HighQuality_SNVs_found_in_rollingDB_catalogue_' + job_number + '.pkl' )<jupyter_output><empty_output><jupyter_text>[3] This section of the notebook launches the script that was created in the previous section of the notebook.<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np
import itertools
import shutil
#list of all SNV DFs for each isolate that was run through MegaPipe in rollingDB
SNV_DF_list = os.listdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB')

isolates_per_job = 500 #number of isolates each job will catalogue
njobs = int( np.ceil( float( len( SNV_DF_list ) ) / float(isolates_per_job) ) ) #number of jobs required if we split for every 500 isolates
print njobs #number of jobs we will submit
first_isolate_index = 0
max_isolate_index = int( len( SNV_DF_list ) )

#run each job using a difference set of indices for isolates that we will catalogue SNVs from
for job_number in np.arange(1 , njobs+1):

    last_isolate_index = first_isolate_index + isolates_per_job
    if last_isolate_index > max_isolate_index: #watch out for the last job
        last_isolate_index = max_isolate_index
    
    rollingDB_SNV_cataloging_job = 'python /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/python_scripts/rollingDB_HQ_SNV_catalogue_scrape.py ' + str(first_isolate_index) + ' ' + str(last_isolate_index) + ' ' + str(job_number)
    
    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/SNV_catalogue_scrape_jobs/')

    job_name = 'SNV_' + str(job_number)

    s = Slurm(job_name , {'partition':'short' , 'n':'1' , 't':'0-01:00:00' , 'mem-per-cpu':'4G' , 'mail-type':'FAIL' , 'mail-user':'<EMAIL>'})

    #submits the job
    job_id = s.run(rollingDB_SNV_cataloging_job)

    print job_name  + ' : ' +  str(job_id)
    
    #setup isolate indices for next job
    first_isolate_index = last_isolate_index<jupyter_output>submitted: Submitted batch job 3299559
submitted: Submitted batch job 3299560<jupyter_text>[4] The purpose of this section of the notebook is generate the master catalogue of (high-quality) SNV calls made in rollingDB by merging the catalogues of SNVs found in seperate batches of isolates.<jupyter_code>import os
import pandas as pd
import numpy as np
import ast
import itertools
import time
import sys
import itertools
#INPUTS
njobs = 68 #number of seperate SNV catalogues for subsets of isolates in rollingDB<jupyter_output><empty_output><jupyter_text>Load in SNV catalogue(s) (High Quality SNV calls across all isolates in RollingDB) We're going to iterate through each DataFrame (catalogue of SNVs detected within each batch of isolates) and create a master DF of all SNVs detected within all the isolates in rollingDB.<jupyter_code>def convert_SNV_DF_rows_to_DICT_elments(row):
    
    '''
    This function converts the row of a SNV catalogue DF to a
    list with two elements: (1) The Reference Position of an SNV detected within a subset of isolates
    and (2) another list containing the Reference Allele, all Alternate Alleles, and Number of Isolates
    with an SNV at this Reference Position
    '''
    
    #create set of alternate alleles present for this SNV
    alt_alleles_SNV_i = set( [row['alt1'] , row['alt2'] , row['alt3']] )
    #purge nan values
    alt_alleles_SNV_i = {x for x in alt_alleles_SNV_i if x==x}
        
    return [ row['pos'] , [row['ref'] , alt_alleles_SNV_i , row['num_isolates']] ]

def convert_SNV_DF_to_DICT(SNV_DF):
    
    '''
    This function takes an SNV DF catalogue for a subset of isolates and converts it into 
    an SNV DICT catalogue using the function above by applying it to each row.
    '''
    
    SNV_dict = SNV_DF.apply( convert_SNV_DF_rows_to_DICT_elments , axis = 1 )
    return dict(list( SNV_dict.values ))<jupyter_output><empty_output><jupyter_text>Convert SNV catalogue(s) DFs to SNV catalogue DICTs and store SNP catalogue info across jobs<jupyter_code>#master set of all reference positions detected within rollingDB (faster to query a set)
SNV_set_Ref_Position = set([])

#master dictionaries of all SNVs detected within rollingDB
SNV_dict_Ref_Allele = {} #keys - reference position; values - ref_allele (string)
SNV_dict_Alt_Allele = {} #keys - reference position; values - alt_alleles (set)
SNV_dict_Num_Isolates = {} #keys - reference position; values - num_isolates (int)

for job_num in np.arange(1 , njobs+1):
    
    #load the SNV DataFrame holding SNVs from the batch of isolates
    SNV_DF_subset = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/seperate_HQ_SNV_DFs/HighQuality_SNVs_found_in_rollingDB_catalogue_' + str(job_num) + '.pkl' )
    
    #convert SNV DataFrame to an SNV Dict
    SNV_dict_subset = convert_SNV_DF_to_DICT(SNV_DF_subset)
    
    #iterate through each SNV in Dict
    for ref_pos in SNV_dict_subset.keys():
        
        #update master SNV catalogue dict with SNV info
        ref , alt_alleles_SNV_i , num_isolates = SNV_dict_subset[ref_pos]
    
        #[A] reference position not stored in master SNV dictionaries
        ##if ref_pos not in SNV_dict_Ref_Allele.keys():
        if ref_pos not in SNV_set_Ref_Position:
            
            SNV_set_Ref_Position.add(ref_pos)
            
            SNV_dict_Ref_Allele[ref_pos] = ref
            SNV_dict_Alt_Allele[ref_pos] = alt_alleles_SNV_i
            SNV_dict_Num_Isolates[ref_pos] = num_isolates
        
        #[B] reference position already stored in master SNV dictionaries
        else:
            
            #update SNV info with new SNV info (ref_allele stays the same)
            SNV_dict_Alt_Allele[ref_pos] = SNV_dict_Alt_Allele[ref_pos].union(alt_alleles_SNV_i) #store any new alleles detected
            SNV_dict_Num_Isolates[ref_pos] = SNV_dict_Num_Isolates[ref_pos] + num_isolates #add number of isolates
            
    #keep track of progress
    print job_num<jupyter_output>1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68<jupyter_text>Convert master SNV dict to SNV df<jupyter_code>#convert dictionaries into pandas series objects
SNV_series_Ref_Allele = pd.Series(SNV_dict_Ref_Allele , name = 'ref')
SNV_series_Alt_Allele = pd.Series(SNV_dict_Alt_Allele , name = 'alt')
SNV_series_Num_Isolates = pd.Series(SNV_dict_Num_Isolates , name = 'num_isolates')

#construct DataFrame from Series objects
SNV_DF = pd.DataFrame([SNV_series_Ref_Allele , SNV_series_Alt_Allele , SNV_series_Num_Isolates]).T

#convert Alternate Allele(s) column to list
SNV_DF['alt'] = SNV_DF.alt.apply(list)

#create a column for reference position
SNV_DF['pos'] = list(SNV_DF.index)

#re-index DataFrame
SNV_DF.reset_index(inplace = True , drop = True)

#store master SNV DataFrame as a pickled object
SNV_DF.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/all_HighQuality_SNVs_found_in_rollingDB_catalogue.pkl' )<jupyter_output><empty_output><jupyter_text>Use to load pickle<jupyter_code>SNV_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/all_HighQuality_SNVs_found_in_rollingDB_catalogue.pkl')
SNV_DF.head()
SNV_DF.tail()
np.shape(SNV_DF)<jupyter_output><empty_output><jupyter_text>[5] The purpose of this section of the notebook is to analyze the catalogue of (high-quality) SNV calls made in rollingDB. For population genetic analyses, we will retain only SNVs that are present in at least 2 isolates. We then construct a skeleton genotype DF (rows = Reference Position , columns = isolates) for all of the remaining SNVs.   Construct a series for the isolate annotation<jupyter_code>HQ_SNV_per_isolate_files = os.listdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB/')
isolate_ids_list = [isolate_id.split('_SNV')[0] for isolate_id in HQ_SNV_per_isolate_files]
isolate_id_series = pd.Series(isolate_ids_list)
isolate_id_series.head() #index = column indices of genotype matrix, values isolate IDs
np.shape(isolate_id_series)<jupyter_output><empty_output><jupyter_text>Remove 2 bad isolates from isolate annotation **Remove:** Delete ID for isolate **IDR1400002533** since Full GZipped VCF looks corrupt and remove from list accordingly!<jupyter_code>isolate_id_series = isolate_id_series[isolate_id_series != 'IDR1400002533']<jupyter_output><empty_output><jupyter_text>**Remove:** Delete ID for **Peru4994** since F2 calculations keeps failing for this isolate (lineage-defining SNP depths df is empty) **Maybe don't need to remove this isolate?**<jupyter_code>isolate_id_series = isolate_id_series[isolate_id_series != 'Peru4994']
isolate_id_series.reset_index(inplace = True, drop = True)
isolate_id_series.head()
len(isolate_id_series)
isolate_id_series.to_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl' )<jupyter_output><empty_output><jupyter_text>Load in catologue of SNVs detected within rollingDB<jupyter_code>#SNV data for all isoaltes in rollingDB 
SNV_DF = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/all_HighQuality_SNVs_found_in_rollingDB_catalogues/all_HighQuality_SNVs_found_in_rollingDB_catalogue.pkl' )
SNV_DF.head()
np.shape(SNV_DF)<jupyter_output><empty_output><jupyter_text>How many SNPs are not biallelic?<jupyter_code>#SNPs that are not biallelic will have more than 1 alternate allele
np.sum( SNV_DF.alt.apply(len) > 1 )
SNV_DF[SNV_DF.alt.apply(len) > 1].head(n = 10)
plt.style.use('ggplot')
plt.rcParams['lines.linewidth']=3
plt.rcParams['axes.facecolor']='1.0'
plt.rcParams['xtick.color']='black'
plt.rcParams['axes.grid']=True
plt.rcParams['axes.edgecolor']='black'
plt.rcParams['grid.color']= '1.0'
plt.rcParams.update({'font.size': 12})
fig , ax = plt.subplots()

n , bins , patches = plt.hist(SNV_DF.num_isolates.astype(float) , bins = 60 , color = 'k' , rwidth = 1.0 , edgecolor='white', linewidth=1.35)

plt.title('Number of Isolates Each SNP is Found in', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.xlabel('Number of Isolates', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.ylabel('Number of SNPs', fontweight = 'bold' , fontsize = 12, color = 'k')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.tick_params(labelcolor = 'k')
ax.set_yscale("log", nonposy='clip')

fig = plt.gcf()
fig.set_size_inches(7.5, 5.25)
fig.tight_layout()

plt.tick_params(axis='y', which='major', labelsize=12 , labelcolor = 'k')
plt.tick_params(axis='x', which='major', labelsize=12 , labelcolor = 'k')

#file_name = '/home/rv76/Farhat_Lab/Figures/Final_5/Scraping_rollingDB/Mean_Coverage_Distribution_across_all_isolates.pdf'
#plt.savefig(file_name, bbox_inches='tight', dpi = 300 , transparent = True)

plt.show()<jupyter_output><empty_output><jupyter_text>How many SNP sites have at least *1* isolate with an alternate allele?<jupyter_code>np.shape(SNV_DF[SNV_DF.num_isolates >= 1])<jupyter_output><empty_output><jupyter_text>How many SNP sites have at least *2* isolates with an alternate allele?<jupyter_code>np.shape(SNV_DF[SNV_DF.num_isolates >= 2])<jupyter_output><empty_output><jupyter_text>How many SNP sites have at least *3* isolates with an alternate allele?<jupyter_code>np.shape(SNV_DF[SNV_DF.num_isolates >= 3])<jupyter_output><empty_output><jupyter_text>How many SNP sites have at least *4* isolates with an alternate allele?<jupyter_code>np.shape(SNV_DF[SNV_DF.num_isolates >= 4])<jupyter_output><empty_output><jupyter_text>How many SNP sites have at least *5* isolates with an alternate allele?<jupyter_code>np.shape(SNV_DF[SNV_DF.num_isolates >= 5])<jupyter_output><empty_output><jupyter_text>Keep only SNVs in which the corresponding *Reference Position* has a SNV in $\ge$ 5 isolates SNVs founds in only 1-4 isolate(s) are dropped<jupyter_code>SNV_DF = SNV_DF[SNV_DF.num_isolates >= 5]
np.shape(SNV_DF)
SNV_DF.head()<jupyter_output><empty_output><jupyter_text>Re-Structure DataFrame to include only the Reference Position, Reference Allele and Alternate Allele(s)<jupyter_code>SNV_annotation_DF = pd.DataFrame(index = range( 0 , np.shape(SNV_DF)[0] ) ) #index for each SNV
SNV_annotation_DF['pos'] = list(SNV_DF.pos) #H37Rv Reference Position
SNV_annotation_DF['ref'] = list(SNV_DF.ref) #reference allele
SNV_annotation_DF['alt'] = list(SNV_DF.alt) #alterante allele(s) list
np.shape(SNV_annotation_DF)
SNV_annotation_DF.head() #index = row indices of genotype matrix , 0 supports ref, 1 supports alt<jupyter_output><empty_output><jupyter_text>This is the DataFrame which we'll use to extract SNVs of interest when iterating over all of the VCF files in rollingDB Output Annotation for SNVs of interest in Genotypes folder<jupyter_code>SNV_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_SNP_annotation.pkl')<jupyter_output><empty_output><jupyter_text>Create a text file with a single column of all Reference Positions that we're interested in based off of DF above. We will use this text file to subset the *full* VCF files to only the positions that Reference Positions that have a SNV in at least 2 isolates.<jupyter_code>SNV_ref_positions = ['NC_000962.3\t' + ref_pos for ref_pos in list( SNV_annotation_DF.pos.astype(str) ) ]

#save each element as a new row in a text file
with open('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/ref_positions_to_select.txt', 'w') as f:
    for item in SNV_ref_positions:
        f.write("%s\n" % item)<jupyter_output><empty_output><jupyter_text>[6] This section of the notebook submits jobs to subset Full VCFs to the 146,874 reference positions that have a SNV in $\ge 5$ isolates in rollingDB.<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np
import itertools
import shutil<jupyter_output><empty_output><jupyter_text>Every isolate in the pool of **33,873** isolates in rollingDB had a reduced VCF file and a corresponding Full VCF file. We will subset each of these Full VCF files and create subset VCF files containing information for only the **146,874** reference positions that we're interested in.<jupyter_code>isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl' )
isolate_id_series.head()
len(isolate_id_series)<jupyter_output><empty_output><jupyter_text>Create directory for each isolate to store output for genotyping<jupyter_code>for tag in isolate_id_series.values:

    #create directory to store subset VCF and genotypes for each isolate
    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/' + tag
    if os.path.exists(isolate_genotype_dir):
        shutil.rmtree(isolate_genotype_dir)
        os.makedirs(isolate_genotype_dir)
    elif not os.path.exists(isolate_genotype_dir):
        os.makedirs(isolate_genotype_dir)<jupyter_output><empty_output><jupyter_text>Subset VCF to reference positions of interest for isolates that had Full VCF files  Create BASH scripts to extract VCFs of N isolates per job<jupyter_code>N_isolates_per_job = 70

njobs = int( np.ceil( float( len( isolate_id_series ) ) / float(N_isolates_per_job) ) ) #number of jobs required if we split for every N isolates
print njobs
bash_scripts = []

#genomic data directory
rolling_DB_dir = '/n/data1/hms/dbmi/farhat/rollingDB/genomic_data/'

#split jobs up into tasks of N
for first_isolate_i_for_job in np.arange(0 , (njobs)*N_isolates_per_job , N_isolates_per_job):

    last_isolate_i_for_job = first_isolate_i_for_job + N_isolates_per_job
    
    commands_list = []
    commands_list.append( 'module load bcftools/1.9' ) #only need to load this once per job

    for isolate_i in range(first_isolate_i_for_job , last_isolate_i_for_job): 

        try:
            tag = isolate_id_series[isolate_i]
            zipped_VCF = rolling_DB_dir + tag + '/pilon/' + tag + '_full.vcf.gz'
            isolate_genotypes_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/' + tag

            commands_list.append( 'cd {}'.format(isolate_genotypes_dir) )
            commands_list.append( 'time bcftools view {} -O b -o current.bcf'.format(zipped_VCF) )
            commands_list.append( 'bcftools index current.bcf' )
            commands_list.append( 'bcftools view current.bcf --regions-file /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/ref_positions_to_select.txt -o {0}_ref_positions_for_genotyping.vcf -O v'.format(tag) )
            commands_list.append( 'rm current.bcf.csi current.bcf' )
            
        except KeyError: #applies only to the last job if it has < N isolates
            continue
        
    bash_scripts.append(commands_list)<jupyter_output><empty_output><jupyter_text>Submit each job to O2<jupyter_code>job_num = 1
for job_i_commands_list in bash_scripts:
    
    #append all commands in a single string to be submitted as a job
    VCF_extract_job = ''
    for command_i in job_i_commands_list:
        VCF_extract_job  = VCF_extract_job  + '\n' + command_i
    
    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/O2_SLURM_logs/VCF_refpos_for_genotyping_extraction/')

    job_name = 'VCF_' + str(job_num)

    s = Slurm(job_name , {'partition':'short' , 'n':'1' , 't':'0-12:00:00' , 'mem-per-cpu':'2G' , 'mail-type':'FAIL' , 'mail-user':'<EMAIL>'})

    #submits the job
    job_id = s.run(VCF_extract_job)

    print job_name
    job_num += 1<jupyter_output>submitted: Submitted batch job 3318878<jupyter_text>Run any remaining isolates through VCF subsetting procedure that don't have a corresponding subset VCF<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np
import itertools
import shutil
isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl' )<jupyter_output><empty_output><jupyter_text>__Check Progress of Jobs:__ Isolates that need to have Full VCF >> Subset VCF for Reference Positions of Interest<jupyter_code>isolate_ids_to_run_VCF_extraction = []
isolate_ids_completed = []

for tag in isolate_id_series.values:

    #create directory to store subset VCF and genotypes for each isolate
    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/' + tag
    
    #current.bcf file
    current_bcf_path = isolate_genotype_dir + '/current.bcf'
    
    #current.bcf.csi file
    current_bcf_csi_path = isolate_genotype_dir + '/current.bcf.csi'
    
    #subset VCF file
    subset_VCF_path = isolate_genotype_dir + '/' + tag + '_ref_positions_for_genotyping.vcf'
    
    #check for the existence of a VCF file, if TRUE job has either started or completed
    if os.path.exists(subset_VCF_path):
        
        #check to see if intermediary files exist, if they do then VCF subset task hasn't completed
        if os.path.exists(current_bcf_path) or os.path.exists(current_bcf_csi_path):
            
            #remove files to start task over
            ## os.remove(current_bcf_path)
            ## os.remove(current_bcf_csi_path)
            ## os.remove(subset_VCF_path)
            
            isolate_ids_to_run_VCF_extraction.append(tag)
            
        #VCF exists, intermediary files gone so task finished
        else: 
            
            isolate_ids_completed.append(tag)
            
    #VCF file doesn't exist, task never started
    else:
        
        isolate_ids_to_run_VCF_extraction.append(tag)
len(isolate_ids_completed)
len(isolate_ids_to_run_VCF_extraction)<jupyter_output><empty_output><jupyter_text>Subset VCF to reference positions of interest for isolates that had Full VCF files<jupyter_code>isolate_id_series = pd.Series(isolate_ids_to_run_VCF_extraction)<jupyter_output><empty_output><jupyter_text>Create BASH scripts to extract VCFs of N isolates per job<jupyter_code>N_isolates_per_job = 65

njobs = int( np.ceil( float( len( isolate_id_series ) ) / float(N_isolates_per_job) ) ) #number of jobs required if we split for every N isolates
print njobs
bash_scripts = []

#genomic data directory
rolling_DB_dir = '/n/data1/hms/dbmi/farhat/rollingDB/genomic_data/'

#split jobs up into tasks of N
for first_isolate_i_for_job in np.arange(0 , (njobs)*N_isolates_per_job , N_isolates_per_job):

    last_isolate_i_for_job = first_isolate_i_for_job + N_isolates_per_job
    
    commands_list = []
    commands_list.append( 'module load bcftools/1.9' ) #only need to load this once per job

    for isolate_i in range(first_isolate_i_for_job , last_isolate_i_for_job): 

        try:
            tag = isolate_id_series[isolate_i]
            zipped_VCF = rolling_DB_dir + tag + '/pilon/' + tag + '_full.vcf.gz'
            isolate_genotypes_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/' + tag

            commands_list.append( 'cd {}'.format(isolate_genotypes_dir) )
            commands_list.append( 'time bcftools view {} -O b -o current.bcf'.format(zipped_VCF) )
            commands_list.append( 'bcftools index current.bcf' )
            commands_list.append( 'bcftools view current.bcf --regions-file /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/ref_positions_to_select.txt -o {0}_ref_positions_for_genotyping.vcf -O v'.format(tag) )
            commands_list.append( 'rm current.bcf.csi current.bcf' )
            
        except KeyError: #applies only to the last job if it has < N isolates
            continue
        
    bash_scripts.append(commands_list)<jupyter_output><empty_output><jupyter_text>Submit each job to O2<jupyter_code>job_num = 1
for job_i_commands_list in bash_scripts:
    
    #append all commands in a single string to be submitted as a job
    VCF_extract_job = ''
    for command_i in job_i_commands_list:
        VCF_extract_job  = VCF_extract_job  + '\n' + command_i
    
    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/O2_SLURM_logs/VCF_refpos_for_genotyping_extraction/')

    job_name = 'VCF_' + str(job_num)

    s = Slurm(job_name , {'partition':'short' , 'n':'1' , 't':'0-12:00:00' , 'mem-per-cpu':'2G' , 'mail-type':'FAIL' , 'mail-user':'<EMAIL>'})

    #submits the job
    job_id = s.run(VCF_extract_job)

    print job_name
    job_num += 1<jupyter_output>submitted: Submitted batch job 3589690
submitted: Submitted batch job 3589692<jupyter_text>[7] This section of the notebook prepares the script that iterates over all of the VCF files in rollingDB and outputs a genotype vector.<jupyter_code>import vcf
import os
import pandas as pd
import numpy as np
isolate_ID = 'Peru4988' #INPUT - isolate ID<jupyter_output><empty_output><jupyter_text>Load SNV-Index annotation (for indexing genotypes)<jupyter_code>SNV_annotation_DF = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_SNP_annotation.pkl' )
SNV_annotation_DF.head()<jupyter_output><empty_output><jupyter_text>Convert SNV annotation into dictionaries with keys = ROW indices of genotype array<jupyter_code>#keys: H37Rv Reference Position , values: Genotype Row Index
ref_pos_dict = dict( zip( list(SNV_annotation_DF.pos) , list(SNV_annotation_DF.index) ) )<jupyter_output><empty_output><jupyter_text>Construct 'skeleton' Numpy Array that will store the supported base for all Reference Positions of interest for a single isolate in rollingDB- **0** for _A_- **1** for _C_- **2** for _G_- **3** for _T_- **9** for Bad Quality SNV calls<jupyter_code>#number of Reference Positions we're interested in (where High Quality SNV calls were made)
number_of_SNVs = len(ref_pos_dict.keys()) 
number_of_SNVs<jupyter_output><empty_output><jupyter_text>In constructing the Genotype Vector, we're going to iterate through all Reference Positions of interest and find the supported base from the VCF file for each position<jupyter_code>#create an empty numpy array to store genotypes
genotypes_array = np.full(number_of_SNVs , 9 , dtype = int)
genotypes_array
np.shape(genotypes_array)<jupyter_output><empty_output><jupyter_text>Iterate through the subset VCF file corresponding to the isolate in rollingDB and fill in genotypes vector<jupyter_code>#genomic data directory
VCF_genotypes_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/'

#directory that stores files for each sequenced isolate
directory_for_sequenced_isolate = VCF_genotypes_dir + isolate_ID + '/'

#position-base dictionaries (order of output in Pilon)
base_order = {0:'A' , 1:'C' , 2:'G', 3:'T'}
base_order_r = {'A':0 , 'C':1 , 'G':2 , 'T':3}

#load in VCF file if one was found
vcf_reader = vcf.Reader( open(directory_for_sequenced_isolate + isolate_ID + '_ref_positions_for_genotyping.vcf' , 'r') )

#iterate through each Variant Call 
for record in vcf_reader:
    
    #check to see if variant is SNP (length of Reference Allele is 1 base and there's only 1 alternate allele and length of alternate allele is 1 base is alternate allele exists)
    if (len(record.REF) == 1) and (len(record.ALT) == 1) and ( (record.ALT[0] == None) or (len(str(record.ALT[0])) == 1) ):

        #check to see if the call is a PASS by Pilon
        if record.FILTER == []:

            ##### Retrieve Relevant information for filtering quality of Base Call #####
            # Mean Base Quality @ locus
            BQ = record.INFO['BQ']
            # Mean Mapping Quality @ locus
            MQ = record.INFO['MQ']
            # Number of Reads w/ Deletion 
            DC = record.INFO['DC']
            # Number of Reads w/ Insertion
            IC = record.INFO['IC']
            # Depth of Valid Reads in Pileup
            VD = record.INFO['DP']

            ### Filtering Criteria
            #---> Mean Base Quality > 20
            #---> Mean Mapping Quality > 30
            #---> No Reads Supporting Insertions
            #---> No Reads Supporting Deletions
            #---> Number of High Quality Reads >= 25
            if (BQ > 20) and (MQ > 30) and (DC == 0) and (IC == 0) and (VD >= 25): 

                # Quality-Percentage - Percentage of As, Cs, Gs, Ts weighted by Q & MQ at locus
                QP = record.INFO['QP']

                if (sum(np.array(QP) >= 90) == 1): #1 Allele is supported by at least 90% of the reads, store information for this Base Call

                    #create a Boolean list for the Bases that have at least 90% of support from reads
                    supported_bases_boolean = list( np.array( QP ) >= 90 )

                    #position of Base [A,C,G,T] with at least 90% of support from reads
                    supported_bases_i = [i for i,x in enumerate(supported_bases_boolean) if x == True] #returns list of length 
                    supported_base = base_order[supported_bases_i[0]]

                    #write an exception in case a SNV (with a reference position that we're not interested in) 'sneaks' past the subsetting and filtering
                    try:
                        #store a '0', '1', '2', or '3' in Genotype Vector (Reference Position Index) referencing the Base with at least 90% support from the reads
                        genotypes_array[ref_pos_dict[record.POS]] = base_order_r[supported_base]
                    except KeyError:
                        pass<jupyter_output><empty_output><jupyter_text>Pickle Genotypes Array for downstream analysis<jupyter_code>genotypes_array
np.save(directory_for_sequenced_isolate + '/' + isolate_ID + '_genotypes' , genotypes_array , allow_pickle = True)<jupyter_output><empty_output><jupyter_text>[8] This section of the notebook launches the script that was created in the previous section of the notebook and submits jobs for batches of isolates to make genotypes arrays for each isolate.<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np
isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl' )
np.shape(isolate_id_series)<jupyter_output><empty_output><jupyter_text>Extract genotypes for each isolate from VCF file containing a subset of Reference Positions from the Full VCF files Create BASH scripts to extract genotypes of N isolates per job<jupyter_code>N_isolates_per_job = 720

njobs = int( np.ceil( float( len( isolate_id_series ) ) / float(N_isolates_per_job) ) ) #number of jobs required if we split for every N isolates
print njobs
bash_scripts = []

#split jobs up into tasks of N
for first_isolate_i_for_job in np.arange(0 , (njobs)*N_isolates_per_job , N_isolates_per_job):

    last_isolate_i_for_job = first_isolate_i_for_job + N_isolates_per_job
    
    commands_list = []

    for isolate_i in range(first_isolate_i_for_job , last_isolate_i_for_job): 

        try:
            tag = isolate_id_series[isolate_i]

            commands_list.append( 'python /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/python_scripts/rollingDB_extract_HQ_SNV_genotypes.py {}'.format(tag) )
            
        except KeyError: #applies only to the last job if it has < N isolates
            continue
        
    bash_scripts.append(commands_list)<jupyter_output><empty_output><jupyter_text>Submit each job to O2<jupyter_code>job_num = 1
for job_i_commands_list in bash_scripts:
    
    #append all commands in a single string to be submitted as a job
    genotypes_extract_job = ''
    for command_i in job_i_commands_list:
        genotypes_extract_job  = genotypes_extract_job  + '\n' + command_i
    
    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/O2_SLURM_logs/make_genotypes_arrays_from_VCF/')

    job_name = 'Geno' + str(job_num)

    s = Slurm(job_name , {'partition':'short' , 'n':'1' , 't':'0-6:00:00' , 'mem-per-cpu':'2G' , 'mail-type':'FAIL' , 'mail-user':'<EMAIL>'})

    #submits the job
    job_id = s.run(genotypes_extract_job)

    print job_name  + ' : ' +  str(job_id)
    job_num += 1<jupyter_output>submitted: Submitted batch job 3793777
submitted: Submitted batch job 3793778<jupyter_text>[9] Inspect the distribution of missing values after Genotype Matrix is constructed (from genotype vectors) and filter out SNPs & Isolates<jupyter_code>%matplotlib inline
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.ticker as ticker
from itertools import compress
from matplotlib.colors import LogNorm
from matplotlib import gridspec
import pickle<jupyter_output><empty_output><jupyter_text>Make sure all isolates have a corresponding _genotypes_ array<jupyter_code>isolate_id_series = pd.read_pickle( '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl' )
np.shape(isolate_id_series)<jupyter_output><empty_output><jupyter_text>__Check for Genotypes Arrays:__ Get isolate IDs that have a corresponding genotypes array<jupyter_code>isolate_ids_with_genotypes = []
isolate_ids_without_genotypes = []

for tag in isolate_id_series.values:

    #create directory to store subset VCF and genotypes for each isolate
    isolate_genotype_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/' + tag
    
    #genotypes_array
    genotypes_array_pickle_file = isolate_genotype_dir + '/' + tag + '_genotypes.npy'
    
    #check for the existence of a genotypes file for this isolate
    if os.path.exists(genotypes_array_pickle_file):
            
        isolate_ids_with_genotypes.append(tag)
            
    else:
        
        isolate_ids_without_genotypes.append(tag)
len(isolate_ids_with_genotypes)
len(isolate_ids_without_genotypes)
del isolate_ids_with_genotypes #clear variable from memory<jupyter_output><empty_output><jupyter_text>Concatenate genotypes arrays into a genotypes matrix for all isolates in rollingDB Store the __genotypes array__ for each isolate into a list<jupyter_code>genotypes_directory = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/Genotypes_per_isolate/'
genotype_arrays_per_isolate_list = []

#iterate through the genotypes arrays for each isolate
isolate_i = 0
for isolate_ID in isolate_id_series:
    
    #load genotypes array
    genotypes_array = np.load(genotypes_directory + isolate_ID + '/' + isolate_ID + '_genotypes.npy')
    
    #store genotypes array for this isolate in a list with all other genotypes arrays
    genotype_arrays_per_isolate_list.append(genotypes_array)
    
    if isolate_i % 1000 == 0:
        print float(isolate_i) / float(len(isolate_id_series))
    isolate_i += 1<jupyter_output>0.0
0.0295220382015
0.059044076403
0.0885661146046
0.118088152806
0.147610191008
0.177132229209
0.206654267411
0.236176305612
0.265698343814
0.295220382015
0.324742420217
0.354264458418
0.38378649662
0.413308534821
0.442830573023
0.472352611224
0.501874649426
0.531396687627
0.560918725829
0.59044076403
0.619962802232
0.649484840433
0.679006878635
0.708528916836
0.738050955038
0.767572993239
0.797095031441
0.826617069642
0.856139107844
0.885661146046
0.915183184247
0.944705222449
0.97422726065<jupyter_text>Stack genotype arrays for all isolates to create one __genotypes matrix__ for all isolates<jupyter_code>genotypes_matrix = np.stack(genotype_arrays_per_isolate_list , axis = -1)
np.shape(genotypes_matrix)
genotypes_matrix<jupyter_output><empty_output><jupyter_text>Load __SNP annotation file__<jupyter_code>SNP_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_SNP_annotation.pkl')
SNP_annotation_DF.head()
np.shape(SNP_annotation_DF)
del genotype_arrays_per_isolate_list #clear variable from memory<jupyter_output><empty_output><jupyter_text>[1] Drop SNP sites that correspond that (1) are located in Mobile Genetic Element (MGE) regions and/or (2) have low Empirical Base Recall (EBR) Scores   (Filter 1)  *FILTER* out SNPs with low EBR scores  Load list of low EBR score Reference positions to Drop (IL mapping/variant calling: sites with empirical score < 0.80)<jupyter_code>with open('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/H37Rv_sites_with_low_EBR_to_drop.pkl', 'rb') as f:
    H37Rv_positions_to_drop = pickle.load(f)
    
#convert to set, faster to check
H37Rv_positions_to_drop = set(H37Rv_positions_to_drop)<jupyter_output><empty_output><jupyter_text>Filter out sites with low EBR scores Drop unique SNPs from the __GENOTYPE MATRIX__ and __SNP ANNOTATION FILE__<jupyter_code>#construct boolean filter
SNPs_to_keep_filter = []
for H37Rv_pos_i in SNP_annotation_DF.pos:
    
    if H37Rv_pos_i not in H37Rv_positions_to_drop:
        
        SNPs_to_keep_filter.append(True)
        
    elif H37Rv_pos_i in H37Rv_positions_to_drop:
        
        SNPs_to_keep_filter.append(False)
#filter Genotypes Matrix
genotypes_matrix = genotypes_matrix[SNPs_to_keep_filter , :]

#filter SNP annotation file
SNP_annotation_DF = SNP_annotation_DF[SNPs_to_keep_filter]
SNP_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)
len(SNPs_to_keep_filter) - sum(SNPs_to_keep_filter) #number of sites dropped by this filter
np.shape(genotypes_matrix)
SNP_annotation_DF.head()
np.shape(SNP_annotation_DF) #Annotation for Rows of Genotypes Matrix<jupyter_output><empty_output><jupyter_text>(Filter 2) *FILTER* out SNPs in MGE regions  Collect list of MGE Reference Positions to Drop<jupyter_code>#csv files with category for each gene
gene_categories_df = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories.csv')
gene_categories_df.head()
#reference genome annotation
reference_genome_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt', '\t').set_index('name')<jupyter_output><empty_output><jupyter_text>Collect positions for **MGEs**<jupyter_code>H37Rv_positions_MGE_genes = []

for gene_id in gene_categories_df[gene_categories_df.Gene_Category == 'Mobile Genetic Element'].gene_id:
    
    #get reference positions tht corresond to each gene
    gene_id_info = reference_genome_annotation.loc[gene_id, :]
    gene_id = gene_id_info.name

    chrom_start = gene_id_info.chromStart
    chrom_end = gene_id_info.chromEnd

    #find the position of the first base relative to H37Rv in 5' -> 3'
    H37Rv_start = min(chrom_start , chrom_end)
    H37Rv_end = max(chrom_start , chrom_end)
    
    #store all corresponding H37Rv Reference Positions to gene_id in list of reference positions to scan
    H37Rv_positions_MGE_genes.append([H37Rv_start+1 , H37Rv_end+1])
len(H37Rv_positions_MGE_genes)<jupyter_output><empty_output><jupyter_text>FILTER out sites in MGE regions<jupyter_code>#construct boolean filter
SNPs_to_keep_filter = []
for H37Rv_pos_i in SNP_annotation_DF.pos:

    bad_pos_i = False

    for problem_region in H37Rv_positions_MGE_genes:

        if (H37Rv_pos_i > problem_region[0]) and (H37Rv_pos_i < problem_region[1]):
            bad_pos_i = True
            break

    if bad_pos_i == False:
        SNPs_to_keep_filter.append(True)
    else:
        SNPs_to_keep_filter.append(False)
#filter Genotypes Matrix
genotypes_matrix = genotypes_matrix[SNPs_to_keep_filter , :]

#filter SNP annotation file
SNP_annotation_DF = SNP_annotation_DF[SNPs_to_keep_filter]
SNP_annotation_DF.reset_index(drop = True , inplace = True) #re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)
len(SNPs_to_keep_filter) - sum(SNPs_to_keep_filter) #number of sites dropped by this filter
np.shape(genotypes_matrix)
SNP_annotation_DF.head()
np.shape(SNP_annotation_DF) #Annotation for Rows of Genotypes Matrix<jupyter_output><empty_output><jupyter_text>(Filter 3) *FILTER* out SNPs that are missing in $> 25\%$ of isolates  Investigate the Number of Missing Calls for each Unique SNP<jupyter_code>plt.style.use('ggplot')
plt.rcParams['lines.linewidth']=3
plt.rcParams['axes.facecolor']='1.0'
plt.rcParams['xtick.color']='black'
plt.rcParams['axes.grid']=True
plt.rcParams['axes.edgecolor']='black'
plt.rcParams['grid.color']= '1.0'
plt.rcParams.update({'font.size': 12})<jupyter_output><empty_output><jupyter_text>Look for SNPs with an abundance of missing values, drop SNPs that have missing values $> 25\%$ of the isolates<jupyter_code>missing_values = genotypes_matrix == 9
missing_values_per_SNP = missing_values.sum(axis = 1) #will range from 0 - 33,873
#impose a threshold for maximum number of missing values per SNP (across all isolates)
max_proportion_missing = 0.25
number_of_isolates = float(np.shape(genotypes_matrix)[1])

max_number_values_missing_per_SNP = np.floor( max_proportion_missing * number_of_isolates )
print max_number_values_missing_per_SNP #that we're going to allow
fig , ax = plt.subplots()

n , bins , patches = plt.hist( missing_values_per_SNP / number_of_isolates , bins = 50 , color = 'k' , rwidth = 1.0 , edgecolor='white', linewidth=1.0)
plt.hist( missing_values_per_SNP[missing_values_per_SNP > max_number_values_missing_per_SNP] / number_of_isolates , bins = bins , color = 'red' , rwidth = 1.0 , edgecolor='white', linewidth=1.0)

#plt.title('Mean Coverage across sequenced isolates in rollingDB', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.xlabel('Proportion of Isolates with Missing Values', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.ylabel('Number of Unique SNPs', fontweight = 'bold' , fontsize = 12, color = 'k')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.tick_params(labelcolor = 'k')
ax.set_yscale("log", nonposy='clip')

fig = plt.gcf()
fig.set_size_inches(6.5, 5.25)
fig.tight_layout()

plt.tick_params(axis='y', which='major', labelsize=12 , labelcolor = 'k')
plt.tick_params(axis='x', which='major', labelsize=12 , labelcolor = 'k')

#file_name = '/home/rv76/Farhat_Lab/Figures/Final_5/Scraping_rollingDB/Mean_Coverage_Distribution_across_all_isolates.pdf'
#plt.savefig(file_name, bbox_inches='tight', dpi = 300 , transparent = True)

plt.show()
print 'Number of Isolates = ' + str( number_of_isolates )<jupyter_output>Number of Isolates = 33873.0<jupyter_text>Range and Average Number of Missing Values per SNP<jupyter_code>np.min(missing_values_per_SNP)
np.mean(missing_values_per_SNP)
np.max(missing_values_per_SNP) 
print 'SNP with most missing values is missing from ' + str( float(np.max(missing_values_per_SNP) ) / number_of_isolates * 100) + '% of isolates'<jupyter_output>SNP with most missing values is missing from 99.7992501402% of isolates<jupyter_text>Drop SNPs that are flagged as 'missing values' in more than 25% of isolates in rollingDB<jupyter_code>print sum( missing_values_per_SNP > max_number_values_missing_per_SNP )<jupyter_output>3220<jupyter_text>Percentage of SNPs with many missing values that we're going to drop from analysis<jupyter_code>number_of_unique_SNPs = float(np.shape(genotypes_matrix)[0])
print 100.0 * ( sum( missing_values_per_SNP > max_number_values_missing_per_SNP ) / number_of_unique_SNPs )<jupyter_output>2.251480593216191<jupyter_text>FILTER out sites with missing calls in $> 25\%$ of isolates Drop unique SNPs with a high abundance of missing values<jupyter_code>np.shape( genotypes_matrix )<jupyter_output><empty_output><jupyter_text>Drop unique SNPs from the __GENOTYPE MATRIX__<jupyter_code>len( missing_values_per_SNP <= max_number_values_missing_per_SNP )
SNPs_to_keep_filter = missing_values_per_SNP <= max_number_values_missing_per_SNP
sum(SNPs_to_keep_filter)
genotypes_matrix = genotypes_matrix[SNPs_to_keep_filter , :]
np.shape(genotypes_matrix) #Genotype Matrix<jupyter_output><empty_output><jupyter_text>Drop unique SNPs with a high abundance of missing values from the __SNP ANNOTATION FILE__<jupyter_code>SNP_annotation_DF = SNP_annotation_DF[SNPs_to_keep_filter]

#re-index new filtered SNP annotation DF (so new index matches indexing of genotypes matrix rows)
SNP_annotation_DF.reset_index(drop = True , inplace = True)
SNP_annotation_DF.head()
np.shape(SNP_annotation_DF) #Annotation for Rows of Genotypes Matrix<jupyter_output><empty_output><jupyter_text>[9.2] Drop low-quality isolates  (Filter 4) *FILTER* out isolates that have missing base calls $> 25\%$ of SNP sites  Investigate the Number of Missing Calls for each Sequenced Isolate Look for isolates with an abundance of missing values, drop isolates that have missing values > 25% of SNP calls<jupyter_code>missing_values = genotypes_matrix == 9
missing_values_per_isolate = missing_values.sum(axis = 0) #will range from 0 - 139,797
#impose a threshold for maximum number of missing values per isolate (across all SNPs)
max_proportion_missing = 0.25
number_of_SNPs = float(np.shape(genotypes_matrix)[0])

max_number_values_missing_per_isolate = np.floor( max_proportion_missing * number_of_SNPs )
print max_number_values_missing_per_isolate #that we're going to allow
fig , ax = plt.subplots()

n , bins , patches = plt.hist( missing_values_per_isolate / number_of_SNPs , bins = 50 , color = 'k' , rwidth = 1.0 , edgecolor='white', linewidth=1.0)
plt.hist( missing_values_per_isolate[missing_values_per_isolate > max_number_values_missing_per_isolate] / number_of_SNPs , bins = bins , color = 'red' , rwidth = 1.0 , edgecolor='white', linewidth=1.0)


#plt.title('Mean Coverage across sequenced isolates in rollingDB', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.xlabel('Proportion of SNPs with Missing Values', fontweight = 'bold' , fontsize = 12, color = 'k')
plt.ylabel('Number of Isolates', fontweight = 'bold' , fontsize = 12, color = 'k')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

ax.tick_params(labelcolor = 'k')
ax.set_yscale("log", nonposy='clip')

fig = plt.gcf()
fig.set_size_inches(7.5, 5.25)
fig.tight_layout()

plt.tick_params(axis='y', which='major', labelsize=12 , labelcolor = 'k')
plt.tick_params(axis='x', which='major', labelsize=12 , labelcolor = 'k')

#file_name = '/home/rv76/Farhat_Lab/Figures/Final_5/Scraping_rollingDB/Mean_Coverage_Distribution_across_all_isolates.pdf'
#plt.savefig(file_name, bbox_inches='tight', dpi = 300 , transparent = True)

plt.show()
print 'Number of unique SNPs = ' + str( number_of_SNPs )<jupyter_output>Number of unique SNPs = 139797.0<jupyter_text>Range and Average Number of Missing Values per unique SNP<jupyter_code>np.min(missing_values_per_isolate)
np.mean(missing_values_per_isolate)
np.max(missing_values_per_isolate) 
print 'Isolate with most missing values has ' + str( float(np.max(missing_values_per_isolate) ) / number_of_SNPs * 100) + '% missing values'<jupyter_output>Isolate with most missing values has 100.0% missing values<jupyter_text>Drop isolates that are flagged as 'missing values' in more than 25% of isolates in rollingDB<jupyter_code>print sum( missing_values_per_isolate > max_number_values_missing_per_isolate )<jupyter_output>1518<jupyter_text>Percentage of isolates with many missing values that we're going to drop from analysis<jupyter_code>number_of_unique_isolates = float(np.shape(genotypes_matrix)[1])
print 100.0 * ( sum( missing_values_per_isolate > max_number_values_missing_per_isolate ) / number_of_unique_isolates )<jupyter_output>4.481445398990346<jupyter_text>FILTER out low-quality isolates with $> 25\%$ missing SNP calls Drop unique isolates with a high abundance of missing values<jupyter_code>np.shape( genotypes_matrix )<jupyter_output><empty_output><jupyter_text>Drop isolates from the __GENOTYPE MATRIX__<jupyter_code>len( missing_values_per_isolate <= max_number_values_missing_per_isolate )
isolates_to_keep_filter = missing_values_per_isolate <= max_number_values_missing_per_isolate
sum(isolates_to_keep_filter)
genotypes_matrix = genotypes_matrix[: , isolates_to_keep_filter]
np.shape(genotypes_matrix) #Genotype Matrix<jupyter_output><empty_output><jupyter_text>Drop isolates with a high abundance of missing values from __ISOLATE ANNOTATION FILE__<jupyter_code>isolate_annotation_series = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes/genotypes_isolate_annotation.pkl')
isolate_annotation_series.head()
np.shape(isolate_annotation_series) #Annotation for Columns of Genotypes Matrix
isolate_annotation_series = isolate_annotation_series[isolates_to_keep_filter]

#re-index new filtered isolate annotation series (so new index matches indexing of genotypes matrix columns)
isolate_annotation_series.reset_index(drop = True , inplace = True)
isolate_annotation_series.head()
np.shape(isolate_annotation_series) #Annotation for Columns of Genotypes Matrix<jupyter_output><empty_output><jupyter_text>Save annotation files and genotypes matrix<jupyter_code>#save isolate annotation file (columns of Genotype Matrix)
isolate_annotation_series.to_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes_Filtered/genotypes_isolate_annotation.pkl')

#save SNP annotation file (rows of Genotype Matrix)
SNP_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes_Filtered/genotypes_SNP_annotation.pkl')

#save Genotypes Matrix
np.save('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes_Filtered/genotypes_matrix' , genotypes_matrix , allow_pickle = True)<jupyter_output><empty_output><jupyter_text>[9] Annotate SNPs w/ H37Rv Locus Info   Retrieve information for H37Rv Reference Genome<jupyter_code>####################################################################################################################################################################################
#load reference genome and reference annotation
reference_genome = '/home/rv76/Farhat_Lab/Reference_Seqs/H37Rv/h37rv.fasta'
for reference_genome in SeqIO.parse(reference_genome, "fasta"):
    reference_genome.seq.alphabet = IUPAC.unambiguous_dna

reference_genome_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/h37rv_genome_summary.txt', '\t').set_index('name')

####### Make a dictionary of all H37Rv reference positions and corresponding genes #######
ReferencePosition_Gene_mapping = {} #keys: H37Rv Reference Positions , values: gene_ids (may be multiple genes for some Reference Positions)

#store a list corresponding to every Reference Position (to store all genes that map to the Reference Position)
for H37Rv_RefPos in range(0 , len(reference_genome.seq) + 1):
    ReferencePosition_Gene_mapping[H37Rv_RefPos] = []

for gene_id_index in range(0 , len(reference_genome_annotation.index)):

    gene_id_info = reference_genome_annotation.ix[gene_id_index , :]
    gene_id = gene_id_info.name

    chrom_start = gene_id_info.chromStart
    chrom_end = gene_id_info.chromEnd

    #find the position of the first base relative to H37Rv in 5' -> 3'
    H37Rv_start = min(chrom_start , chrom_end)
    H37Rv_end = max(chrom_start , chrom_end)

    #store all corresponding H37Rv Reference Positions to gene_id in dictionary
    for H37Rv_RefPos in range(H37Rv_start+1 , H37Rv_end+1):

        ReferencePosition_Gene_mapping[H37Rv_RefPos].append(gene_id) #append gene_id to list already in dict
####################################################################################################################################################################################<jupyter_output>/home/rv76/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:18: DeprecationWarning: 
.ix is deprecated. Please use
.loc for label based indexing or
.iloc for positional indexing

See the documentation here:
http://pandas.pydata.org/pandas-docs/stable/indexing.html#ix-indexer-is-deprecated<jupyter_text>Annotate SNPs To annotate SNPs, we get gene_ids that the SNPs 'land in' and the gene-relative coordinate.<jupyter_code>#load SNP annotation file (rows of Genotype Matrix)
SNP_annotation_DF = pd.read_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes_Filtered/genotypes_SNP_annotation.pkl')

SNP_annotation_DF.head()
np.shape(SNP_annotation_DF)<jupyter_output><empty_output><jupyter_text>Gene Categories<jupyter_code>#get Gene Categories
gene_categories = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/gene_categories/gene_categories.csv').set_index('name')

#create dictionary, keys: H37Rv locus, values: gene category
gene_categories_map = dict( zip( list(gene_categories.gene_id) , list(gene_categories.Gene_Category) ) )<jupyter_output><empty_output><jupyter_text>Gene Names Parse through Mycobrowser (tuberculist) to collect gene names<jupyter_code>mycobrowser_table = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/H37Rv/Mycobrowser_Release3/Mycobacterium_tuberculosis_H37Rv_txt_v3.txt' , sep = '\t')

#retreive columns for locus tag & gene name
mycobrowser_table = mycobrowser_table.loc[: , ['Locus' , 'Name']]

#create dictionary, keys: H37Rv loucs, values: gene category
gene_names_map = dict( zip( list(mycobrowser_table.Locus) , list(mycobrowser_table.Name) ) )<jupyter_output><empty_output><jupyter_text>Iterate through each SNP to annotate Some SNPs are present in more than 1 gene, have to account for this<jupyter_code>#create dictionaries for each columsn of information
locus_tag_dict = {} #H37Rv id
genomic_coord_dict = {} #coordinate on gene, 5' -> 3'
gene_categ_dict = {} #gene Category
gene_name_dict = {} #gene Name

for SNP_index in SNP_annotation_DF.index:

    #reference position, ref allele and alt allele for SNP
    SNP_info = SNP_annotation_DF.loc[SNP_index , :]

    #find the gene that SNP occurs on
    ref_seq_position = SNP_info.pos

    #check list corresponding to H37Rv coordinate to see if there are any genes associated with RefPosition
    if len(ReferencePosition_Gene_mapping[ref_seq_position]) > 0:
        
        #create lists to store annotation information for SNP and associated gene(s)
        gene_intergenic_id_list = []
        gene_relative_coord_list = []
        gene_categ_list = []
        gene_name_list = []

        #iterate through all genes that ReferencePosition is mapped to (i.e. SNP might correspond to 2 genes)
        for gene_intergenic_id in ReferencePosition_Gene_mapping[ref_seq_position]:

            #find genomic coordinate of SNP relative to gene (subtract 1 since reference seq starts counting at 1)
            gene_relative_coord = (ref_seq_position - 1) - min( reference_genome_annotation.loc[gene_intergenic_id , 'chromStart'] , reference_genome_annotation.loc[gene_intergenic_id , 'chromEnd'] )

            gene_intergenic_id_list.append( gene_intergenic_id )
            gene_relative_coord_list.append( gene_relative_coord )
            
            #use gene id to get category and names from mappings
            if gene_intergenic_id in gene_categories_map.keys():
                gene_categ_list.append( gene_categories_map[gene_intergenic_id] )
            else:
                gene_categ_list.append( 'None' )
                
            if gene_intergenic_id in gene_names_map.keys():
                gene_name_list.append( gene_names_map[gene_intergenic_id] ) 
            else:
                gene_name_list.append( 'None' )

        #store all information for SNP in relevant dictionaries
        locus_tag_dict[SNP_index] = gene_intergenic_id_list
        genomic_coord_dict[SNP_index] = gene_relative_coord_list
        gene_categ_dict[SNP_index] = gene_categ_list
        gene_name_dict[SNP_index] = gene_name_list
            
    #H37Rv RefSeq position was not found in dictionary so Reference position must be intergenic    
    elif len(ReferencePosition_Gene_mapping[ref_seq_position]) == 0:

        #store all information for SNP in relevant dictionaries
        locus_tag_dict[SNP_index] = ['intergenic']
        genomic_coord_dict[SNP_index] = [np.nan]
        gene_categ_dict[SNP_index] = ['None']
        gene_name_dict[SNP_index] = ['None']

    #keep track of progress of loop
    num_SNPs = float( np.shape(SNP_annotation_DF)[0] )
    if SNP_index % np.ceil(0.20 * num_SNPs) == 0:
        print float(SNP_index) / num_SNPs<jupyter_output>0.0
0.200004291938
0.400008583875
0.600012875813
0.80001716775<jupyter_text>Update SNP annotation DF<jupyter_code>#append columns for info from annotation
SNP_annotation_DF['gene_id'] = pd.Series(locus_tag_dict)
SNP_annotation_DF['gene_pos'] = pd.Series(genomic_coord_dict)
SNP_annotation_DF['gene_category'] = pd.Series(gene_categ_dict)
SNP_annotation_DF['gene_name'] = pd.Series(gene_name_dict)
np.shape( SNP_annotation_DF )
SNP_annotation_DF.head()
SNP_annotation_DF.tail()
#update SNP annotation file (rows of Genotype Matrix) with extra SNP information
SNP_annotation_DF.to_pickle('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/rolling_DB_scrape/Genotypes_Filtered/genotypes_SNP_annotation_with_gene_info.pkl')<jupyter_output><empty_output><jupyter_text>How to investigate SNPs that appeared on more than 1 gene<jupyter_code>SNP_annotation_DF[ SNP_annotation_DF.gene_name.apply(len) > 1 ].head()<jupyter_output><empty_output><jupyter_text>How to investigate SNPs that appeared within a particular gene<jupyter_code>SNP_annotation_DF[['Rv0667' in gene_list for gene_list in SNP_annotation_DF.gene_id]].head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09237457732385619, 'fraction_numerical': 0.02248292043337244, 'mean_word_length': 4.052367338400391, 'pattern_counts': {'":': 0, '<': 317, '<?xml version=': 0, '>': 365, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 42, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '80', 'dir': 'jupyter-structured-clean-dedup', 'id': '541096', 'license': 'Naumen', 'n_tokens_mistral': 25712, 'n_tokens_neox': 23200, 'n_words': 5983, 'path': 'SNPs from Public Samples & tSNE Homoplasy Visualization/(A) Scraping Public WGS Samples for Constructing Genotypes Matrix.ipynb', 'repo_name': 'farhat-lab/in-host-Mtbc-dynamics', 'size': '72531'}