<jupyter_start><jupyter_text>Examples<jupyter_code>def split_hannah(input):
    return [line.split('\t')[2] for line in input.split('\n')]
wuhan6 = """CMS001_042_Ra_S23	bin10	NODE_10_length_2444_cov_49.183354	1702	36384680	46.77792961 seg1=PB2
CMS001_042_Ra_S23	bin11	NODE_11_length_2437_cov_68.228814	2361	36384680	64.88994819 seg2=PB1
CMS001_042_Ra_S23	bin12	NODE_12_length_2238_cov_65.551134	2072	36384680	56.94704475 seg3=PA
CMS001_042_Ra_S23	bin13	NODE_15_length_1855_cov_99.279528	2584	36384680	71.01890136 seg4=NP
CMS001_042_Ra_S23	bin3	NODE_28_length_1501_cov_37.367978	787	36384680	21.62998273 seg5=gp
CMS001_042_Ra_S23	bin4	NODE_131_length_861_cov_90.163265	1027	36384680	28.22616552 seg6=hypoth"""

wuhan6_ish_contigs = split_hannah(wuhan6)


whidbey_ish = """CMS002_018b_Rb_S129_L004	bin1	NODE_4_length_2203_cov_29.747883	906	89837146	10.08491521 PA
CMS002_018b_Rb_S129_L004	bin15	NODE_10_length_884_cov_107.267658	1237	89837146	13.76935995 NP
CMS002_018b_Rb_S129_L004	bin15	NODE_9_length_939_cov_28.932715	361	89837146	4.018382329 NP
CMS002_018b_Rb_S129_L004	bin2	NODE_2_length_2430_cov_26.274968	876	89837146	9.750977619 PB1"""

whidbey_ish_contigs = split_hannah(whidbey_ish)


wm3_ish = """CMS001_027_Ra_S16	bin14	NODE_6_length_2203_cov_26.203669	930	29422908	31.6080246 PA
CMS001_027_Ra_S16	bin16	NODE_16_length_1779_cov_20.285546	605	29422908	20.56220955 NP
CMS001_027_Ra_S16	bin17	NODE_3_length_2447_cov_21.006751	882	29422908	29.97664269 PB1
CMS001_027_Ra_S16	bin7	NODE_4_length_2420_cov_26.034571	1026	29422908	34.87078843 PB2"""

wm3_ish_contigs = split_hannah(wm3_ish)


def sort_by(x_list, idx_list):
     return [x for (idx,x) in sorted(zip(idx_list,x_list), key=lambda pair: pair[0])]
def fish_dark_matter(bait_contigs, overlap_cutoff=0.9, aggregation=None):
    bait_clusters = [contig_to_cluster[c] for c in bait_contigs]
    neighborhood = nbhd(bait_clusters, 3)
    df = df_from_subset(neighborhood)
    samples_with_bait = df.loc[:, bait_clusters].sum(axis = 1) == len(bait_clusters)
    n_samples_with_bait = df.loc[:, bait_clusters].sum().mean()
    clusters_containing_bait = df.loc[samples_with_bait].sum() >= (overlap_cutoff * n_samples_with_bait)
    clusters_not_overflowing_bait = (df.sum() * overlap_cutoff) <= n_samples_with_bait

    overlapping_clusters = df.columns[clusters_containing_bait & clusters_not_overflowing_bait]
    
    if aggregation:
        df = df_from_subset(neighborhood, aggregation=aggregation)
    df = df[overlapping_clusters]
    df = df.loc[df.sum(axis = 1) > 0,:]
    
    return df

base_path='/Users/evogytis/Dropbox/Jupyter_notebooks/Biohub/California_mosquitoes/s3_bucket/contigs/'

import os,glob
from Bio import SeqIO
def printContigs(marked):
    base_path='/Users/evogytis/Dropbox/Jupyter_notebooks/Biohub/California_mosquitoes/s3_bucket/contigs/'
    output=[]
    for fpath in glob.glob(os.path.join(base_path,'*','*.fasta')):
        sname=os.path.basename(fpath)
        
        contigs=SeqIO.parse(open(fpath,'r'),'fasta')
        for seq in contigs:
#             print(marked,seq.id)
            if seq.id in marked:
                output.append('>%s|%s\n%s'%(os.path.basename(os.path.dirname(fpath)),seq.id,seq.seq))
    return '\n'.join(output)+'\n'
#     print('>%s|%s\n%s'%(os.path.basename(os.path.dirname(fpath)),seq.id,seq.seq))

def fast_fish(bait_contigs, overlap_cutoff=0.9, aggregation=None):
    bait_clusters = [contig_to_cluster[c] for c in bait_contigs]
    neighborhood = nbhd(bait_clusters, 3)
    df = df_from_subset(neighborhood)
    samples_with_bait = df.loc[:, bait_clusters].sum(axis = 1) == len(bait_clusters)
    n_samples_with_bait = df.loc[:, bait_clusters].sum().mean()
    clusters_containing_bait = df.loc[samples_with_bait].sum() >= (overlap_cutoff * n_samples_with_bait)
    clusters_not_overflowing_bait = (df.sum() * overlap_cutoff) <= n_samples_with_bait

    overlapping_clusters = df.columns[clusters_containing_bait & clusters_not_overflowing_bait]
    
    if aggregation:
        df = df_from_subset(neighborhood, aggregation=aggregation)
    df = df[overlapping_clusters]
    df = df.loc[df.sum(axis = 1) > 0,:]
    
    return [cid for cid in df]
            
def display_dark_matter(df,report_seq=False,fltr=None,toFile=None):
    overlapping_clusters = df.columns
    
    if fltr==None:
        fltr=lambda x: True
    
    print('Found ', len(overlapping_clusters), ' candidate segments.')
    lengths = [get_representative(clusters[cluster_id]).length for cluster_id in overlapping_clusters]
    print('Segment lengths: ', sorted(lengths, reverse=True))
    
    if report_seq==True:
        for cluster_id in overlapping_clusters:
#             if toFile:
#                 toFile.write('%s\n'%('>%s\n%s'%(contig.contig,contig.)))
            out=printContigs([contig.contig for contig in clusters[cluster_id] if fltr(contig)])
            if toFile:
                toFile.write(out)
            else:
                print(out)
    
    for cluster in sort_by(overlapping_clusters, lengths):
        print('---')
        print_cluster(cluster)
#     if (df > 0).mean().mean() < 1:
#         sns.clustermap(df)

def segment_correlations(abundances):
    df = pd.DataFrame(np.corrcoef(abundances.T), columns=abundances.columns, index=abundances.columns)
    sns.clustermap(df)
# From one wuhan 6 contig
display_dark_matter(fish_dark_matter(wuhan6_ish_contigs[:1], 0.8))
# wm3 candidates
display_dark_matter(fish_dark_matter(wm3_ish_contigs))
awuhan6_seed = 'NODE_4_length_2235_cov_10.352641'

rep_seq=False
fltr=None
toFile=None
# toFile=open('/Users/evogytis/Downloads/your_fasta.fasta','w')

display_dark_matter(fish_dark_matter([awuhan6_seed]),report_seq=rep_seq,fltr=fltr,toFile=toFile)
toFile.close()<jupyter_output>Found  9  candidate segments.
Segment lengths:  [6634, 2449, 2442, 2301, 2245, 2123, 1842, 1534, 869]
---
Cluster:  23854
Rep contig:  CMS002_045c_Rb_S185_L004 NODE_14_length_869_cov_33.226010
Rep Length:  869
N samples:  15
Contig Lengths:  857,857,857,853,860,664,857,534,857,863,869,856,856,861,857
---
Cluster:  4999
Rep contig:  CMS002_045f_Rb_S189_L004 NODE_10_length_1534_cov_13.974605
Rep Length:  1534
N samples:  15
Contig Lengths:  1494,1505,1514,1507,1496,1001,1476,899,1505,1510,1510,1515,1511,1534,1513
---
Cluster:  2828
Rep contig:  CMS002_045e_Rb_S188_L004 NODE_7_length_1842_cov_33.126346
Rep Length:  1842
N samples:  15
Contig Lengths:  1838,1823,1815,1839,1830,1079+689,1819,1804,1817,1832,1835,1838,1842,1841,1827
---
Cluster:  1813
Rep contig:  CMS002_026e_Rb_S150_L004 NODE_5_length_2123_cov_21.675953
Rep Length:  2123
N samples:  16
Contig Lengths:  2091,2085,2107,2093,2123,1404,793,2056,2040,2087,2084,2093,2093,2089,2079,2096
---
Cluster:  1513
Rep contig:  CMS002_026c_R[...]<jupyter_text>If we contract the search to an exact match, we find just six segments.<jupyter_code>display_dark_matter(fish_dark_matter([awuhan6_seed], 1.0))<jupyter_output><empty_output><jupyter_text>Adding Species<jupyter_code>metadata = pd.read_csv('../data/CMS001_CMS002_MergedAnnotations_190325.csv')
metadata.head()
metadata.columns
species = {sample: species for sample, species in zip(metadata['NewIDseqName'], metadata['compute_species'])}
genus = {sample: species for sample, species in zip(metadata['NewIDseqName'], metadata['compute_genus'])}<jupyter_output><empty_output><jupyter_text>Likely, many of the Culex contigs are host which made it through filtering, because the genomes for the Culex species is much worse than that for Aedes.<jupyter_code>genus_combo = pd.Series(["-".join(set([genus[member.sample] for member in clusters[cluster_id]]))
                      for cluster_id in clusters]).value_counts()
genus_combo.plot.bar(color='C0')
species_combo = pd.Series(["-".join(set([species[member.sample] for member in clusters[cluster_id]]))
                      for cluster_id in clusters]).value_counts()
species_combo.plot.bar(color='C0', figsize=(20,5))<jupyter_output><empty_output><jupyter_text>Graph<jupyter_code>from cyjupyter import Cytoscape
from networkx.readwrite.json_graph import cytoscape_data
components = [G.subgraph(c) for c in nx.connected_components(G)]
subset = nbhd(wuhan6_exemplar_clusters, 1)
subg = G.subgraph(subset)
nx.draw(subg)
Cytoscape(data=cytoscape_data(subg))
def is_sample(node):
    return 1 if node.startswith('CMS') else 0
nx.draw(subg, node_size=50, node_color = [is_sample(node) for node in subg.nodes()])
Cytoscape(data=cytoscape_data(aw6_subgraph), layout_name='circle')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11202346041055719, 'fraction_numerical': 0.133841642228739, 'mean_word_length': 5.381736526946108, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '400080', 'license': 'MIT', 'n_tokens_mistral': 4151, 'n_tokens_neox': 3418, 'n_words': 537, 'path': 'analysis/darkmatter/Multisegmented Virus Hunting.ipynb', 'repo_name': 'wolfQK/california-mosquito-study', 'size': '8525'}