samp="M1"

#command for only obtaining the context '__CG_'
!grep "[A-Z][A-Z]CG[A-Z]" <methratio_out_M1.txt> methratio_out_M1_CG.txt 
#3x coverage
!awk '{if ($8 >= 3) print $1,"methratio","CpG",$2,$2+1,$5,".","."}' <methratio_out_M1_CG.txt> filt_methratio_out_M1_CG.gff 
!tr ' ' "\t" <filt_methratio_out_M1_CG.gff> filt3_M1.gff

#make igv
!awk '{if ($8 >= 3) print $1,$2-1,$2+1,"CpG",$5}' <methratio_out_M1_CG.txt> filt_methratio_out_M1_CG.igv 
!tr ' ' "\t" <filt_methratio_out_M1_CG.igv> filt3_M1.igv

#make igv
!grep "[A-Z][A-Z]CG[A-Z]" <methratio_out_M3.txt> methratio_out_M3_CG.txt 
!awk '{if ($8 >= 3) print $1,$2-1,$2+1,"CpG",$5}' <methratio_out_M3_CG.txt> filt_methratio_out_M3_CG.igv 
!tr ' ' "\t" <filt_methratio_out_M3_CG.igv> filt3_M3.igv

#make igv with join column
!awk '{if ($8 >= 3) print $1"_"$2-1,$1,$2-1,$2+1,"CpG",$5,$8}' <methratio_out_M3_CG.txt> filt_methratio_out_M3_CGj.txt 
!tr ' ' "\t" <filt_methratio_out_M3_CGj.txt> filt3_M3j.txt

#make igv with join column
!awk '{if ($8 >= 3) print $1"_"$2-1,$1,$2-1,$2+1,"CpG",$5,$8}' <methratio_out_M1_CG.txt> filt_methratio_out_M1_CGj.txt 
!tr ' ' "\t" <filt_methratio_out_M1_CGj.txt> filt3_M1j.txt

!echo "M3ID M3scaffold M3start M3end M3feat M3ratio M3coverage" >> m3head

!cat m3head filt_methratio_out_M3_CGj.txt > hdfilt_methratio_out_M3_CGj.txt

!tr ' ' "\t" <hdfilt_methratio_out_M1_CGj.txt> filt3_M1j.txt

!tr ' ' "\t" <hdfilt_methratio_out_M3_CGj.txt> filt3_M3j.txt

smp="M1"
#Need to manually change id in awk line.

smp="M3"

smp="T1D3"

smp="T1D5"

smp="T3D3"

smp="T3D5"

SELECT * FROM [sr320@washington.edu].[filt3_M1]m1
 join
  [sr320@washington.edu].[filt3_M3]m3
  on
  m1.M1ID=m3.M3ID
  left join
    [sr320@washington.edu].[CG_Ensembl_featurejoin]cg
    on
    m1.M1ID=cg.ChrStart

#confirmed will work in theory? meaning joined to feature

SELECT * FROM [sr320@washington.edu].[filt3M1_M3_CGfeature]
  where M1ratio > 0.5
  and
  M3ratio > 0.5

!python /Users/sr320/sqlshare-pythonclient/tools/fetchdata.py -s "SELECT * FROM [sr320@washington.edu].[filt3M1_M3_CGfeature] where M1ratio > 0.5 and M3ratio > 0.5" -f tsv -o /Volumes/web/cnidarian/M1_M3_05.txt

!head /Volumes/web/cnidarian/M1_M3_05.txt

#above was done with join on feature tracks
#lets do with just simple Joins on samples

#lets do heatmaps
#need to get Methylation levels for all 6 samples (common loci)

!python /Users/sr320/sqlshare-pythonclient/tools/fetchdata.py -s "SELECT M1ID,M1ratio,M3ratio,T1D3ratio,T1D5ratio,T3D3ratio,T3D5ratio FROM [sr320@washington.edu].[filt3_M1]m1 join [sr320@washington.edu].[filt3_M3]m3 on m1.M1ID=m3.M3ID join [sr320@washington.edu].[filt3_T1D3]t1d3 on m1.M1ID=t1d3.T1D3ID join [sr320@washington.edu].[filt3_T1D5]t1d5 on m1.M1ID=t1d5.T1D5ID join [sr320@washington.edu].[filt3_T3D3]t3d3 on m1.M1ID=t3d3.T3D3ID join [sr320@washington.edu].[filt3_T3D5]t3d5 on m1.M1ID=t3d5.T3D5ID" -f tsv -o /Volumes/web/cnidarian/filt3_BiGoLarvae_common.txt

!python /Users/sr320/sqlshare-pythonclient/tools/fetchdata.py -s "SELECT M1ratio,M3ratio,T1D3ratio,T1D5ratio,T3D3ratio,T3D5ratio FROM [sr320@washington.edu].[filt3_M1]m1 join [sr320@washington.edu].[filt3_M3]m3 on m1.M1ID=m3.M3ID join [sr320@washington.edu].[filt3_T1D3]t1d3 on m1.M1ID=t1d3.T1D3ID join [sr320@washington.edu].[filt3_T1D5]t1d5 on m1.M1ID=t1d5.T1D5ID join [sr320@washington.edu].[filt3_T3D3]t3d3 on m1.M1ID=t3d3.T3D3ID join [sr320@washington.edu].[filt3_T3D5]t3d5 on m1.M1ID=t3d5.T3D5ID" -f tsv -o /Volumes/web/cnidarian/filt3_BiGoLarvae_common2.txt

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import scipy.spatial.distance as distance
import scipy.cluster.hierarchy as sch

from pandas import *

# read data from data file into a pandas DataFrame  
larDB = read_table("http://eagle.fish.washington.edu/cnidarian/filt3_BiGoLarvae_common3.txt", # name of the data file
            #sep=",", # what character separates each column?
            na_values=["", " "]) # what values should be considered "blank" values?

# helper for cleaning up axes by removing ticks, tick labels, frame, etc.
def clean_axis(ax):
    """Remove ticks, tick labels, and frame from axis"""
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])
    for sp in ax.spines.values():
        sp.set_visible(False)

# cluster
clusters = sch.linkage(pairwise_dists,method='complete')

!python /Users/sr320/sqlshare-pythonclient/tools/fetchdata.py -s "SELECT M1ID,M1ratio,M3ratio,T1D3ratio,T1D5ratio,T3D3ratio,T3D5ratio FROM [sr320@washington.edu].[filt3_M1]m1 join [sr320@washington.edu].[filt3_M3]m3 on m1.M1ID=m3.M3ID join [sr320@washington.edu].[filt3_T1D3]t1d3 on m1.M1ID=t1d3.T1D3ID join [sr320@washington.edu].[filt3_T1D5]t1d5 on m1.M1ID=t1d5.T1D5ID join [sr320@washington.edu].[filt3_T3D3]t3d3 on m1.M1ID=t3d3.T3D3ID join [sr320@washington.edu].[filt3_T3D5]t3d5 on m1.M1ID=t3d5.T3D5ID where [M1coverage] >= '5' and [M3coverage] >= '5' and [T1D3coverage] >= '5' and [T1D5coverage] >= '5' and [T3D3coverage] >= '5' and [T3D5coverage] >= '5'" -f tsv -o /Volumes/web/cnidarian/filt5_BiGoLarvae_common.txt

!python /Users/sr320/sqlshare-pythonclient/tools/fetchdata.py -s "SELECT M1ID,M1ratio,M3ratio,T1D3ratio,T1D5ratio,T3D3ratio,T3D5ratio FROM [sr320@washington.edu].[filt3_M1]m1 join [sr320@washington.edu].[filt3_M3]m3 on m1.M1ID=m3.M3ID join [sr320@washington.edu].[filt3_T1D3]t1d3 on m1.M1ID=t1d3.T1D3ID join [sr320@washington.edu].[filt3_T1D5]t1d5 on m1.M1ID=t1d5.T1D5ID join [sr320@washington.edu].[filt3_T3D3]t3d3 on m1.M1ID=t3d3.T3D3ID join [sr320@washington.edu].[filt3_T3D5]t3d5 on m1.M1ID=t3d5.T3D5ID where [M1coverage] >= '5' and [M3coverage] >= '5' and [T1D3coverage] >= '5' and [T1D5coverage] >= '5' and [T3D3coverage] >= '5' and [T3D5coverage] >= '5' and M1ratio = 0 and T1D3ratio = 0 and T1D5ratio = 0 and M3ratio = 0 and T3D3ratio = 0 and T3D5ratio = 0" -f tsv -o /Volumes/web/cnidarian/filt5_BiGoLarvae_allzero.txt

#lets consider taking 674 1 lineaage with cv <0.2 and get info. this is on 11,000 that match across 6

mkdir igv

sp="M3"

#command for only obtaining the context '__CG_'
#!grep "[A-Z][A-Z]CG[A-Z]" </Volumes/web/cnidarian/BiGo_larvae_3/methratio_out_{sp}.txt> {sp}_methratio_CG.txt
#5x coverage
!awk '{if ($8 >= 5) print $1,$2-1,$2+1,"CpG",$5}' <M1_methratio_CG.txt> M1_methratio_CG5x.igv 
!tr ' ' "\t" <M1_methratio_CG5x.igv> M1_CG5x.igv 

#command for only obtaining the context '__CG_'
!grep "[A-Z][A-Z]CG[A-Z]" </Volumes/web/cnidarian/BiGo_larvae_3/methratio_out_{sp}.txt> {sp}_methratio_CG.txt

sp="T1D3"

#command for only obtaining the context '__CG_'
!grep "[A-Z][A-Z]CG[A-Z]" </Volumes/web/cnidarian/BiGo_larvae_3/methratio_out_{sp}.txt> {sp}_methratio_CG.txt

sp="T1D5"

sp="M3"

sp="T3D3"

sp="T3D5"