%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'
df_raw= read_url(urlgenes, labelgenes,requery=REQUERY_URL)

display(df_raw)

#setVariables(loglevel=getSetting('SHOW_DEBUG'))
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'
PLINKFILE=DATADIR + '/cs10-wgs7ds-fmis5maf2exchet5hwe5qual100k' #600k
#draw_heatmap=['LOC115696986', 'LOC115697762']  
vi_jaccard_scores=dict()
gene2range=dict()
projects=['PRJNA310948', 'PRJNA734114', 'PRJNA866500', 'PRJNA575581']
import itertools

drawtrees=['PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581'] #,'PRJNA310948-PRJNA575581','PRJNA734114-PRJNA575581']

projcombination= list(itertools.combinations(projects,4)) + list(itertools.combinations(projects,1)) + list(itertools.combinations(projects,2)) + list(itertools.combinations(projects,3))

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,reg,requery=REQUERY_URL)
    #display(df_snp)
    
    for prjcomb in projcombination:
        prjcomb=list(prjcomb)
        print(prjcomb)
        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.loc[:, df_hasgroup.loc['NCBI BioProject'].isin([TYPE_PROP,'NCBI BioProject'] + prjcomb)]
        df_hasgroup=df_hasgroup.reset_index(drop=True)
        #print('has group')
        #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'}
        heatmap=False #row['name'] in draw_heatmap
        #if heatmap:
        #vi_jaccard_scores['wgs7ds_600k']=do_phylo(df_hasgroup, 'cluster_' + reg,K=2,convert_group=convert_group,plink_file=PLINKFILE,plink_extract=[])
        # group by project
        comblabel="-".join(prjcomb)
       
        drawtree= len(prjcomb)==1 or comblabel in drawtrees
        pval=None
        if len(prjcomb)>1:
            vi_jaccard_scores['wgs7ds_600k_proj_' + comblabel]=do_phylo(df_hasgroup, 'cluster_byprjn_' + comblabel,K=2,cluster_group='NCBI BioProject',plink_file=PLINKFILE,recalc=True,drawtree=drawtree)
            pval=vi_jaccard_scores['wgs7ds_600k_proj_' + comblabel][6]
            if not drawtree and pval is not None and pval >0.05:
                vi_jaccard_scores['wgs7ds_600k_proj_' + comblabel]=do_phylo(df_hasgroup, 'cluster_byprjn_' + comblabel,K=2,cluster_group='NCBI BioProject',plink_file=PLINKFILE,recalc=True,drawtree=True)

        drawtree=drawtree or (pval is not None and pval >0.05)
        vi_jaccard_scores['wgs7ds_600k_cult_' + comblabel]=do_phylo(df_hasgroup, 'cluster_bycult_' + comblabel,K=2,cluster_group='cultivar_group',convert_group=convert_group,plink_file=PLINKFILE,recalc=True,drawtree=drawtree)

            
    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]:
            #print("\nvi=" + str(j[0]))
            #print("\njaccard=" + str(j[1]))
            vi=j[0]
            nvi=j[1]
            J=str(j[2])
            #j=(vi, mapGroupClust2Jaccard, mapGroup2Sample,  mapCluster2Sample)
            chistat=j[5]
            chip=j[6]
    listreg.append([reg,vi,nvi,chistat,chip,J])

import pandas as pd
pdviall=pd.DataFrame(listreg,columns=['region','VI','Normalized VI','chi2 statistics','chi2 pvalue','Jaccard'])
display(pdviall)
pdviall.to_csv('wgs7ds_600k_clusters_vi_all.tsv',sep="\t")

['PRJNA310948', 'PRJNA734114', 'PRJNA866500', 'PRJNA575581']

rm: cannot remove ‘cluster_byprjn_PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581*’: No such file or directory

group PRJNA310948
group PRJNA734114
group PRJNA866500
group PRJNA575581

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading 'meta' package (version 6.5-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
Warning message:
In for (i in seq_along(cenv$extra)) { :
  closing unused connection 3 (cluster_byprjn_PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581.bed)

'Variation of Information=0.9176026205214576  Normalized VI=0.48534659180689876'

'Chisquare statistics=None  pvalue=None'

python draw-tree-3.py cluster_byprjn_PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581.newick

rm: cannot remove ‘cluster_bycult_PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581*’: No such file or directory

group Hemp
group Drug

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading 'meta' package (version 6.5-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
Warning message:
In for (i in seq_along(cenv$extra)) { :
  closing unused connection 3 (cluster_bycult_PRJNA310948-PRJNA734114-PRJNA866500-PRJNA575581.bed)

#pdvi_all_byproj=pd.concat([pdvi1,pdvi2, pdvi3])
#display(pdvi_all_byproj)

show_clusters=['LOC115697762','LOC115696986']

prjcomb=['PRJNA310948','PRJNA575581']
        
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,reg,requery=REQUERY_URL)
        display(df_snp)
        df_snp_prjn=df_snp.set_index('property',drop=False)
        df_snp_prjn=df_snp_prjn.loc[:, df_snp_prjn.loc['NCBI BioProject'].isin([TYPE_PROP,'NCBI BioProject'] + prjcomb)]
        df_snp_prjn=df_snp_prjn.reset_index(drop=True)
        display(df_snp_prjn)
        do_clustering(df_snp_prjn,row['name'],'cultivar_group',convert_group=convert_group,K=2,heatmap=False,plink_file=PLINKFILE)
        break

NC_044378.1_30980884_30982768 LOC115697762 cannabidiolic acid synthase

group Drug
group Hemp

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading 'meta' package (version 6.5-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.0346707348843731  Normalized VI=0.26369637102311905'

'Chisquare statistics=None  pvalue=None'

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'])

import itertools

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

[('LOC115717841',),
 ('LOC115718836',),
 ('LOC115720716',),
 ('LOC115697880',),
 ('LOC115696909',),
 ('LOC115696884',),
 ('LOC115697762',),
 ('LOC115696987',),
 ('LOC115696986',),
 ('LOC115697126',),
 ('LOC115697019',),
 ('LOC115717841', 'LOC115718836'),
 ('LOC115717841', 'LOC115720716'),
 ('LOC115717841', 'LOC115697880'),
 ('LOC115717841', 'LOC115696909'),
 ('LOC115717841', 'LOC115696884'),
 ('LOC115717841', 'LOC115697762'),
 ('LOC115717841', 'LOC115696987'),
 ('LOC115717841', 'LOC115696986'),
 ('LOC115717841', 'LOC115697126'),
 ('LOC115717841', 'LOC115697019'),
 ('LOC115718836', 'LOC115720716'),
 ('LOC115718836', 'LOC115697880'),
 ('LOC115718836', 'LOC115696909'),
 ('LOC115718836', 'LOC115696884'),
 ('LOC115718836', 'LOC115697762'),
 ('LOC115718836', 'LOC115696987'),
 ('LOC115718836', 'LOC115696986'),
 ('LOC115718836', 'LOC115697126'),
 ('LOC115718836', 'LOC115697019'),
 ('LOC115720716', 'LOC115697880'),
 ('LOC115720716', 'LOC115696909'),
 ('LOC115720716', 'LOC115696884'),
 ('LOC115720716', 'LOC115697762'),
 ('LOC115720716', 'LOC115696987'),
 ('LOC115720716', 'LOC115696986'),
 ('LOC115720716', 'LOC115697126'),
 ('LOC115720716', 'LOC115697019'),
 ('LOC115697880', 'LOC115696909'),
 ('LOC115697880', 'LOC115696884'),
 ('LOC115697880', 'LOC115697762'),
 ('LOC115697880', 'LOC115696987'),
 ('LOC115697880', 'LOC115696986'),
 ('LOC115697880', 'LOC115697126'),
 ('LOC115697880', 'LOC115697019'),
 ('LOC115696909', 'LOC115696884'),
 ('LOC115696909', 'LOC115697762'),
 ('LOC115696909', 'LOC115696987'),
 ('LOC115696909', 'LOC115696986'),
 ('LOC115696909', 'LOC115697126'),
 ('LOC115696909', 'LOC115697019'),
 ('LOC115696884', 'LOC115697762'),
 ('LOC115696884', 'LOC115696987'),
 ('LOC115696884', 'LOC115696986'),
 ('LOC115696884', 'LOC115697126'),
 ('LOC115696884', 'LOC115697019'),
 ('LOC115697762', 'LOC115696987'),
 ('LOC115697762', 'LOC115696986'),
 ('LOC115697762', 'LOC115697126'),
 ('LOC115697762', 'LOC115697019'),
 ('LOC115696987', 'LOC115696986'),
 ('LOC115696987', 'LOC115697126'),
 ('LOC115696987', 'LOC115697019'),
 ('LOC115696986', 'LOC115697126'),
 ('LOC115696986', 'LOC115697019'),
 ('LOC115697126', 'LOC115697019'),
 ('LOC115717841', 'LOC115718836', 'LOC115720716'),
 ('LOC115717841', 'LOC115718836', 'LOC115697880'),
 ('LOC115717841', 'LOC115718836', 'LOC115696909'),
 ('LOC115717841', 'LOC115718836', 'LOC115696884'),
 ('LOC115717841', 'LOC115718836', 'LOC115697762'),
 ('LOC115717841', 'LOC115718836', 'LOC115696987'),
 ('LOC115717841', 'LOC115718836', 'LOC115696986'),
 ('LOC115717841', 'LOC115718836', 'LOC115697126'),
 ('LOC115717841', 'LOC115718836', 'LOC115697019'),
 ('LOC115717841', 'LOC115720716', 'LOC115697880'),
 ('LOC115717841', 'LOC115720716', 'LOC115696909'),
 ('LOC115717841', 'LOC115720716', 'LOC115696884'),
 ('LOC115717841', 'LOC115720716', 'LOC115697762'),
 ('LOC115717841', 'LOC115720716', 'LOC115696987'),
 ('LOC115717841', 'LOC115720716', 'LOC115696986'),
 ('LOC115717841', 'LOC115720716', 'LOC115697126'),
 ('LOC115717841', 'LOC115720716', 'LOC115697019'),
 ('LOC115717841', 'LOC115697880', 'LOC115696909'),
 ('LOC115717841', 'LOC115697880', 'LOC115696884'),
 ('LOC115717841', 'LOC115697880', 'LOC115697762'),
 ('LOC115717841', 'LOC115697880', 'LOC115696987'),
 ('LOC115717841', 'LOC115697880', 'LOC115696986'),
 ('LOC115717841', 'LOC115697880', 'LOC115697126'),
 ('LOC115717841', 'LOC115697880', 'LOC115697019'),
 ('LOC115717841', 'LOC115696909', 'LOC115696884'),
 ('LOC115717841', 'LOC115696909', 'LOC115697762'),
 ('LOC115717841', 'LOC115696909', 'LOC115696987'),
 ('LOC115717841', 'LOC115696909', 'LOC115696986'),
 ('LOC115717841', 'LOC115696909', 'LOC115697126'),
 ('LOC115717841', 'LOC115696909', 'LOC115697019'),
 ('LOC115717841', 'LOC115696884', 'LOC115697762'),
 ('LOC115717841', 'LOC115696884', 'LOC115696987'),
 ('LOC115717841', 'LOC115696884', 'LOC115696986'),
 ('LOC115717841', 'LOC115696884', 'LOC115697126'),
 ('LOC115717841', 'LOC115696884', 'LOC115697019'),
 ('LOC115717841', 'LOC115697762', 'LOC115696987'),
 ('LOC115717841', 'LOC115697762', 'LOC115696986'),
 ('LOC115717841', 'LOC115697762', 'LOC115697126'),
 ('LOC115717841', 'LOC115697762', 'LOC115697019'),
 ('LOC115717841', 'LOC115696987', 'LOC115696986'),
 ('LOC115717841', 'LOC115696987', 'LOC115697126'),
 ('LOC115717841', 'LOC115696987', 'LOC115697019'),
 ('LOC115717841', 'LOC115696986', 'LOC115697126'),
 ('LOC115717841', 'LOC115696986', 'LOC115697019'),
 ('LOC115717841', 'LOC115697126', 'LOC115697019'),
 ('LOC115718836', 'LOC115720716', 'LOC115697880'),
 ('LOC115718836', 'LOC115720716', 'LOC115696909'),
 ('LOC115718836', 'LOC115720716', 'LOC115696884'),
 ('LOC115718836', 'LOC115720716', 'LOC115697762'),
 ('LOC115718836', 'LOC115720716', 'LOC115696987'),
 ('LOC115718836', 'LOC115720716', 'LOC115696986'),
 ('LOC115718836', 'LOC115720716', 'LOC115697126'),
 ('LOC115718836', 'LOC115720716', 'LOC115697019'),
 ('LOC115718836', 'LOC115697880', 'LOC115696909'),
 ('LOC115718836', 'LOC115697880', 'LOC115696884'),
 ('LOC115718836', 'LOC115697880', 'LOC115697762'),
 ('LOC115718836', 'LOC115697880', 'LOC115696987'),
 ('LOC115718836', 'LOC115697880', 'LOC115696986'),
 ('LOC115718836', 'LOC115697880', 'LOC115697126'),
 ('LOC115718836', 'LOC115697880', 'LOC115697019'),
 ('LOC115718836', 'LOC115696909', 'LOC115696884'),
 ('LOC115718836', 'LOC115696909', 'LOC115697762'),
 ('LOC115718836', 'LOC115696909', 'LOC115696987'),
 ('LOC115718836', 'LOC115696909', 'LOC115696986'),
 ('LOC115718836', 'LOC115696909', 'LOC115697126'),
 ('LOC115718836', 'LOC115696909', 'LOC115697019'),
 ('LOC115718836', 'LOC115696884', 'LOC115697762'),
 ('LOC115718836', 'LOC115696884', 'LOC115696987'),
 ('LOC115718836', 'LOC115696884', 'LOC115696986'),
 ('LOC115718836', 'LOC115696884', 'LOC115697126'),
 ('LOC115718836', 'LOC115696884', 'LOC115697019'),
 ('LOC115718836', 'LOC115697762', 'LOC115696987'),
 ('LOC115718836', 'LOC115697762', 'LOC115696986'),
 ('LOC115718836', 'LOC115697762', 'LOC115697126'),
 ('LOC115718836', 'LOC115697762', 'LOC115697019'),
 ('LOC115718836', 'LOC115696987', 'LOC115696986'),
 ('LOC115718836', 'LOC115696987', 'LOC115697126'),
 ('LOC115718836', 'LOC115696987', 'LOC115697019'),
 ('LOC115718836', 'LOC115696986', 'LOC115697126'),
 ('LOC115718836', 'LOC115696986', 'LOC115697019'),
 ('LOC115718836', 'LOC115697126', 'LOC115697019'),
 ('LOC115720716', 'LOC115697880', 'LOC115696909'),
 ('LOC115720716', 'LOC115697880', 'LOC115696884'),
 ('LOC115720716', 'LOC115697880', 'LOC115697762'),
 ('LOC115720716', 'LOC115697880', 'LOC115696987'),
 ('LOC115720716', 'LOC115697880', 'LOC115696986'),
 ('LOC115720716', 'LOC115697880', 'LOC115697126'),
 ('LOC115720716', 'LOC115697880', 'LOC115697019'),
 ('LOC115720716', 'LOC115696909', 'LOC115696884'),
 ('LOC115720716', 'LOC115696909', 'LOC115697762'),
 ('LOC115720716', 'LOC115696909', 'LOC115696987'),
 ('LOC115720716', 'LOC115696909', 'LOC115696986'),
 ('LOC115720716', 'LOC115696909', 'LOC115697126'),
 ('LOC115720716', 'LOC115696909', 'LOC115697019'),
 ('LOC115720716', 'LOC115696884', 'LOC115697762'),
 ('LOC115720716', 'LOC115696884', 'LOC115696987'),
 ('LOC115720716', 'LOC115696884', 'LOC115696986'),
 ('LOC115720716', 'LOC115696884', 'LOC115697126'),
 ('LOC115720716', 'LOC115696884', 'LOC115697019'),
 ('LOC115720716', 'LOC115697762', 'LOC115696987'),
 ('LOC115720716', 'LOC115697762', 'LOC115696986'),
 ('LOC115720716', 'LOC115697762', 'LOC115697126'),
 ('LOC115720716', 'LOC115697762', 'LOC115697019'),
 ('LOC115720716', 'LOC115696987', 'LOC115696986'),
 ('LOC115720716', 'LOC115696987', 'LOC115697126'),
 ('LOC115720716', 'LOC115696987', 'LOC115697019'),
 ('LOC115720716', 'LOC115696986', 'LOC115697126'),
 ('LOC115720716', 'LOC115696986', 'LOC115697019'),
 ('LOC115720716', 'LOC115697126', 'LOC115697019'),
 ('LOC115697880', 'LOC115696909', 'LOC115696884'),
 ('LOC115697880', 'LOC115696909', 'LOC115697762'),
 ('LOC115697880', 'LOC115696909', 'LOC115696987'),
 ('LOC115697880', 'LOC115696909', 'LOC115696986'),
 ('LOC115697880', 'LOC115696909', 'LOC115697126'),
 ('LOC115697880', 'LOC115696909', 'LOC115697019'),
 ('LOC115697880', 'LOC115696884', 'LOC115697762'),
 ('LOC115697880', 'LOC115696884', 'LOC115696987'),
 ('LOC115697880', 'LOC115696884', 'LOC115696986'),
 ('LOC115697880', 'LOC115696884', 'LOC115697126'),
 ('LOC115697880', 'LOC115696884', 'LOC115697019'),
 ('LOC115697880', 'LOC115697762', 'LOC115696987'),
 ('LOC115697880', 'LOC115697762', 'LOC115696986'),
 ('LOC115697880', 'LOC115697762', 'LOC115697126'),
 ('LOC115697880', 'LOC115697762', 'LOC115697019'),
 ('LOC115697880', 'LOC115696987', 'LOC115696986'),
 ('LOC115697880', 'LOC115696987', 'LOC115697126'),
 ('LOC115697880', 'LOC115696987', 'LOC115697019'),
 ('LOC115697880', 'LOC115696986', 'LOC115697126'),
 ('LOC115697880', 'LOC115696986', 'LOC115697019'),
 ('LOC115697880', 'LOC115697126', 'LOC115697019'),
 ('LOC115696909', 'LOC115696884', 'LOC115697762'),
 ('LOC115696909', 'LOC115696884', 'LOC115696987'),
 ('LOC115696909', 'LOC115696884', 'LOC115696986'),
 ('LOC115696909', 'LOC115696884', 'LOC115697126'),
 ('LOC115696909', 'LOC115696884', 'LOC115697019'),
 ('LOC115696909', 'LOC115697762', 'LOC115696987'),
 ('LOC115696909', 'LOC115697762', 'LOC115696986'),
 ('LOC115696909', 'LOC115697762', 'LOC115697126'),
 ('LOC115696909', 'LOC115697762', 'LOC115697019'),
 ('LOC115696909', 'LOC115696987', 'LOC115696986'),
 ('LOC115696909', 'LOC115696987', 'LOC115697126'),
 ('LOC115696909', 'LOC115696987', 'LOC115697019'),
 ('LOC115696909', 'LOC115696986', 'LOC115697126'),
 ('LOC115696909', 'LOC115696986', 'LOC115697019'),
 ('LOC115696909', 'LOC115697126', 'LOC115697019'),
 ('LOC115696884', 'LOC115697762', 'LOC115696987'),
 ('LOC115696884', 'LOC115697762', 'LOC115696986'),
 ('LOC115696884', 'LOC115697762', 'LOC115697126'),
 ('LOC115696884', 'LOC115697762', 'LOC115697019'),
 ('LOC115696884', 'LOC115696987', 'LOC115696986'),
 ('LOC115696884', 'LOC115696987', 'LOC115697126'),
 ('LOC115696884', 'LOC115696987', 'LOC115697019'),
 ('LOC115696884', 'LOC115696986', 'LOC115697126'),
 ('LOC115696884', 'LOC115696986', 'LOC115697019'),
 ('LOC115696884', 'LOC115697126', 'LOC115697019'),
 ('LOC115697762', 'LOC115696987', 'LOC115696986'),
 ('LOC115697762', 'LOC115696987', 'LOC115697126'),
 ('LOC115697762', 'LOC115696987', 'LOC115697019'),
 ('LOC115697762', 'LOC115696986', 'LOC115697126'),
 ('LOC115697762', 'LOC115696986', 'LOC115697019'),
 ('LOC115697762', 'LOC115697126', 'LOC115697019'),
 ('LOC115696987', 'LOC115696986', 'LOC115697126'),
 ('LOC115696987', 'LOC115696986', 'LOC115697019'),
 ('LOC115696987', 'LOC115697126', 'LOC115697019'),
 ('LOC115696986', 'LOC115697126', 'LOC115697019')]

vi_jaccard_scores=dict()
#draw_heatmap=['LOC115720716_LOC115696884_LOC115697126','LOC115696884_LOC115697126_LOC115697019','LOC115718836_LOC115697762_LOC115696987']
#draw_heatmap=['LOC115696909_LOC115696987_LOC115696986','LOC115718836_LOC115696987_LOC115696986']
for row in genecombination:

    # query for properties, use one gene
    gene1=row[0]
    gene1range=gene2range[gene1]
    urlloc='https://icgrc.info/api/user/variant/' + gene1range[0] + ':' + gene1range[1] + '-' + gene1range[2] + '?snpds=7ds&ref=cs10' 
    reg=gene1range[0] + '_' + gene1range[1] + '_' + gene1range[2]  
    df_snp= read_url(urlloc,reg,requery=REQUERY_URL)

    df_snp_prjn=df_snp.set_index('property',drop=False)
    df_snp_prjn=df_snp_prjn.loc[:, df_snp_prjn.loc['NCBI BioProject'].isin([TYPE_PROP,'NCBI BioProject'] + prjcomb)]
    df_snp_prjn=df_snp_prjn.reset_index(drop=True)
    #display(df_snp_prjn)

    reg="_".join(row) 
    plink_extract=[]
    for combs in list(row):
        plink_extract.append(gene2range[combs])
    print(plink_extract)
    vi_jaccard_scores[reg]=do_clustering(df_snp_prjn,reg,'cultivar_group',convert_group=convert_group,K=2,heatmap=False,plink_file=PLINKFILE,plink_extract=plink_extract)
    
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]
            chip=j[6]
    listreg.append([reg,vi,nvi,chistat,chip,J])

pdvi2=pd.DataFrame(listreg,columns=['region','VI','Normalized VI','chi2 statistics','chi2 pvalue','Jaccard'])
display(pdvi2)
pdvi2.to_csv(labelgenes + '_clusters_vi_prjn_' + str(G) +'.tsv',sep="\t")

[('NC_044374.1', '52617968', '52619389', 'LOC115717841')]
group Drug
group Hemp

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading 'meta' package (version 6.5-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=2.290316631308452  Normalized VI=0.9697874454976867'

chisquare_contingency
xs=['Drug', 'Hemp']
ys=['1', '2']
[[18 40]
 [ 6  3]]

'chi2 statistic:     2.8924'

'p-value:            0.088996'

'degrees of freedom: 1'

'expected frequencies:'

array([[20.7761194, 37.2238806],
       [ 3.2238806,  5.7761194]])

'Chisquare statistics=2.89243954008138  pvalue=0.0889960735982675'

pdvi_all=pd.concat([pdvi1, pdvi2, pdvi3])
display(pdvi_all)

#combinations=['LOC115720716_LOC115696909', 'LOC115718836_LOC115720716_LOC115697762', 'LOC115718836_LOC115697762_LOC115697019']

# 2 projects
combinations=['LOC115720716_LOC115696884_LOC115696986','LOC115696987_LOC115696986_LOC115697019','LOC115718836_LOC115696884_LOC115697019','LOC115718836_LOC115696884'] 

vi_jaccard_scores=dict()
for reg in combinations:
    genes=reg.split("_")
    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]  
    df_snp= read_url(urlloc,reg,requery=REQUERY_URL)
    df_snp_prjn=df_snp.set_index('property',drop=False)
    df_snp_prjn=df_snp_prjn.loc[:, df_snp_prjn.loc['NCBI BioProject'].isin([TYPE_PROP,'NCBI BioProject'] + prjcomb)]
    df_snp_prjn=df_snp_prjn.reset_index(drop=True)
    
    plink_extract=[]
    for combs in genes:
        plink_extract.append(gene2range[combs])
    print(plink_extract)
    vi_jaccard_scores[reg]=do_clustering(df_snp_prjn,reg+'_2prjn','cultivar_group',convert_group=convert_group,K=2,heatmap=True,plink_file=PLINKFILE,plink_extract=plink_extract,recalc=True)

[('NC_044375.1', '15324213', '15325849', 'LOC115720716'), ('NC_044378.1', '29635083', '29636878', 'LOC115696884'), ('NC_044378.1', '34056091', '34061726', 'LOC115696986')]

rm: cannot remove ‘cluster_LOC115720716_LOC115696884_LOC115696986_2prjn*’: No such file or directory

group Drug
group Hemp

Loading required package: MASS
Loading required package: abind
Loading required package: epitools
Loading required package: meta
Loading 'meta' package (version 6.5-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.4914487084384396  Normalized VI=0.9897832859299888'

chisquare_contingency
xs=['Drug', 'Hemp']
ys=['1', '2']
[[38 20]
 [ 4  5]]

'chi2 statistic:     0.71538'

'p-value:            0.39766'

'degrees of freedom: 1'

'expected frequencies:'

ICGRC Omics API Demo (cannabinoid genes SNP clusters and Hemp/Drug label matching, datasource independent)¶