<jupyter_start><jupyter_text>This notebook tests the pipeline for calling differing Base Calls (Single Nucleotide Variant) between two samples from a pair of samples (*Longitudinal* or *Replicate* pair). This notebook also submits the jobs to run this variant extraction pipeline for each sample pair in our study.<jupyter_code>from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))<jupyter_output><empty_output><jupyter_text>Prepare the script<jupyter_code>%matplotlib inline

import vcf
import os
import pandas as pd
import numpy as np
from itertools import compress
import ast
import itertools
import time
import sys
import pickle
import shutil

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<jupyter_output><empty_output><jupyter_text>*Functions* to extract Base Calls (that have different alterante allele frequencies) from a pair of VCF files<jupyter_code>def Get_Base_Calls_from_VCF(parent_output_dir, tag):

    '''
    This function takes as input, the isolate tag, then extracts SINGLE base calls from
    the VCF file generated by Pilon. The function returns a DataFrame for a single isolate
    that contains [1] PASS calls that support the reference allele, [2] PASS calls that support
    an alternate allele and [3] AMB calls that support both the reference and alternate alleles.
    Note: Many base calls won't have an Alternate Allele.
    '''

    #directory where VCF file for sequenced isolated is stored & filename
    VCF_file = parent_output_dir + tag + '/pilon/' + tag + '.vcf'

    vcf_reader = vcf.Reader(open(VCF_file , 'r'))

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

    #create dictionaries to store information for each call
    ref_bases = {}
    qry_bases = {}
    ref_positions = {}
    quality_scores = {}
    SNP_type = {}
    PASS_filter = {}
    INFO_for_call = {}
    alt_allele_freq = {}
    valid_depth = {}

    #store the DEPTH of all base & variant calls (depth not stored for major structural variants)
    coverage_for_all_base_and_small_indel_calls = []

    #indexer for dataframe containing SNPs
    index = 0

    #iterate through each record from VCF file
    for record in vcf_reader:

        ######## store the DEPTH of all base & variant calls (depth not stored for major structural variants) ########
        try:
            coverage_for_all_base_and_small_indel_calls.append( record.INFO['DP'] )

        except KeyError: #major structural variant with no depth value
            pass 

        ######## Classify SNP-type according to Pilon Flags ########

        ### Filtering Criteria
        #---> Check that Filter is either PASS ([]) or AMBIGUOUS (['Amb']) , no 'DEL' or 'LOWCOV' flags from Pilon

        ### Call supporting REFERENCE
        #if there's no alternate allele and Filter is PASS; PASS call supporting Ref allele
        if (len(record.ALT) == 1) and (record.ALT == [None]) and (record.FILTER == []):

            alt_allele = 'Z' #no alternate allele supported by reads (might change through filtering process further down if small fraction of reads support an alternate)
            ref_allele = record.REF #allele on H37Rv
            SNP_type_for_Call = 'Ref_PASS'

        ### Call supporting ALTERNATE 
        #if there's 1 alternate allele, and Filter is PASS; PASS call supporting Alt allele
        elif (len(record.ALT) == 1) and (record.ALT != [None]) and (record.FILTER == []):

            alt_allele = str( record.ALT[0] ) #alternate allele supported by reads
            ref_allele = record.REF #allele on H37Rv
            SNP_type_for_Call = 'Alt_PASS'

        ### Call ambiguous between REFERNENCE and ALTERNATE
        #if there's 1 alternate allele, and Filter is AMB; AMB call supporting supporting neither Alt nor Ref
        elif (len(record.ALT) == 1) and (record.ALT != [None]) and (record.FILTER == ['Amb']):

            alt_allele = str( record.ALT[0] ) #alternate allele supported by reads
            ref_allele = record.REF #allele on H37Rv
            SNP_type_for_Call = 'Amb'

        ### Bad Quality Call, don't store information for this Base Call, move onto next
        else:

            continue

        ### Filtering Criteria (for Reference Confirmation or SNP)
        #---> Make sure that the length of all any alleles is 1 (No Variant or 1 Base) ; if variant, then variant is a SNP and not an INDEL
        if (len(ref_allele) == 1) and (len(alt_allele) == 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):

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

                ### Filtering Criteria
                #---> Make sure that Reads support at most 2 alleles (either confirm the Reference or provide evidence for an alternate allele)
                #-----> If there are 2 alleles (make sure 1 of them is the Reference Allele)

                if (sum(np.array(QP) == 0) == 3): #1 Allele supported by Reads (either Reference or an 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
                    
                    #Reference Allele is NOT supported by at least 1% of the reads (2 alleles supported but neither is the reference), don't store information for this Base Call
                    else:
                        continue

                #more than 2 alleles supported by reads, don't store information for this Base Call
                else: 
                    continue 

                #After extensive filtering and categorization, store all of the pertinent information about the Base Call
                ref_bases[index] = ref_allele
                qry_bases[index] = alt_allele
                ref_positions[index] = record.POS
                PASS_filter[index] = record.FILTER
                quality_scores[index] = record.QUAL
                SNP_type[index] = SNP_type_for_Call
                INFO_for_call[index] = record.INFO
                alt_allele_freq[index] = alt_allele_frequency
                valid_depth[index] = record.INFO['DP']

                index += 1

    #convert dictionaries to series
    ref_bases = pd.Series(ref_bases)
    qry_bases = pd.Series(qry_bases)
    ref_positions = pd.Series(ref_positions)
    quality_scores = pd.Series(quality_scores)
    SNP_type = pd.Series(SNP_type)
    PASS_filter = pd.Series(PASS_filter)
    INFO_for_call = pd.Series(INFO_for_call)
    alt_allele_freq = pd.Series(alt_allele_freq)
    valid_depth = pd.Series(valid_depth)
    tag_series = pd.Series([tag]*len(ref_bases.index) , index = ref_bases.index) #psuedo-series to store the tag id of the isolate

    #create DataFrame to hold all base calls for a given sample
    Variant_Call_DF = pd.DataFrame()
    Variant_Call_DF['ref_base'] = ref_bases
    Variant_Call_DF['alt_base'] = qry_bases
    Variant_Call_DF['ref_position'] = ref_positions
    Variant_Call_DF['quality'] = quality_scores
    Variant_Call_DF['SNP_type'] = SNP_type
    Variant_Call_DF['PASS_filter'] = PASS_filter
    Variant_Call_DF['INFO'] = INFO_for_call
    Variant_Call_DF['alt_AF'] = alt_allele_freq
    Variant_Call_DF['depth'] = valid_depth
    Variant_Call_DF['tag'] = tag_series

    ### Filtering Criteria (high coverage)
    #---> Filter out Base Calls with total # of reads aligning above the 99th percentile of covered positions in sample
    ninety_ninth_percentile_coverage = np.percentile( np.array( coverage_for_all_base_and_small_indel_calls ) , 99 )
    Variant_Call_DF = Variant_Call_DF[Variant_Call_DF.depth < ninety_ninth_percentile_coverage] #keep base calls that dont have unusually high coverage

    #reset index after dropping base calls with unusually high coverage
    Variant_Call_DF.reset_index(inplace = True , drop = True)

    #save REFERENCE POSITIONS only for filtered (good) single base call variants as a pickle in same directory as VCF file
    Variant_Call_DF.loc[: , 'ref_position'].to_pickle(parent_output_dir + tag + '/' + tag + '_ref_positions_for_filtered_base_calls.pkl')

    ##SIDE NOTE: use "pd.read_pickle(parent_output_dir + tag + '/' + tag + '_filtered_base_calls.pkl')" to unpickle series object

    return Variant_Call_DF #DataFrame for base calls for a single isolate
def Get_different_Base_Calls_between_isolates(parent_output_dir, tag_list):
    
    '''
    This function uses the previously defined function to filter the VCF files corresponding to multiple
    isolates from the same isolate pair. The input is a list of tags (fastq ids) for isolates obtained for a given isolate pair
    Once the VCF files have been seperately filtered for good calls (PASS supporting Reference Allele or PASS supporting 
    Alternate Allele or ALT supporting two Alleles; and no unusually high coverage) this function then compares the 
    single base calls between the multiple isolates. The output is a DataFrame that contains GOOD base calls
    for both isolates but for which the Alternate Allee Frequencies DIFFERED.
    '''

    #store VCF output for both isolates, key: tag , value: DF for variant output
    Variant_Call_DFs = {}

    #iterate through both isolates that came from subject
    for tag in tag_list:

        #call function to get variants for isolate and store in dict
        Variant_Call_DFs[tag] = Get_Base_Calls_from_VCF(parent_output_dir, tag)

    #concatenate DataFrames for each isolate into 1 DataFrame
    Variant_Call_Master_DF = pd.concat(Variant_Call_DFs.values())

    #reset index of master dataframe
    Variant_Call_Master_DF.reset_index(inplace = True, drop = True)

    #sort values in order of reference position
    Variant_Call_Master_DF.sort_values(by = 'ref_position' , ascending = True , inplace = True)

    ############################################################################################

    #create a list to store all of the indices that should be dropped
    base_calls_to_drop = [] #indices in dataframe to drop later

    #### FILTER methodology ####
    ## if a call for a reference position is PASS-REFERENCE with alt AF of 0.0 for both isolates, then drop base call (for each isolate) (all reads from both sequenced isolates support the same allele)
    PASS_REF_filter = [(Ref_PASS and zero_alt_AF) for Ref_PASS , zero_alt_AF in  zip( list( Variant_Call_Master_DF.SNP_type == 'Ref_PASS' ) , list( Variant_Call_Master_DF.alt_AF == 0.0 ) ) ]

    ## if a call for a reference position is PASS-ALTERNATE with alt AF of 1.0 for both isolates, then drop base call (for each isolate) (all reads from both sequenced isolates support the same allele)
    PASS_ALT_filter = [(Alt_PASS and one_alt_AF) for Alt_PASS , one_alt_AF in  zip( list( Variant_Call_Master_DF.SNP_type == 'Alt_PASS' ) , list( Variant_Call_Master_DF.alt_AF == 1.0 ) ) ]

    #subset to base calls that support criteria for each pair matching criteria
    PASS_REF_base_calls = Variant_Call_Master_DF[PASS_REF_filter]
    PASS_ALT_base_calls = Variant_Call_Master_DF[PASS_ALT_filter]

    #find matching base calls that support Reference allele entirely
    matching_PASS_REF_base_calls_filter = PASS_REF_base_calls.duplicated(subset = ['ref_base' , 'alt_base' , 'ref_position' , 'SNP_type' , 'alt_AF'] , keep = False)
    matching_base_call_indices_to_drop = list( PASS_REF_base_calls[matching_PASS_REF_base_calls_filter].index )
    base_calls_to_drop = base_calls_to_drop + matching_base_call_indices_to_drop

    #find matching base calls that support Alternate allele entirely
    matching_PASS_ALT_base_calls_filter = PASS_ALT_base_calls.duplicated(subset = ['ref_base' , 'alt_base' , 'ref_position' , 'SNP_type' , 'alt_AF'] , keep = False)
    matching_base_call_indices_to_drop = list( PASS_ALT_base_calls[matching_PASS_ALT_base_calls_filter].index )
    base_calls_to_drop = base_calls_to_drop + matching_base_call_indices_to_drop

    #drop base calls from Master variant DF that match between isolates AND support 1 allele completely
    Variant_Call_Master_DF.drop(labels = base_calls_to_drop , axis = 0 , inplace = True)

    ############################################################################################

    #### FILTER methodology ####
    ## from the remaining base calls, find base calls for which we have a reference position for BOTH isolates (calls in both samples were good quality)
    ## that is, we will have a PASS call or AMB call (w/o low coverage or deletion flags) for a reference position for BOTH isolates
    Variant_Call_Master_DF = Variant_Call_Master_DF[Variant_Call_Master_DF.duplicated(subset = 'ref_position' , keep = False)]

    #reset index of filtered SNP-variant Master DataFrame
    Variant_Call_Master_DF.reset_index(inplace = True, drop = True)
    
    ############################################################################################
    
    #### FILTER methodology ####
    ## filter out paired Base Calls that have a difference in Alternate Allele Frequency of less than 5%
    
    alt_AF_isolate_A = Variant_Call_Master_DF.loc[range(0 , np.shape(Variant_Call_Master_DF)[0] , 2) , 'alt_AF']
    alt_AF_isolate_B = Variant_Call_Master_DF.loc[range(1 , np.shape(Variant_Call_Master_DF)[0] , 2) , 'alt_AF']

    alt_AF_diff_btwn_paired_isolates = abs(alt_AF_isolate_A.values - alt_AF_isolate_B.values)

    isolate_A_Base_Call_indices_small_change_alt_AF = list(alt_AF_isolate_A[alt_AF_diff_btwn_paired_isolates < 0.05].index)
    isolate_B_Base_Call_indices_small_change_alt_AF = list(alt_AF_isolate_B[alt_AF_diff_btwn_paired_isolates < 0.05].index)

    Variant_Call_Indices_SMALL_Alt_AF_Diff = isolate_A_Base_Call_indices_small_change_alt_AF + isolate_B_Base_Call_indices_small_change_alt_AF

    #drop paired Base Calls w/ corresponding change in Alterante Allele Frequency < 5%
    Variant_Call_Master_DF.drop(Variant_Call_Indices_SMALL_Alt_AF_Diff , axis = 0 , inplace = True)
    
    #reset index of final filtered SNP-variant Master DataFrame
    Variant_Call_Master_DF.reset_index(inplace = True, drop = True)
    
    return Variant_Call_Master_DF #DataFrame for good base calls that differed between isolates from same subject or replicate pair<jupyter_output><empty_output><jupyter_text>Call Functions to call SNPs between a (*longitudinal* or *replicate*) pair of isolates  *For Testing* (longitudinal pairs)<jupyter_code>#INPUT: each sub-directory has the pilon-outputed VCF file
parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output/'

#INPUT: Sample annotation file path for longitudinal or replicate samples (already filtered with Kraken and F2)
sample_annotation_file_path = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Longitudinal_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv'

#INPUT: either subject ID (longitudinal) or isolate pair ID (replicate)
isolate_pair_ID = '1960'

#INPUT: parent directory where all SNPs
SNP_call_parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/longitudinal_SNPs/all_SNPs_between_longitudinal_pairs/'<jupyter_output><empty_output><jupyter_text>*For Testing* (replicate pairs)<jupyter_code>#INPUT: each sub-directory has the pilon-outputed VCF file
parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output_REPLICATES/'

#INPUT: Sample annotation file path for longitudinal or replicate samples (already filtered with Kraken and F2)
sample_annotation_file_path = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Replicate_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv'

#INPUT: either subject ID (longitudinal) or isolate pair ID (replicate)
isolate_pair_ID = 'replicate_pair_67'

#INPUT: parent directory where all SNPs
SNP_call_parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/replicate_SNPs/all_SNPs_between_replicate_pairs/'<jupyter_output><empty_output><jupyter_text>*For Script*<jupyter_code>#INPUT: each sub-directory has the pilon-outputed VCF file
parent_output_dir = sys.argv[1]

#INPUT: Sample annotation file path for longitudinal or replicate samples (already filtered with Kraken and F2)
sample_annotation_file_path = sys.argv[2]

#INPUT: either subject ID (longitudinal) or isolate pair ID (replicate)
isolate_pair_ID = sys.argv[3]

#INPUT: parent directory where all SNPs 
SNP_call_parent_output_dir = sys.argv[4] 
#Sample Annotation & Tag Information for Isolate Pair
##############################################################################################################################

#import annotation file for all isolates and corresponding tags
sample_annotation = pd.read_csv(sample_annotation_file_path , sep = ',').set_index('patient_id')

#subset to both of the isolates taken from the isolate pair
sequenced_isolates_annot_for_sample_pair = sample_annotation.loc[isolate_pair_ID , :] #annotation (only) for both isolates from the sample pair

#get list of tags for both longitudinal isolates from isolate pair
tags = []

#iterate across both isolates
for isolate_i in [0 , 1]:
    
    #subset to 1 of the isolates taken from the isolate pair
    sequenced_isolate_annot = sequenced_isolates_annot_for_sample_pair.iloc[isolate_i , :] #annotation for 1 of the isolates from isolate pair
    
    #get relevant information for specific isolate
    population = sequenced_isolate_annot.population
    isolate_pair_ID = sequenced_isolate_annot.name
    tag = sequenced_isolate_annot.tag
    
    #append tag to tags list to get Variants From and compare other other isolates in tag list
    tags.append(tag)
    
#make necessary directory for isolate pair (to store variation between longitudinal isolates)
SNP_call_output_dir = SNP_call_parent_output_dir + population + '_' + isolate_pair_ID
    
if os.path.exists(SNP_call_output_dir):
    shutil.rmtree(SNP_call_output_dir)
    os.makedirs(SNP_call_output_dir)
elif not os.path.exists(SNP_call_output_dir):
    os.makedirs(SNP_call_output_dir)
    
#BASE CALLS
##############################################################################################################################
    
#get 'good' calls that differ between isolates for the same isolate pair (input is the tags corresponding to the isolates from isolate pair)
Base_Call_Master_DF = Get_different_Base_Calls_between_isolates(parent_output_dir, tags)

#append psuedo-columns to note population & isolate pair ID
population_series = pd.Series( [population] * np.shape(Base_Call_Master_DF)[0] , index = Base_Call_Master_DF.index )
isolate_pair_ID_series = pd.Series( [isolate_pair_ID] * np.shape(Base_Call_Master_DF)[0] , index = Base_Call_Master_DF.index )

Base_Call_Master_DF['population'] = population_series
Base_Call_Master_DF['patient_id'] = isolate_pair_ID_series

##pickle CSV file containing base calls that differ between paired isolates
Base_Call_Master_DF.to_pickle(SNP_call_output_dir + '/' + 'base_calls_different_between_isolates.pkl')<jupyter_output><empty_output><jupyter_text>Submit Jobs - Longitudinal Isolates<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np<jupyter_output><empty_output><jupyter_text>Import annotation file for all sample pairs<jupyter_code>#annotation file for all isolates and corresponding tags
sample_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Longitudinal_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv' , sep = ',').set_index('patient_id')<jupyter_output><empty_output><jupyter_text>Find which subjects have 2 valid VCF files for both longitudinal isolates<jupyter_code>successful_JankyPipe_run = []

for isolate_i in range(0 , np.shape(sample_annotation)[0]):

    #get the tag ID for the fastq files (same as ID for fastq files)
    tag = sample_annotation.tag[isolate_i]

    #where pilon VCF and lineage information will be stored [LAB FOLDER]
    output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output/' + tag
    
    #check to see 'Lineage Call' folder exists in the output directory (last thing that is run in JankyPipe)
    if os.path.exists(output_dir + '/fast-lineage-caller/'):
        
        successful_JankyPipe_run.append('yes')
        
    else:
        
        successful_JankyPipe_run.append('no')
        
sample_annotation['successful_JankyPipe_run'] = successful_JankyPipe_run
sample_annotation[sample_annotation.successful_JankyPipe_run == 'no']
#subjects for which we do not have 2 good VCF files for
bad_isolate_pair_IDs = set( list( sample_annotation[sample_annotation.successful_JankyPipe_run == 'no'].index ) ) 
print list(bad_isolate_pair_IDs)<jupyter_output>[]<jupyter_text>Get list of subjects (longitudinal isolate pairs) to get SNP variants from<jupyter_code>isolate_pair_ID_list = list( set( list( sample_annotation.index ) ) - bad_isolate_pair_IDs )
print isolate_pair_ID_list
len(isolate_pair_ID_list)<jupyter_output><empty_output><jupyter_text>Submit a job for each subject to get the SNPs<jupyter_code>parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output/'
sample_annotation_file_path = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Longitudinal_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv'
SNP_call_parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/longitudinal_SNPs/all_SNPs_between_longitudinal_pairs/'

for isolate_pair_ID in isolate_pair_ID_list:
    
    longitudinal_SNP_variants_job = 'python /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/python_scripts/BaseCall_variant_between_isolate_pair_collection.py ' + parent_output_dir + ' ' + sample_annotation_file_path + ' ' + isolate_pair_ID + ' ' + SNP_call_parent_output_dir

    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/SLURM_logs/BaseCalls_between_isolate_pairs_SLURM_logs/')

    job_name = isolate_pair_ID

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

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

    print job_name  + ' : ' +  str(job_id)<jupyter_output>submitted: Submitted batch job 2663382<jupyter_text>Submit Jobs - Replicate Isolates<jupyter_code>import os
import pandas as pd
from slurmpy import Slurm
import numpy as np<jupyter_output><empty_output><jupyter_text>Import annotation file for all sample pairs<jupyter_code>#annotation file for all isolates and corresponding tags
sample_annotation = pd.read_csv('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Replicate_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv' , sep = ',').set_index('patient_id')<jupyter_output><empty_output><jupyter_text>Find which replicate pairs have 2 valid VCF files for both replicate isolates<jupyter_code>successful_JankyPipe_run = []

for isolate_i in range(0 , np.shape(sample_annotation)[0]):

    #get the tag ID for the fastq files (same as ID for fastq files)
    tag = sample_annotation.tag[isolate_i]

    #where pilon VCF and lineage information will be stored [LAB FOLDER]
    output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output_REPLICATES/' + tag
    
    #check to see 'Lineage Call' folder exists in the output directory (last thing that is run in JankyPipe)
    if os.path.exists(output_dir + '/fast-lineage-caller/'):
        
        successful_JankyPipe_run.append('yes')
        
    else:
        
        successful_JankyPipe_run.append('no')
        
sample_annotation['successful_JankyPipe_run'] = successful_JankyPipe_run
sample_annotation[sample_annotation.successful_JankyPipe_run == 'no']
#subjects for which we do not have 2 good VCF files for
bad_isolate_pair_IDs = set( list( sample_annotation[sample_annotation.successful_JankyPipe_run == 'no'].index ) ) 
print list(bad_isolate_pair_IDs)<jupyter_output>[]<jupyter_text>Get list of subjects (longitudinal isolate pairs) to get SNP variants from<jupyter_code>isolate_pair_ID_list = list( set( list( sample_annotation.index ) ) - bad_isolate_pair_IDs )
print isolate_pair_ID_list
len(isolate_pair_ID_list)<jupyter_output><empty_output><jupyter_text>Submit a job for each subject to get the SNPs<jupyter_code>parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/JankyPipe/output_REPLICATES/'
sample_annotation_file_path = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/CSV_files/sample_annotation_files/Replicate_fastq_path_names_and_JankyPipe_tags_filtered_with_Kraken_and_F2.csv'
SNP_call_parent_output_dir = '/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/pickled_files/variant_calling/replicate_SNPs/all_SNPs_between_replicate_pairs/'

for isolate_pair_ID in isolate_pair_ID_list:
    
    replicate_SNP_variants_job = 'python /n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/python_scripts/BaseCall_variant_between_isolate_pair_collection.py ' + parent_output_dir + ' ' + sample_annotation_file_path + ' ' + isolate_pair_ID + ' ' + SNP_call_parent_output_dir

    #directory where you want output + error files
    os.chdir('/n/data1/hms/dbmi/farhat/Roger/inhost_TB_dynamics_project/SLURM_logs/BaseCalls_between_isolate_pairs_SLURM_logs/')

    job_name = isolate_pair_ID

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

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

    print job_name  + ' : ' +  str(job_id)<jupyter_output>submitted: Submitted batch job 2663672
{'fraction_non_alphanumeric': 0.09437086092715231, 'fraction_numerical': 0.006690093255845384, 'mean_word_length': 3.994431319608505, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 83, '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': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '475547', 'license': 'Naumen', 'n_tokens_mistral': 9670, 'n_tokens_neox': 8943, 'n_words': 2487, 'path': '(E) Process Replicate and Longitudinal Samples – Collect SNV Calls with Differing Alt AF between Sample Pairs.ipynb', 'repo_name': 'farhat-lab/in-host-Mtbc-dynamics', 'size': '29628'}