<jupyter_start><jupyter_text>Results - smMIPs Sequencing and Data Analysis  Tools<jupyter_code>#!/usr/bin/env python3
import numpy as np
import pandas as pd
import seaborn as sns
sns.set(style='white')
sns.set_context("talk")
import matplotlib.pyplot as plt
import glob
import os
import requests
from pyliftover import LiftOver
lo = LiftOver('hg19', 'hg38')
li = LiftOver('hg38', 'hg19')<jupyter_output><empty_output><jupyter_text>Pull in Input Files<jupyter_code>samples = pd.read_csv('../data/validation_samples/sample_dataframe.txt', sep='\t')
samples_QC = pd.read_csv('../data/validation_samples/sequencing_quality_check.txt', sep='\t')
overlap_with_smmips = pd.read_csv('../output/variant_overlap.tsv', sep='\t')
WEX_variants = pd.read_csv('../data/original_sequencing/VCF_exome.txt', sep='\t')<jupyter_output><empty_output><jupyter_text>Initial Quality Check<jupyter_code>print('Mean: ', samples_QC['total tags weighted'].mean())
mean = samples_QC['total tags weighted'].mean()
print('Standard Deviation: ', samples_QC['total tags weighted'].std())
std = samples_QC['total tags weighted'].std()
ineligible_samples = []
for i,row in samples_QC.iterrows():
    if row['total tags weighted'] < mean - 1*std or row['total unique mips captured'] < 600000:
        print(row['count'] + '_' +  row['Type'])
        if row['Type'] == 'tumor':
            ineligible_samples.append(row['count'] + '_' +  row['Type'])
print('Number of elgible Samples: ', len(samples_QC) - len(ineligible_samples))
# If tumor is ineligible, eliminate variants for subsequent analysis
overlap_with_smmips['Passed QC'] = ''
for i,row in overlap_with_smmips.iterrows():
    if row['sample'] + '_' +  'tumor' in ineligible_samples:
        overlap_with_smmips.loc[[i], 'Passed QC'] = 'no'
    else:
        overlap_with_smmips.loc[[i], 'Passed QC'] = 'yes'
print('Number of Eligible Individuals: ', len(overlap_with_smmips[overlap_with_smmips['Passed QC'] == 'yes'][['sample']].drop_duplicates()))
print('Number of Eligible Variants: ',overlap_with_smmips[overlap_with_smmips['Passed QC'] == 'yes'].groupby('sample').size().sum())
overlap_with_smmips.to_csv('../output/variant_overlap2.tsv', sep='\t')<jupyter_output><empty_output><jupyter_text>Accuracy profile of smMIPs CIViC panel when compared to WEX   Pull in VCF Files<jupyter_code>smmips_variants = pd.DataFrame()
smmips_variants = pd.DataFrame()
for item in glob.glob('../data/smmips_sequencing/*T*.snp.vcf'):
    name = item.split('/')[3].split('_')[0] + '_' + item.split('/')[3].split('_')[1]
    current = pd.read_csv(item, sep='\t', comment='#', header=None).filter(items=[0,1,1,3,4,9])
    current['sample'] = name
    smmips_variants = smmips_variants.append(current)
for item in glob.glob('../data/smmips_sequencing/*T*.indel.vcf'):
    name = item.split('/')[3].split('_')[0] + '_' + item.split('/')[3].split('_')[1]
    current = pd.read_csv(item, sep='\t', comment='#', header=None).filter(items=[0,1,1,3,4,9])
    current['sample'] = name
    current[1] = current[1] + 1
    current[3] = current[3].str[1:].replace('', '-')
    current[4] = current[4].str[1:].replace('', '-')
    smmips_variants = smmips_variants.append(current)
smmips_variants.columns = ['chrom', 'start', 'stop', 'reference', 'variant', 'CIViC Panel VAF', 'sample']
smmips_variants = smmips_variants.reset_index()
CIViC_smMIPs_average = []
for i,row in smmips_variants.iterrows():
    if len(str(row['CIViC Panel VAF']).split(':')) > 5:
        coverage = str(row['CIViC Panel VAF']).split(':')[6]
        CIViC_smMIPs_average.append(int(coverage))
print('Average coverage for smMIPs variants is:', np.mean(CIViC_smMIPs_average))
print('Standard deviation for coverage of smMIPs variants is:', np.std(CIViC_smMIPs_average))
for i,row in smmips_variants.iterrows():
    VAF = float(str(row['CIViC Panel VAF']).split(':')[-1])*100
    smmips_variants.loc[[i], 'CIViC Panel VAF'] = VAF
overlap_with_smmips = overlap_with_smmips.merge(smmips_variants, on=['chrom', 'start','sample'], how='left')
overlap_with_smmips['CIViC Panel VAF'] = overlap_with_smmips['CIViC Panel VAF'].replace(np.nan, 0)
overlap_with_smmips['CIViC Panel VAF'] = overlap_with_smmips['CIViC Panel VAF'].astype('float')
overlap_with_smmips = overlap_with_smmips.drop(['reference_y', 'stop_x', 'stop_y', 'variant_y', 'index'], axis=1)
print('Total eligible variants: ', len(overlap_with_smmips[overlap_with_smmips['Passed QC'] == 'yes']))
print('Total eligible overlap with smMIPs: ', len(overlap_with_smmips[(overlap_with_smmips['CIViC Panel VAF'] > 0) & (overlap_with_smmips['Passed QC'] == 'yes')]))
overlap_with_smmips = overlap_with_smmips.drop(['Unnamed: 0'], axis=1)<jupyter_output><empty_output><jupyter_text>Build Waterfall Dataframe<jupyter_code>samples_waterfall = overlap_with_smmips[overlap_with_smmips['Passed QC'] == 'yes']
# Make dataframe for WaterFall plot (Figure 1)
samples_waterfall = samples_waterfall.filter(items=['sample', 'gene_name', 'amino_acid', 'VAF', 'CIViC Panel VAF'])
samples_waterfall['Validated'] = (samples_waterfall['CIViC Panel VAF'] > 0).astype('int')

samples_waterfall.to_csv('../data/smMIPs_waterfall/waterfall_dataframe.tsv', sep='\t')<jupyter_output><empty_output><jupyter_text>Variant allele frequency correlation between  smMIPs CIViC panel and exome/genome sequencing<jupyter_code>palette_color ={"yes":"#182C61","no":"#D9B08C",True:"#182C61", False:"#D9B08C", 
                'OSCC':'#2F4454', 'SCLC':'#2E151B', 'CRC':'#DA7B93', 'AML':'#376E6F', 'HL':'#1C3334',
               '100ng':'#2F4454', '250ng':'#2E151B', '500ng':'#DA7B93',
               '0X':'#D7E2E2','0-100X':'#87A8A9','100-250X':'#376E6F','250-500X':'#1C3738','500-1,000X':'#0E1C1C', '>1,000X':'#000000'}
print('Variants obtained from samples that failed CIViC smMIPs capture panel:', len(overlap_with_smmips[overlap_with_smmips['Passed QC'] == 'no']))
sample_info = pd.read_csv('../data/validation_samples/sample_dataframe.txt', sep='\t')
correlation = overlap_with_smmips[(overlap_with_smmips['CIViC Panel VAF'] != 0) & (overlap_with_smmips['VAF'] != 0)]
print('Variants in correlation:', len(correlation))
pd.DataFrame.corr(correlation.filter(items=['CIViC Panel VAF', 'VAF']),method='pearson')
overlap_with_smmips_data = overlap_with_smmips.merge(sample_info, right_on='Sample', left_on='sample')
overlap_with_smmips_data['Status'] = overlap_with_smmips_data['CIViC Panel VAF'] > 0
plt.figure(figsize=(5,5))
overlap_with_smmips_data = overlap_with_smmips_data.rename(columns={'Passed QC':'Sequencing Status'})
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data, alpha=0.8, hue='Sequencing Status', palette=palette_color)
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_PassedQC.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
overlap_with_smmips_data_PASSED = overlap_with_smmips_data[overlap_with_smmips_data['Sequencing Status'] == 'yes']
overlap_with_smmips_data_PASSED.drop(56, inplace=True)
overlap_with_smmips_data_PASSED[overlap_with_smmips_data_PASSED['Sample'] == 'CRC_5']
plt.figure(figsize=(5,5))
overlap_with_smmips_data_PASSED = overlap_with_smmips_data_PASSED.rename(columns={'Status':'Validation Status'})

sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, hue='Validation Status', palette=palette_color)
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_Status.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
plt.figure(figsize=(5,5))
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, hue='sample')
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_Sample.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
plt.figure(figsize=(5,5))
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, hue='Matched Normal', palette=palette_color)
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_MatchedNormal.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
plt.figure(figsize=(5,5))
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, hue='Mass', palette=palette_color)
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_Mass.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
def avg_coverage (row):
    if row['Coverage'] == 0 :
        return '0X'
    if row['Coverage'] > 0 and row['Coverage'] < 100:
        return '0-100X'
    if row['Coverage'] >= 100 and row['Coverage'] < 250:
        return '100-250X'
    if row['Coverage'] >= 250 and row['Coverage'] < 500:
        return '250-500X'
    if row['Coverage'] >= 500 and row['Coverage'] < 1000:
        return '500-1,000X'
    if row['Coverage'] >= 1000:
        return '>1,000X'
overlap_with_smmips_data_PASSED['Average Coverage'] = overlap_with_smmips_data_PASSED.apply(lambda row: avg_coverage (row),axis=1)
plt.figure(figsize=(5,5))
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, palette=palette_color, hue='Average Coverage', 
                hue_order=['0X','0-100X','100-250X','250-500X','500-1,000X', '>1,000X'])
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.savefig('../data/Figures/VAF_correlation_Coverage.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
plt.figure(figsize=(5,5))
sns.scatterplot(x="VAF", y="CIViC Panel VAF", data=overlap_with_smmips_data_PASSED, alpha=0.8, hue='Tumor Type', palette=palette_color)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.xlabel('Original VAF')
plt.ylabel('CIViC smMIPs VAF')
plt.savefig('../data/Figures/VAF_correlation_TumorType.pdf', bbox_inches='tight', dpi=400)
# plt.show()
plt.close()
import scipy.stats as ss
n = 100
p = .01
hh = ss.binom(n, p)
hh.pmf(3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1793786669231509, 'fraction_numerical': 0.034433009521977496, 'mean_word_length': 6.895216400911162, 'pattern_counts': {'":': 2, '<': 31, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '152485', 'license': 'MIT', 'n_tokens_mistral': 4412, 'n_tokens_neox': 4078, 'n_words': 641, 'path': 'smMIPs_Manuscript/notebooks/3. smMIP Sequencing and Data Analytics.ipynb', 'repo_name': 'griffithlab/civic-panel', 'size': '10397'}