<jupyter_start><jupyter_text>Read in data<jupyter_code>gitdir = "/Users/Hanna/Desktop/MyBox/aa_DeRisi/mosquito/california-mosquito-study/"
identity_qcov_threshold = 0.9
## Load metadata
metadata = pd.read_csv(gitdir + "/data/metadata/CMS001_CMS002_MergedAnnotations.csv", header=0)
idseq_data = pd.read_csv(gitdir + "/data/metadata/idseq_metadata.csv", header=0)
idseq_data = idseq_data[~idseq_data['sample'].str.contains('Placeholder')]
## Load contig data
contig_calls_all = pd.read_csv(gitdir + "/data/s3/contig_calls.tsv", sep="\t", header=0)
contig_calls_decontam = pd.read_csv(gitdir + "/data/s3/contig_calls_decontam.tsv", sep="\t", header=0)
## Calculate nonhost reads taking hexapoda and contaminating reads into account
metadata_cols = ["NewIDseqName","ska_genus", "ska_species", "collected_by"]
idseq_cols = ["sample","nonhost_reads", "total_reads"]
# Calculate reads removed during decontamination and hexapoda subtraction
contam_counts = pd.DataFrame(contig_calls_all.groupby('sample').read_count.sum() - 
                             contig_calls_decontam.groupby('sample').read_count.sum() )
# Sum reads to hexapoda contigs by sample
contigs_hex = contig_calls_all[(contig_calls_all.hexapoda == True) & (contig_calls_all.curated == False)]
hex_counts = contigs_hex.groupby("sample")["read_count"].sum().reset_index()
# Collect metadata, idseq nonhost reads, and additional hexapoda reads
metadata = pd.concat([ metadata[metadata_cols].rename(columns={"NewIDseqName":"sample"}).set_index("sample"), 
           idseq_data[idseq_cols].rename(columns={"nonhost_reads":"idseq_nonhost_reads"}).set_index("sample"), 
           hex_counts.rename(columns={"read_count":"hex_read_count"}).set_index("sample"),
           contam_counts.rename(columns={'read_count': 'contaminating_and_hex_reads'})], 
                     axis=1, sort=True)
#Check that hexapoda counts are < all contaminating_and_hex counts for all samples
#print(metadata[~(metadata['contaminating_and_hex_reads'] >= metadata['hex_read_count'])]) #shld only contain waters
# Subtract contaminating+hex reads from idseq_nonhost_reads to get final nonhost reads
metadata['hex_read_count'].fillna(0,inplace=True)
metadata['contaminating_and_hex_reads'].fillna(0,inplace=True)
metadata['nonhost_reads'] = metadata['idseq_nonhost_reads'] - metadata['contaminating_and_hex_reads']
# Reorganizing columns
metadata['sample'] = metadata.index
metadata.reset_index(inplace=True)
metadata.drop(columns='index',inplace=True)
# Write final nonhost reads to file
metadata[['sample','nonhost_reads']].to_csv(
    gitdir + "/data/metadata/nonhost_reads_decontam_nohexapoda.tsv", index=False, sep="\t")
# Update contig stats with nonhost reads
contig_stats = pd.merge(contig_calls_decontam, metadata, how="left", on="sample")
# Check no hexapoda remain
print( contig_stats[(contig_stats.hexapoda == True) & (contig_stats.curated == False)].shape[0] == 0 ) #shld be True
# Calculate read proportion and identity/query-coverage threshold
contig_stats['read_prop'] = contig_stats['read_count'] / contig_stats['nonhost_reads']
# Calculate identity/query-coverage threshold
contig_stats = contig_stats.assign(identity_qcov=(contig_stats["identity"]/100 
                                                  * contig_stats["align_length"]/contig_stats["contig_length"]
                                                 ).apply(lambda x: min(x, 1)))
contig_stats.loc[contig_stats['curated'] == True, 'identity_qcov'] = np.nan<jupyter_output>True<jupyter_text>Viruses<jupyter_code>contig_stats_viral = contig_stats[(contig_stats["taxon_group"]=="Viruses") & contig_stats["curated"]==True]
# Information about Baltimore classification of virus family groups
viral_family_groups = pd.read_csv(gitdir + "/data/virus_family_groups.csv", header=0)
viral_family_groups = viral_family_groups.loc[:, ~viral_family_groups.columns.str.startswith('Unnamed')]
viral_family_groups = viral_family_groups.append(
    {"family":"Chuviridae",
    "baltimore_group":"V",
    "genome_description":"",
    "number_of_segments_conventional":""}
    , ignore_index = True)

# Convert virus json into data frame
with open (gitdir +"/data/darkmatter/virus.json", "r") as f:
    viral_json = pd.DataFrame(json.load(f)).T
    viral_json['poly_group'] = viral_json.index

# Check which viral families not represented
print(set(viral_json.family.unique()) - set (viral_family_groups.family.unique())) #should be empty

# Merge viral info with contig stats
viral_info = viral_json.merge(viral_family_groups,on='family')
viral_info.provisional_name.fillna(viral_info.name, inplace=True)
viral_info.rename(columns={'provisional_name':'sci_name'},inplace=True)
contigs_viral = contig_stats_viral.merge(viral_info[['poly_group','family','baltimore_group',
                                                     'genome_description','number_of_segments_conventional',
                                                     'sci_name']].
                                              astype({'poly_group': 'float64'}),
                                              on='poly_group', how = 'left')
contigs_viral = contigs_viral.assign(group="Virus")<jupyter_output>set()<jupyter_text>Non-viral contigs<jupyter_code>def get_rows_taxid(df, taxid, taxid_colname="taxid"):
    if (not isinstance(df, pd.DataFrame)):
        if isinstance(taxid, str):
            taxid = ncbi.get_name_translator([taxid])[taxid][0]
        return (taxid in ncbi.get_lineage(df))
    outdf = df[df[taxid_colname].apply(get_rows_taxid, taxid=taxid)]
    return (outdf) 

def group_at_higher_tax(df, taxonomic_group, family_name, taxid_colname="taxid", family_colname="family"):
    groups = {}
    for x in taxonomic_group:
        if isinstance(x, str):
            groups[x] = ncbi.get_name_translator([x])[x][0]
        else:
            groups[ncbi.get_taxid_translator(x)[x]] = x
    family_assignments = {}
    for x in df["taxid"].unique():
        lin = ncbi.get_lineage(x)
        family_assignments[x] = family_name
        for key, taxid_x in groups.items():
            if taxid_x in lin:
                family_assignments[x] = key
                break
    df[family_colname] = df["taxid"].apply(lambda x: family_assignments[x])
    return (df)

def get_summary_table (df, colnames, metric):
    df = df.groupby(colnames)[metric].sum().reset_index()
    if not isinstance(metric, list):
        metric = [metric]
    sort_order = colnames+metric
    sort_order.remove("sample")
    return (df.sort_values(by=sort_order, ascending=False))
# Subset to contigs that are not viral, and are strong NCBI hits based on identity/qcov threshold
contig_stats_nonviral = contig_stats[ (contig_stats["taxon_group"]!="Viruses") &
                                    (contig_stats["identity_qcov"]>=identity_qcov_threshold) ]
## WOLBACHIA
wolbachia_taxid = 952
wolbachia_contigs = get_rows_taxid(contig_stats_nonviral, wolbachia_taxid);
## METAZOANS
metazoan_contigs = contig_stats_nonviral[contig_stats_nonviral["taxon_group"]=="Metazoa"];
# In the 'family' column, add information about the higher taxonomic grouping for taxids
metazoan_contigs = group_at_higher_tax(metazoan_contigs, 
                                       taxonomic_group=["Nematoda"], 
                                       family_name="Other metazoa", taxid_colname="taxid", family_colname="family")
## CHORDATA
chordate_contigs = get_rows_taxid(contig_stats_nonviral, taxid="Chordata");
# In the 'family' column, add information about the higher taxonomic grouping for taxids
chordate_contigs = group_at_higher_tax(chordate_contigs, 
                                       taxonomic_group=["Pecora", "Aves", "Carnivora", "Rodentia", "Leporidae"], 
                                       family_name="Other chordates", taxid_colname="taxid", family_colname="family")
## OTHER EUKARYOTES
eukaryotic_contigs = contig_stats_nonviral[contig_stats_nonviral["taxon_group"]=="Eukaryota"];
# In the 'family' column, add information about the higher taxonomic grouping for taxids
eukaryotic_contigs = group_at_higher_tax(eukaryotic_contigs, 
                                         taxonomic_group=["Trypanosomatidae", "Apicomplexa", "Microsporidia"], 
                                         family_name="Other eukaryotes", taxid_colname="taxid", family_colname="family")
# Combine nonviral
contigs_nonviral = pd.concat([wolbachia_contigs.assign(group="Wolbachia").assign(family="Wolbachia"), 
                            metazoan_contigs.assign(group="Metazoa"),
                            chordate_contigs.assign(group="Chordates"),
                            eukaryotic_contigs.assign(group="Other Eukaryotes")], sort=False)
# Convert taxids to human_readable
contigs_nonviral = contigs_nonviral.assign(sci_name=contigs_nonviral["taxid"].
                                           apply(lambda x: ncbi.get_taxid_translator([x])[x]))
contigs_nonviral.loc[contigs_nonviral['group']=='Wolbachia', ['sci_name']] = 'Wolbachieae'<jupyter_output><empty_output><jupyter_text>Merging nonviral and viral<jupyter_code>contigs_out = pd.concat([contigs_nonviral, contigs_viral], sort=False)
# Count total number of contigs in each broad group
numbers = contigs_out.groupby("group")["group"].count()
print(numbers)
# Summary table to combine reads to all contigs from the same sample-taxon pair
summary_nonviral = get_summary_table(contigs_out[contigs_out["group"]!="Virus"], 
                            colnames=["ska_genus","ska_species","collected_by","sample","taxid","group","family","sci_name"],
                                     metric=["read_count", "read_prop"])
summary_viral = get_summary_table(contigs_out[contigs_out["group"]=="Virus"], 
                            colnames=["ska_genus","ska_species","collected_by","sample","taxid","group","family","sci_name",
                                      "poly_group","baltimore_group","genome_description","number_of_segments_conventional"],
                            metric=["read_count", "read_prop"])
contigs_out_summary = pd.concat([summary_nonviral, summary_viral], sort=False, ignore_index = True)
# Rename column for consistency
contigs_out_summary.rename(columns={"read_count":"reads"}, inplace = True)
# Check all reads are accounted for in the summarizing
contigs_out['read_count'].sum() == contigs_out_summary['reads'].sum() #Shld be True
# Write to file
contigs_out_summary.to_csv(gitdir + "/figures/fig3/all_contigs_df_new.tsv", sep="\t", index=False)
# Write contig names to file for eukaryotes-in-bloodmeals analysis
category_names = ['Apicomplexa','Trypanosomatidae']
for catname in category_names:
    df = contigs_out[(contigs_out['family']==catname) & (contigs_out['sample'].str.contains('CMS001'))]
    out = df['sample'] + "~" + df['contig_name']
    filename = "contig_names_CMS001_" + catname + ".txt"
    print(catname + ": " + str(len(out)) + " contigs")
    out.to_csv(gitdir + "/figures/fig4/" + filename, sep="\n", index=False)<jupyter_output>Apicomplexa: 9 contigs
Trypanosomatidae: 71 contigs<jupyter_text>Comparing old vs. new<jupyter_code>def whatsin(df):
    return [{x:df[x].unique()} for x in df.columns.tolist()]

def get_desc_contigs(contig_df, txid):
    taxa = ncbi.get_descendant_taxa(txid, intermediate_nodes=True)
    taxa.append(txid)
    df = contig_df[contig_df["taxid"].apply(lambda x: x in taxa)]
    return df
# Checking column names and values in old vs. new
old_df = pd.read_csv(gitdir + "/figures/fig3/all_contigs_df.tsv", sep="\t", header=0)
for mycol in old_df.columns:
    print("For column: " + str(mycol))
    a = set(old_df[mycol].unique())
    b = set(contigs_out_summary[mycol].unique())
    print("old df uniquely contains: " + str(a-b))
    print("new df uniquely contains: " + str(b-a))
    print("------")
txid = 5654
category_name = 'Trypanosomatidae'
new_level = 'family'
old_level = 'family'

def compare_new_to_old(txid, category_name, new_level = 'family', old_level = 'family'):
    df = pd.concat([ get_desc_contigs(contigs_out,txid).groupby('sample')['read_prop'].sum().rename('new_contigs_by_taxid'),
         contigs_out[(contigs_out[new_level]==category_name)].groupby('sample')['read_prop'].sum().rename('new_contigs_by_category'),
         contigs_out_summary[contigs_out_summary[new_level]==category_name].groupby('sample')['read_prop'].sum().rename('new_summary_by_category'),
         old_df[old_df[old_level]==category_name].groupby('sample')['read_prop'].sum().rename('old_by_category') ],
         axis=1, sort=False)
    return df
compare_new_to_old(5654,'Trypanosomatidae')
compare_new_to_old(5794,'Apicomplexa')
compare_new_to_old(6231,'Nematoda')
compare_new_to_old(629,'Onchocercidae')
compare_new_to_old(6029,'Microsporidia')
contigs_out[contigs_out.taxon_group == "Viruses"].groupby("family")["family"].count()
contigs_out[contigs_out.taxon_group != "Viruses"].groupby("family")["family"].count()
contigs_out[contigs_out.taxon_group != "Viruses"].groupby("group")["group"].count()
contigs_out.family.unique()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1326530612244898, 'fraction_numerical': 0.005217124443762468, 'mean_word_length': 4.50464527027027, 'pattern_counts': {'":': 9, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '77041', 'license': 'MIT', 'n_tokens_mistral': 4429, 'n_tokens_neox': 4185, 'n_words': 685, 'path': 'scripts/process_contigs_singleton_fig.ipynb', 'repo_name': 'wolfQK/california-mosquito-study', 'size': '13034'}