%matplotlib inline
%load_ext autoreload
   
import os
import requests
from IPython.display import display, Image, FileLink
#from jupyter_datatables import init_datatables_mode
#from ipydatagrid import DataGrid
from itables import init_notebook_mode
init_notebook_mode(all_interactive=True)
#init_notebook_mode(connected=True)
import itables.options as itablesopt
# display table options
#itablesopt.lengthMenu = [2, 5, 10, 20, 50, 100, 200, 500]
#itablesopt.maxBytes = 10000 , 0 to disable
#itablesopt.maxRows
#itablesopt.maxColumns

#%load_ext jupyter_require
#%requirejs d3 https://d3js.org/d3.v5.min
#init_datatables_mode()

%autoreload 2
import omics_api
from omics_api import *

REQUERY_URL=False
DEBUGGNG=False
myunit='percent_of_dw'

if DEBUGGNG:
    setVariables(loglevel=getSetting('SHOW_DEBUG'),keep_unit=myunit)
else:
    setVariables(loglevel=getSetting('SHOW_ERRORONLY'),keep_unit=myunit)
#setVariables(loglevel=getSetting('SHOW_ERRORONLY'))

"set variables: \nkeep_unit=percent_of_dw\nr_path=Rscript\nplink_path=plink\nLOG_LEVEL=1\nunit_converter={('percent_of_dw', 'ug_per_gdw'): 10000.0, ('ug_per_gdw', 'percent_of_dw'): 0.0001}"

urlgenes='https://icgrc.info/api/user/gene/cannabid,cannabic,cannabiv,cannabin?annot=cs10&accession=all' 
labelgenes='cannab_live'
df_raw= read_url(urlgenes, labelgenes,requery=REQUERY_URL)

display(df_raw)

#setVariables(loglevel=getSetting('SHOW_DEBUG'))
setVariables(loglevel=getSetting('SHOW_ERRORONLY'))
DATADIR='data'
#PLINKFILE=DATADIR + '/cs10-merge-wgs7ds-21trich-plink'
#PLINKFILE=DATADIR + '/cs10-wgs7ds-genomicsdb-allsnps.isec_cds.snpeff.plink'
#PLINKFILE=DATADIR + '/cs10-wgs7ds-genomicsdb-gatk4170.allsnpsonly-bcftools.filltags.fmis7maf05exchet9hwe9.plink'
FASTAFILE=DATADIR+'/cs10refseq.fna'
#draw_heatmap=['LOC115696986', 'LOC115697762']  
vi_jaccard_scores=dict()
gene2range=dict()
for idx,row in df_raw.iterrows():
    reg=row['contig'] + '_' + str(row['fmin']) + '_' + str(row['fmax'])  
    gene2range[row['name']]=(row['contig'],str(row['fmin']),str(row['fmax']),row['name'])
    print(reg + ' ' + row['name'] + ' ' + row['gb_keyword'])
    urlloc='https://icgrc.info/api/user/variant/' + row['contig'] + ':' + str(row['fmin']) + '-' + str(row['fmax']) + '?snpds=7ds&ref=cs10' 
    df_snp= read_url(urlloc,labelgenes +'_'+reg,requery=REQUERY_URL)
    #display(df_snp)
    df_hasgroup=df_snp.set_index('property',drop=False)
    df_hasgroup=df_hasgroup.loc[:, df_hasgroup.loc['cultivar_group'].isin([TYPE_PROP,'cultivar_group','Hemp','Drug','BLDT','NLDT','Hemp-type','Drug-type','Drug type feral','Type I','Type III'])]
    df_hasgroup=df_hasgroup.reset_index(drop=True)
    #display(df_hasgroup)
    convert_group={'BLDT':'Drug','NLDT':'Drug','Hemp-type':'Hemp','Drug-type':'Drug','Drug type feral':'Drug','Type I':'Drug','Type III':'Hemp','Drug':'Drug','Hemp':'Hemp'}
    heatmap=False #just compute VI's  #row['name'] in draw_heatmap
    #if heatmap:
    
    if False: # use plink
        vi_jaccard_scores[row['name'] + ' ' + reg + ' ' + row['gb_keyword']]=do_phylo(df_hasgroup, 'cluster_' + reg,K=2,convert_group=convert_group,plink_file=PLINKFILE,plink_extract=[gene2range[row['name']]],heatmap=heatmap,fastapath=FASTAFILE)
    else: # use live query
        vi_jaccard_scores[row['name'] + ' ' + reg + ' ' + row['gb_keyword']]=do_phylo(df_hasgroup, 'cluster_' + reg,K=2,convert_group=convert_group,heatmap=heatmap,fastapath=FASTAFILE)
    #break 
    
listreg=[]
for reg in vi_jaccard_scores.keys():
    print("\n" + reg)
    j=vi_jaccard_scores[reg]
    vi=None
    nvi=None
    chistat=None
    chipval=None
    J=None
    if not j is None: 
        if j[0]:
            #print("\nvi=" + str(j[0]))
            #print("\njaccard=" + str(j[1]))
            vi=j[0]
            nvi=j[1]
            J=str(j[2])
            #j=(vi,nvi, mapGroupClust2Jaccard, mapGroup2Sample,  mapCluster2Sample)
            chistat=j[5]
            chipval=j[6]
    listreg.append([reg,vi,nvi,chistat,chipval,J])

import pandas as pd
pdvi1=pd.DataFrame(listreg,columns=['region','VI','Normalized VI','chisq stat','chisq pvalue','Jaccard'])
display(pdvi1)
pdvi1.to_csv(labelgenes + '_clusters_vi_' + str(1) +'.tsv',sep="\t")

"set variables: \nkeep_unit=percent_of_dw\nr_path=Rscript\nplink_path=plink\nLOG_LEVEL=1\nunit_converter={('percent_of_dw', 'ug_per_gdw'): 10000.0, ('ug_per_gdw', 'percent_of_dw'): 0.0001}"

NC_044374.1_52617968_52619389 LOC115717841 cannabidiolic acid synthase-like

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=1.48283343309618  Normalized VI=0.8675224287712751'

chisquare_contingency
xs=['Hemp', 'Drug']
ys=['1', '2']
[[53 86]
 [77  7]]

'chi2 statistic: 59.54571072493299  p-value: 1.1948628366335583e-14  degrees of freedom: 1  expected frequencies: [[81.03139013 57.96860987]\n [48.96860987 35.03139013]]'

'Chisquare statistics=59.54571072493299  pvalue=1.1948628366335583e-14'

NC_044375.1_15287394_15289081 LOC115718836 tetrahydrocannabinolic acid synthase-like

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=1.492752381195814  Normalized VI=0.8909203982324857'

# 'LOC115697762' best
# 'LOC115696986' worst
show_clusters=['LOC115697762','LOC115696986']
convert_group={'BLDT':'Drug','NLDT':'Drug','Hemp-type':'Hemp','Drug-type':'Drug','Drug type feral':'Drug','Type I':'Drug','Type III':'Hemp'}
gene2range=dict()
for idx,row in df_raw.iterrows():
    gene2range[row['name']]=(row['contig'],str(row['fmin']),str(row['fmax']),row['name'])
    if row['name'] in show_clusters:
        reg=row['contig'] + '_' + str(row['fmin']) + '_' + str(row['fmax'])  
        gene2rangei=[(row['contig'],str(row['fmin']),str(row['fmax']),row['name'])]
        print(reg + ' ' + row['name'] + ' ' + row['gb_keyword'])
        urlloc='https://icgrc.info/api/user/variant/' + row['contig'] + ':' + str(row['fmin']) + '-' + str(row['fmax']) + '?snpds=7ds&ref=cs10' 
        df_snp= read_url(urlloc,labelgenes +'_' + reg,requery=False)
        #do_clustering(df_snp,row['name'],'cultivar_group',convert_group=convert_group,K=2,heatmap=True,plink_file=PLINKFILE,plink_extract=gene2rangei)
        do_clustering(df_snp,row['name'],'cultivar_group',convert_group=convert_group,K=2,heatmap=True,fastapath=FASTAFILE)

NC_044378.1_30980884_30982768 LOC115697762 cannabidiolic acid synthase

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=0.6534187172841676  Normalized VI=0.5335581095538917'

chisquare_contingency
xs=['Drug', 'Hemp']
ys=['1', '2']
[[83  1]
 [ 8 36]]

'chi2 statistic: 87.46009203152062  p-value: 8.599759683363778e-21  degrees of freedom: 1  expected frequencies: [[59.71875 24.28125]\n [31.28125 12.71875]]'

'Chisquare statistics=87.46009203152062  pvalue=8.599759683363778e-21'

'bedtools getfasta -fi data/cs10refseq.fna -bed cluster_LOC115697762.btbed > cluster_LOC115697762.bed.ref'

plink --bfile cluster_LOC115697762 --recode tab --out cluster_LOC115697762_ped --allow-extra-chr

'cluster_LOC115697762_ped.ped.tsv  generated'

python draw_tree_matrix.py cluster_LOC115697762_ped.ped.tsv cluster_LOC115697762.newick

import itertools

G=23
genecombination=list(itertools.combinations(df_raw['name'].tolist(),2))
genecombination=genecombination + list(itertools.combinations(df_raw['name'].tolist(),3))
display(pd.DataFrame(genecombination))

vi_jaccard_scores=dict()
gcount=1
for row in genecombination:

    # query for properties
    regions=[]
    for g in row:
        gene1range=gene2range[g]
        regions.append(gene1range[0] + ':' + gene1range[1] + '-' + gene1range[2])
    regions=",".join(regions)
    
    urlloc='https://icgrc.info/api/user/variant/' + regions + '?snpds=7ds&ref=cs10' 
    #reg=gene1range[0] + '_' + gene1range[1] + '_' + gene1range[2]  
    reg="_".join(row) 
    df_snp= read_url(urlloc,labelgenes + '_' + reg,requery=REQUERY_URL)

    reg="_".join(row) 
    plink_extract=[]
    for combs in list(row):
        plink_extract.append(gene2range[combs])
    print(str(gcount) + '/' + str(len(genecombination)) + '  ' + str(plink_extract))
    #vi_jaccard_scores[reg]=do_clustering(df_snp,reg,'cultivar_group',convert_group=convert_group,K=2,heatmap=False,plink_file=PLINKFILE,plink_extract=plink_extract)
    vi_jaccard_scores[reg]=do_clustering(df_snp,reg,'cultivar_group',convert_group=convert_group,K=2,heatmap=False)
    gcount+=1
    #if gcount>2:
    #    break
listreg=[]
for reg in vi_jaccard_scores.keys():
    print("\n" + reg)
    j=vi_jaccard_scores[reg]
    vi=None
    nvi=None
    J=None
    if not j is None: 
        if j[0]:
            vi=j[0]
            nvi=j[1]
            J=str(j[2])
            chistat=j[5]
            chipval=j[6]
    listreg.append([reg,vi,nvi,chistat,chipval,J])            


pdvi23=pd.DataFrame(listreg,columns=['region','VI','Normalized VI','chisq stat','chisq pvalue','Jaccard'])
display(pdvi23)
pdvi23.to_csv(labelgenes + '_clusters_vi_' + str(G) +'.tsv',sep="\t")

#pdvi3=pd.DataFrame(listreg,columns=['region','VI','Normalized VI','chisq stat','chisq pvalue','Jaccard'])
#display(pdvi3)
#pdvi3.to_csv(labelgenes + '_clusters_vi_' + str(G) +'.tsv',sep="\t")

1/220  [('NC_044374.1', '52617968', '52619389', 'LOC115717841'), ('NC_044375.1', '15287394', '15289081', 'LOC115718836')]

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=1.551982422160159  Normalized VI=0.8918556335950214'

chisquare_contingency
xs=['Drug', 'Hemp']
ys=['1', '2']
[[27 57]
 [37  7]]

'chi2 statistic: 29.125541125541126  p-value: 6.783675403421297e-08  degrees of freedom: 1  expected frequencies: [[42. 42.]\n [22. 22.]]'

'Chisquare statistics=29.125541125541126  pvalue=6.783675403421297e-08'

2/220  [('NC_044374.1', '52617968', '52619389', 'LOC115717841'), ('NC_044375.1', '15324213', '15325849', 'LOC115720716')]

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=1.451272523144821  Normalized VI=0.8588787419836541'

chisquare_contingency
xs=['Drug', 'Hemp']
ys=['1', '2']
[[26 58]
 [39  5]]

#draw_heatmap=['LOC115720716_LOC115696884_LOC115697126','LOC115696884_LOC115697126_LOC115697019','LOC115718836_LOC115697762_LOC115696987']
#draw_heatmap=['LOC115696909_LOC115696987_LOC115696986','LOC115718836_LOC115696987_LOC115696986']
#pdvi_all=pd.concat([pdvi1, pdvi2, pdvi3])


pdvi_all=pd.concat([pdvi1, pdvi23])
printdisplay(MESSAGE_INFO,pdvi_all,'123-cannloci-hempdrup_clusters')

123-cannloci-hempdrup_clusters
['region', 'VI', 'Normalized VI', 'chisq stat', 'chisq pvalue', 'Jaccard']

# filtered
#combinations=['LOC115720716_LOC115696909', 'LOC115718836_LOC115720716_LOC115697762', 'LOC115718836_LOC115697762_LOC115697019']
#live
#combinations=['LOC115696987_LOC115696986','LOC115718836_LOC115696884', 'LOC115696884_LOC115697126']

# live2
combinations=['LOC115718836_LOC115696884','LOC115696884_LOC115697126','LOC115718836_LOC115696909_LOC115696884','LOC115696909_LOC115696884_LOC115697126','LOC115696987_LOC115696986','LOC115696909_LOC115696987_LOC115696986']

vi_jaccard_scores=dict()
for reg in combinations:
    genes=reg.split("_")
    print(reg)
    display(df_raw[df_raw['name'].isin(genes)])
    
    #gene1range=gene2range[genes[0]]
    #urlloc='https://icgrc.info/api/user/variant/' + gene1range[0] + ':' + gene1range[1] + '-' + gene1range[2] + '?snpds=7ds&ref=cs10' 
    #reg=gene1range[0] + '_' + gene1range[1] + '_' + gene1range[2]  
    #plink_extract=[]
    #for combs in genes:
    #    plink_extract.append(gene2range[combs])
    #print(plink_extract)
    
    regions=[]
    for g in genes:
        gene1range=gene2range[g]
        regions.append(gene1range[0] + ':' + gene1range[1] + '-' + gene1range[2])
    regions=",".join(regions)
    
    urlloc='https://icgrc.info/api/user/variant/' + regions + '?snpds=7ds&ref=cs10' 
    df_snp= read_url(urlloc,labelgenes + '_' + reg,requery=False)

    #vi_jaccard_scores[reg]=do_clustering(df_snp,reg,'cultivar_group',convert_group=convert_group,K=2,heatmap=True,plink_file=PLINKFILE,plink_extract=plink_extract)
    vi_jaccard_scores[reg]=do_clustering(df_snp,reg,'cultivar_group',convert_group=convert_group,K=2,heatmap=True,fastapath=FASTAFILE)

LOC115718836_LOC115696884

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading required package: metadat
Loading 'meta' package (version 7.0-0).
Type 'help(meta)' for a brief overview.
Readers of 'Meta-Analysis with R (Use R!)' should install
older version of 'meta' package: https://tinyurl.com/dt4y5drs

'Variation of Information=0.11507166270421601  Normalized VI=0.11716232178800956'

intersect==0 Drug 2

'Chisquare statistics=None  pvalue=None'

'bedtools getfasta -fi data/cs10refseq.fna -bed cluster_LOC115718836_LOC115696884.btbed > cluster_LOC115718836_LOC115696884.bed.ref'

plink --bfile cluster_LOC115718836_LOC115696884 --recode tab --out cluster_LOC115718836_LOC115696884_ped --allow-extra-chr

'cluster_LOC115718836_LOC115696884_ped.ped.tsv  generated'

python draw_tree_matrix.py cluster_LOC115718836_LOC115696884_ped.ped.tsv cluster_LOC115718836_LOC115696884.newick

LOC115696884_LOC115697126

ICGRC Omics API Demo (cannabinoid genes SNP clusters and Hemp/Drug label matching, live query all SNPs)¶