<jupyter_start><jupyter_text>Exploratory Data AnalysisFirst pass at analysing the output of blasting AML RNA-Seq data against the three tryptase sequences from Jonathon.<jupyter_code>import pysam
import glob
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
def get_files(folder):
    
    return glob.glob('{folder}/*.filtered.bam'.format(folder=folder))
list_of_aml_results = get_files('../aml_results/')
list_of_control_results = get_files('../control_results/')
aml_file_names = [x.split('/')[len(x.split('/'))-1] for x in list_of_aml_results]
control_file_names = [x.split('/')[len(x.split('/'))-1] for x in list_of_control_results]
df_aml = pd.DataFrame(index=aml_file_names)
df_control = pd.DataFrame(index=control_file_names)
df_aml['file_location'] = list_of_aml_results
df_control['file_location'] = list_of_control_results
df_aml['source'] = 'AML'
df_control['source'] = 'control'
df_aml.head()
df_control.head()
df = df_aml.append(df_control)
df.count()
def get_read_count(df):

    return int(pysam.view("-c", df['file_location']))
df['alignment_count'] = df.apply(get_read_count, axis=1)
df.head()
def get_transcript_count(df, transcript_name):
    
    """
    Return the count of the transcript i.e how many of each of the three are in there.
    Note - Does not look at the quality of the alignment.
    
    """
    
    sam_file_location = df['file_location']
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    count = 0
    
    for read in samfile:
        
        if read is not None and read.reference_name == transcript_name:
            
            count = count + 1
            
    return count
df['alpha_wt_count'] = df.apply(get_transcript_count, axis=1, args=['Alpha_GEX_64k_HEX'])
df['alpha_dup_count'] = df.apply(get_transcript_count, axis=1, args=['Alpha_GEX_79k_dup_FAM'])
df['beta_count'] = df.apply(get_transcript_count, axis=1, args=['BETA_new_GEX_FAM'])
df.head()
def get_zero_edit_distance_count(df, transcript_name):
    
    """
    For a goven transcript e.g. Alpha_GEX_64k_HEX get how many of the matches have an NM tag of zero
    i.e. the match was exact.
    
    """
    
    sam_file_location = df['file_location']
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    count = 0
    
    for read in samfile:
        
        if read is not None and read.reference_name == transcript_name:
            
            edit_distance = read.get_tag('NM')
            
            if edit_distance == 0:
                
                count = count + 1
            
    return count
df['alpha_wt_zero_edit_count'] = df.apply(get_zero_edit_distance_count, axis=1, args=['Alpha_GEX_64k_HEX'])
df['alpha_dup_zero_edit_count'] = df.apply(get_zero_edit_distance_count, axis=1, args=['Alpha_GEX_79k_dup_FAM'])
df['beta_zero_edit_count'] = df.apply(get_zero_edit_distance_count, axis=1, args=['BETA_new_GEX_FAM'])
df.head()
def get_transcript_read_count_filtered(df, transcript_name, start, end):
    
    """
    Count hits which cover the bit of the reference we are interested in
    
    """
    
    sam_file_location = df['file_location']
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    iter = samfile.fetch(transcript_name, start, end)
    
    count =0
    
    for read in iter:
        
        if read.reference_start <=3 and read.reference_end >= 41:
            
            count = count +1
            
    return count
df['alpha_read_covers_snps_count'] = df.apply(get_transcript_read_count_filtered,
                                              axis=1,
                                              args=['Alpha_GEX_64k_HEX', 0,45])
df['alpha_dup_read_covers_snps_count'] = df.apply(get_transcript_read_count_filtered,
                                              axis=1,
                                              args=['Alpha_GEX_79k_dup_FAM', 0,45])

df['beta_read_covers_snps_count'] = df.apply(get_transcript_read_count_filtered,
                                              axis=1,
                                              args=['BETA_new_GEX_FAM', 0,45])
df.sort_values(by='alpha_dup_read_covers_snps_count', ascending=False).head()
def get_transcript_read_count_filtered_exact(df, transcript_name, start, end):
    
    """
    Count hits which cover the bit of the reference we are interested in and which are exact.
    
    That is do they cross position 0 - 44 of the transcript ( given by transcipt_name) and \
    have an edit distance from the reference of 0. 
    
    Reads which cross this location will cross all four SNPs needed to tell the 3 transcripts apart.
    
    """
    
    sam_file_location = df['file_location']
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    iter = samfile.fetch(transcript_name, start, end)
    
    count =0
    
    for read in iter:
        
        if read.reference_start <=3 and read.reference_end >= 41:
            
            edit_distance = read.get_tag('NM')
            
            if edit_distance == 0:
                
                count = count + 1
                
    return count
df['alpha_read_covers_snps_count_exact'] = df.apply(get_transcript_read_count_filtered_exact,
                                              axis=1,
                                              args=['Alpha_GEX_64k_HEX', 0,45])

df['alpha_dup_read_covers_snps_count_exact'] = df.apply(get_transcript_read_count_filtered_exact,
                                              axis=1,
                                              args=['Alpha_GEX_79k_dup_FAM', 0,45])

df['beta_read_covers_snps_count_exact'] = df.apply(get_transcript_read_count_filtered_exact,
                                              axis=1,
                                              args=['BETA_new_GEX_FAM', 0,45])
# Quick check to see if any such reads exist

df.sort_values(by='alpha_dup_read_covers_snps_count_exact', ascending=False).head()<jupyter_output><empty_output><jupyter_text>Sanity CheckPrint out some of the reads. I then did a manual BLAST and alignment to check the code was working as expected<jupyter_code>def get_transcript_read_count_filtered_exact_print(sam_file_location, transcript_name, start, end):
    
    """
    Same as get_transcript_read_count_filtered_exact() except we just print out the reads instead of \
    counting them.
    
    
    """
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    iter = samfile.fetch(transcript_name, start, end)
    
    for read in iter:
        
        if read.reference_start <=3 and read.reference_end >= 41:
            
            edit_distance = read.get_tag('NM')
            
            if edit_distance == 0:
                
                print (read)
#Pick a random bam file with some read alignments in this area

get_transcript_read_count_filtered_exact_print('../aml_results/aml_aabspliced_SRR5626188.sam.sorted.bam',
                                               'Alpha_GEX_64k_HEX',
                                               0,
                                               45 )
df.describe()
aml_total_hits = df[df['source']=='AML'].sum(axis=0)
aml_total_hits
control_total_hits = df[df['source']=='control'].sum(axis=0)
control_total_hits
df.head()
def get_bases_at_pos(df, transcript, pos):
    """
    For a given position in the reference what are the bases in the reads which align at \
    that position.
    
    NB - only counts bases that are properly aligned (see same flag 163) so the number of bases \
    may be less than the count of reads aligned fully over the area.
    
    """
    
    sam_file_location = df['file_location']
    
    samfile = pysam.AlignmentFile(sam_file_location, "rb")
    
    iter = samfile.pileup(contig=transcript)
    
    for column in iter:
        
        #print (column.reference_name, column.reference_pos)
        
        if column.reference_pos == pos:
            
            base_list = []
              
            for read in column.pileups:

                if read.is_del == 0:

                    base_list.append(read.alignment.query_sequence[read.query_position])

                    #print(read)

                    #print (read.query_position)

                    #print (read.alignment.query_sequence[read.query_position])

                else:

                    base_list.append('_')
                    
            return ''.join(base_list)
df['alpha_wt_pos3_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_64k_HEX', 3])
df['alpha_dp_pos3_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_79k_dup_FAM', 3])
df['beta_pos3_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['BETA_new_GEX_FAM', 3])
df.head(5)
df['alpha_wt_pos25_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_64k_HEX', 25])
df['alpha_dp_pos25_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_79k_dup_FAM', 25])
df['beta_pos25_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['BETA_new_GEX_FAM', 25])
df['alpha_wt_pos40_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_64k_HEX', 40])
df['alpha_dp_pos40_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_79k_dup_FAM', 40])
df['beta_pos40_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['BETA_new_GEX_FAM', 40])
df['alpha_wt_pos41_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_64k_HEX', 41])
df['alpha_dp_pos41_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['Alpha_GEX_79k_dup_FAM', 41])
df['beta_pos41_pileup'] = df.apply(get_bases_at_pos, axis=1, args=['BETA_new_GEX_FAM', 41])
df.head(50)
df[['file_location', 'source','alpha_wt_pos3_pileup',
       'alpha_dp_pos3_pileup', 'alpha_wt_pos25_pileup',
       'alpha_dp_pos25_pileup',  'alpha_wt_pos40_pileup',
       'alpha_dp_pos40_pileup',  'alpha_wt_pos41_pileup',
       'alpha_dp_pos41_pileup']].head(50)
def combine_pileups(df,columns):
    
    final_string = ''
    
    for column in columns:
        
        if df[column] is not None:
        
            final_string = final_string + df[column]
        
    return final_string
df['pos3_pileup'] = df.apply(combine_pileups, axis=1, args=[['alpha_wt_pos3_pileup','alpha_dp_pos3_pileup']])
df['pos25_pileup'] = df.apply(combine_pileups, axis=1, args=[['alpha_wt_pos25_pileup','alpha_dp_pos25_pileup']])
df['pos40_pileup'] = df.apply(combine_pileups, axis=1, args=[['alpha_wt_pos40_pileup','alpha_dp_pos40_pileup']])
df['pos41_pileup'] = df.apply(combine_pileups, axis=1, args=[['alpha_wt_pos41_pileup','alpha_wt_pos41_pileup']])
df[['file_location', 'source','pos3_pileup','pos25_pileup','pos40_pileup', 'pos41_pileup']].tail()
def max_base(df, column):
    
    base_dict = {}
    
    for base in df[column]:
        
        if base not in base_dict:
            
            base_dict[base]=1
            
        else:
            
            base_dict[base] = base_dict[base] +1
            
    total_count = sum(base_dict.values()
            
    return base_dict
df['pos3_call'] = df.apply(max_base, axis=1, args=['pos3_pileup'])
df['pos25_call'] = df.apply(max_base, axis=1, args=['pos25_pileup'])
df['pos40_call'] = df.apply(max_base, axis=1, args=['pos40_pileup'])
df['pos41_call'] = df.apply(max_base, axis=1, args=['pos41_pileup'])
df[['file_location', 'source','pos3_call','pos25_call','pos40_call','pos41_call']].head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10993528074164771, 'fraction_numerical': 0.021164946650341087, 'mean_word_length': 2.7369281045751634, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '44092', 'license': 'MIT', 'n_tokens_mistral': 4052, 'n_tokens_neox': 3792, 'n_words': 722, 'path': 'analysis/.ipynb_checkpoints/AML_analysis-checkpoint.ipynb', 'repo_name': 'NCBI-Hackathons/Tryptase_Duplication', 'size': '11434'}