# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: 'Python 3.7.7 64-bit (''plot'': conda)'
#     language: python
#     name: python37764bitplotconda2c267f18819b461ab1125c6b70f3e730
# ---

# !conda env list

import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib as mpl
# %matplotlib inline
# %config InlineBackend.figure_format="retina"

plt.style.use("fivethirtyeight")
mpl.rcParams['font.family'] = 'Avenir'
plt.rcParams['font.size'] = 18
plt.rcParams['axes.linewidth'] = 2

sns.set_style("whitegrid", {'axes.grid' : False})
sns.set_style("ticks", {"xtick.major.size": 8, "ytick.major.size": 8})

sns.__version__

import os
os.chdir('/Users/pengl7/Downloads/WGS/compare-variants/NIST/')

# %ls -lth

# the long_cleared_QD has removed the reads with DP=0 compared to long_cleared.csv
# it also has correct calculation of QD (QUAL/DP)
# it has 339850 WITH 10 couts less than long_cleared.csv
df = pd.read_csv("long_cleared_QD.csv")

print(df.shape)
print()
df.head()

df.dtypes

# for the sake of plotting, keep the TYPE as float instead of categorical
#df["TYPE"] = df["TYPE"].astype("category")
df["UNMATCH"] = df["UNMATCH"].astype("category")

df.describe()

df.to_csv("long_cleared_QD.csv", index=None)

# +
# some choices of plot paramters

# using mulitple="stack"
# ax=sns.displot(df, x="QUAL", hue="UNMATCH", kind="kde", multiple="stack", log_scale=True, legend=False)
#  missleading in the edge

# using fill=True
# ax=sns.displot(df, x="QUAL", hue="UNMATCH", kind="kde", fill=True, log_scale=True, legend=False)
# color is too light
# -

ax=sns.displot(df, x="QUAL", hue="UNMATCH", kind="kde", multiple="stack", log_scale=True, legend=False)
plt.legend(labels=['UnMatch', 'Match'])
ax.set(xlabel="QUAL (log)")

ax=sns.displot(df, x="QUAL", hue="UNMATCH", kind="kde", log_scale=True, legend=False)
plt.legend(labels=['UnMatch', 'Match'])
ax.set(xlabel="QUAL (log)")

ax=sns.displot(df, x="QUAL", hue="UNMATCH", kind="kde", fill=True, log_scale=True, legend=False)
plt.legend(labels=['UnMatch', 'Match'])
ax.set(xlabel="QUAL (log)")

sns.displot(df, x="AF", hue="UNMATCH", kind="kde", fill=True, legend=False)
plt.legend(loc='upper left', labels=['UnMatch', 'Match'])

sns.displot(df, x="AF", hue="UNMATCH", legend=False, kind="kde")
plt.legend(loc='upper left', labels=['UnMatch', 'Match'])

# +
# get the variants with AF < 0.05

filt_af = df["AF"]< 0.05
maf_df = df[filt_af]
print(maf_df.shape)

# +
# calculate the number of unmatch in these 897 

maf_df["UNMATCH"].value_counts()
# -

sns.displot(maf_df, x="AF", hue="UNMATCH", kind="kde", legend=False)
plt.legend(labels=['UnMatch', 'Match'])

844/(844+53)

sns.displot(df, x="GQ", hue="UNMATCH", kind="kde", fill=True, legend=False)
plt.legend(loc='upper left', labels=['UnMatch', 'Match'])

ax=sns.displot(df, x="FS", hue="UNMATCH", kind="kde", fill=True, legend=False)
ax.set(xlim=(0,20))
plt.legend(labels=['UnMatch', 'Match'])

np.log(60)

ax=sns.displot(df, x="MQ", hue="UNMATCH", kind="kde", fill=True,legend=False)
ax.set(xlim=(55,65))
plt.legend(loc="upper left", labels=['UnMatch', 'Match'])

df["MQ"].max()

df["MQ"].min()

ax = sns.displot(df, x="DP", hue="UNMATCH", kind="kde", fill=True, legend=False)
ax.set(xlim=(0, 100))
plt.legend(labels=['UnMatch', 'Match'])

print(len(df[df["DP"]==0]))
print(len(df[df["DP"]<=10]))

ax = sns.displot(df, x="QD", hue="UNMATCH", kind="kde", fill=True, legend=False)
ax.set(xlim=(0, 100))
plt.legend(labels=['UnMatch', 'Match'])

ax=sns.displot(df, x="TYPE", hue="UNMATCH", multiple="stack", legend=False)
plt.legend(labels=['UnMatch', 'Match'])
ax.set_xticklabels(["", 'SNP', "", 'INDEL'])

df["TYPE"] = df["TYPE"].astype("int64")

ax=sns.displot(df, x="TYPE", hue="UNMATCH", multiple="stack", legend=False, discrete=True, shrink=0.3)
plt.legend(labels=['UnMatch', 'Match'])
ax.set_xticklabels(['', 'SNP', '',  'INDEL'])

# ## Check the quality of the whole vcf files

# +
# draw venn diagram
# -

import matplotlib_venn as venn

LONG VS REF ON chr1
Found 305039 sites common to both files.
Found 34553 sites only in main file.
Found 9939 sites only in second file.
Found 290 non-matching overlapping sites.
After filtering, kept 339882 out of a possible 339882 Sites


def cal_overlap(c, x_uniq, y_uniq):
    perX = 100*c / (c+x_uniq)
    print("1st overlap: ", perX)
    perY = 100*c / (c+y_uniq)
    print("2nd overlap: ", perY)


cal_overlap(305039,34553,9939)

# +
color1 = "darkorange"
color2 = "teal"

plt.title('comparison of long read NIST with its strandard reference (chr1)')
v1=venn.venn2(subsets = (34843, 9939, 305039), set_labels = ('Long read NIST', 'Standard NIST'), 
         set_colors=(color1, color2))
# v1.get_label_by_id('10').set_text('89.8%')
# v1.get_label_by_id('01').set_text('96.9%')
# -

# ## Check the whole call sets

REF VS LONG
Found 3901271 sites common to both files.
Found 129076 sites only in main file.
Found 593117 sites only in second file.
Found 8854 non-matching overlapping sites.
After filtering, kept 4039201 out of a possible 4039201 Sites

cal_overlap(3901271,593117,129076)

plt.title('comparison of long read NIST with its strandard reference')
v1=venn.venn2(subsets = (593117,129076,3901271), set_labels = ('Long read NIST', 'Standard NIST'), 
         set_colors=(color1, color2))
# v1.get_label_by_id('10').set_text('89.8%')
# v1.get_label_by_id('01').set_text('96.9%')

REF vs SHORT
Found 3850385 sites common to both files.
Found 108615 sites only in main file.
Found 832304 sites only in second file.
Found 80201 non-matching overlapping sites.
After filtering, kept 4039201 out of a possible 4039201 Sites

color3="blue"
cal_overlap(3850385,832304,108615)

plt.title('comparison of short read NIST with its strandard reference')
v2=venn.venn2(subsets = (832304,108615,3850385), set_labels = ('Short read NIST', 'Standard NIST'), 
         set_colors=(color3, color2))

# compare short and long
LONG VS SHORT
Found 4154024 sites common to both files.
Found 235884 sites only in main file.
Found 495532 sites only in second file.
Found 113334 non-matching overlapping sites.
After filtering, kept 4503242 out of a possible 4503242 Sites

#  cal_overlap(common, left_uniq, right_uniq)
cal_overlap(4154024, 235884, 495532)

plt.title('comparison of long read NIST with short read NIST')
v4=venn.venn2(subsets = (235884, 495532,4154024), set_labels = ('Long read NIST', 'Short read NIST'), 
         set_colors=(color1, color3))


# +
# calcualte the three sets
# -

# # calculate precison, recall of the sequencing platform

# ## Calculate metrics:
#
# Ref vs long: 
#
# REF VS LONG
# Found 3901271 sites common to both files.
# Found 129076 sites only in main file.
# Found 593117 sites only in second file.
# Found 8854 non-matching overlapping sites.
# After filtering, kept 4039201 out of a possible 4039201 Sites
#
# TP = the overlapping: i.e, common
# FN = REF_uniqe
# FP = Long_unique
# TN = 0
#
#
# recall = 100* TP/(TP+FN) 
# precision = 100*TP/(TP+FP)
# accuracy = 100* (TP+FN)/(total)
#
#

# +

def cal_metrics(overlapping, ref_uniq, sample_uniq, TN=0):
    recall = 100* overlapping/(overlapping+ref_uniq) 
    precision = 100*overlapping/(overlapping+sample_uniq)
    accuracy = 100* (overlapping+ref_uniq)/(overlapping+sample_uniq+ref_uniq+TN)
    print(f'Recall is {recall:.2f}%')
    print(f'Precision is {precision:.2f}%')
    print(f'Accuracy is {accuracy:.2f}%')


# -

cal_metrics(3901271, 129076, 593117)



# ## Do filtering on FS and QD

# +
# It is better fo remove FS > 60 for SNPs and FS > 200 for Indel
# Also remove those with Guality by Depth 2
# OR remove those with DP < 10 
# or AF < 0.05
# or INDELS
# -

# the FS scores in our data is 
# the value fo FS don't use ""
filt1= (df["FS"] > 60)&(df["TYPE"]==0)
filt2 = (df["FS"] > 200)&(df["TYPE"]==1)
filt3 = (df["QD"] < 2)
filt4 = (df["DP"] < 10)
filt5 = (df["AF"] < 0.05)
filt6 = (df["AF"] < 0.01)
filt7 = (df["MQ"] < 40)

print("SNPs with FS > 60: ", len(df[filt1]))
print("INDEL with FS > 200: ", len(df[filt2]))
print("variants with QD < 2: ", len(df[filt3]))
print("variants with DP < 10: ", len(df[filt4]))
print("variants with AF < 0.05: ", len(df[filt5]))
print("variants with AF < 0.01: ", len(df[filt6]))
print("variants with MQ < 40: ", len(df[filt7]))

df.shape

df_filt_FS = df[(~filt1) & (~filt2)]
print(df_filt_FS.shape)

df_filt_FS_QD = df_filt_FS[df_filt_FS["QD"] > 2]
print(df_filt_QD.shape)

df_filt_FS_QD_AF = df_filt_FS_QD[df_filt_FS_QD["AF"] >= 0.01]
print(df_filt_FS_QD_AF.shape)

df_filt_FS_QD_AF_MQ = df_filt_FS_QD_AF[df_filt_FS_QD_AF["MQ"] >= 40]
print(df_filt_FS_QD_AF_MQ.shape)

df_filt_FS_QD_AF_MQ_DP = df_filt_FS_QD_AF_MQ[df_filt_FS_QD_AF_MQ["DP"] >= 10]
print(df_filt_FS_QD_AF_MQ_DP.shape)

df_filt_FS_QD_AF_DP_MQ = df_filt_FS_QD_AF_DP[df_filt_FS_QD_AF_DP["MQ"] >= 40]
print(df_filt_FS_QD_AF_DP_MQ.shape)

# remove INDEL
df_filt_FS_QD_AF_DP_MQ_SNP = df_filt_FS_QD_AF_DP_MQ[df_filt_FS_QD_AF_DP_MQ["TYPE"]!=1]
print(df_filt_FS_QD_AF_DP_MQ_SNP.shape)


# ## Calculate the accuracy after each step of filtering
#
# print(f'pi = {x:.2f}')

def cal_cal_2(df, column="UNMATCH"):
    match = df[column].value_counts().iloc[0]
    unmatch = df[column].value_counts().iloc[1]
    accuracy = 100*match/(match + unmatch)
    print(f'{accuracy:.2f}%')  


cal_cal_2(df)

for item in [df, df_filt_FS, df_filt_FS_QD, df_filt_FS_QD_AF, df_filt_FS_QD_AF_MQ, df_filt_FS_QD_AF_MQ_DP, df_filt_FS_QD_AF_DP_MQ_SNP]:
    cal_cal_2(item)

# %ls -lth

after_filt_df = pd.read_csv("sum of variants filtering.csv")

StarcoderdataJupyterScriptsDedupFiltered