<jupyter_start><jupyter_text>This script downloads metagenomic datasets useds by Pasoli et al. 2019 from their lab repository. This fastas unzipped take up roughly 2.3T. gyrb reads are filtered out using a blastn search by a fasta containing ASVs from the amplicon datasets and reads from the gyrb gtdbtk reference fasta. Reads are output into an ASV table where the columns represent the sequences and the rows represent sample names, and the values represent the number of samples that particular read was found. This enables the metagenomic data to be merged with the amplicon data. Although it should be noted that the read numbers between the two types of datasets are incomparable, and should be used to indicate presene/absence in various sample types.Inputs: list of metagenomic datasets https, Bt-ASV fastaOutput: ASV table of metagenomic reads to be merged with amplicon datasetsDependencies: Biopython(v1.77), BLAST(v2.9.0)<jupyter_code>import os
import pandas as pd
import glob
import numpy as np
from Bio import SeqIO
from Bio import Seq<jupyter_output><empty_output><jupyter_text>Overview of inputs/outpusdue to the size of the data file only the results files will be permanently saved.<jupyter_code>#root_dir = '/Volumes/APE_MacPro_External_2/alex/captive_ape_microbiome'
root_dir = '/Volumes/AHN/captive_ape_microbiome/'
os.chdir(root_dir)

analysis_dir = 'results/gyrb/processing/metagenomic_samples/'
os.chdir(analysis_dir)

metagenomic_https_file = 'metadata_http_to_dnload.txt'
metagenomic_samples_file = 'metagenomic_metadata_all_samples.txt'
ASV_fasta =  'Bt_ASV_gtdbtk.fasta'

os.system('mkdir data') #
os.system('mkdir sample_files')
os.system('mkdir blastout')
os.system('mkdir parse_blast')
os.system('mkdir sseq_adj')
os.system('mkdir final_outputs')
#generate table of metagenomic datasets 
metagenomic_datasets = pd.read_csv(metagenomic_https_file,header=None)
metagenomic_datasets.columns = ['https']
metagenomic_datasets['filename'] = metagenomic_datasets['https'].apply(lambda x: x.split('/')[-1])
metagenomic_datasets['folder'] = metagenomic_datasets['filename'].apply(lambda x: x.split('.')[0])
metagenomic_datasets.head()
metagenomic_datasets.to_csv('metadata_datasets.txt',sep='\t',index=False)
metagenomic_datasets.head()<jupyter_output><empty_output><jupyter_text>Download datasets<jupyter_code>def dnload_file(http): 
    """
    downloads and expands metagenomic dataset from Passoli repo, generates lists of sample files in the data folder
    """
    filename = http.split('/')[-1]
    folder = filename.split('.')[0]
    print(http,filename,folder)
    
    os.system(f'wget {http} --directory-prefix data/')
    os.system(f'tar -vxf data/{filename} -C data/ ')
    os.system(f'ls data/{folder}/*.fa > sample_files/{folder}_files.txt')  

#test download for one dataset
AsnicarF_http = metagenomic_datasets['https'].iloc[0]
dnload_file(AsnicarF_http)
#list datasets already downloaded
datasets_downloaded = [f.split('.')[0] for f in os.listdir('data/') if f.endswith('bz2')]
print('datasets already downloaded')
print(datasets_downloaded)

#determine datasets that need to be downloads
datasets_to_do = metagenomic_datasets[~metagenomic_datasets['folder'].isin(datasets_downloaded)]
print('\ndatasets left to download')
print(datasets_to_do)
for http in datasets_to_do['https']:
    print(http)
    #dnload_file(http) #uncomment to download<jupyter_output>datasets already downloaded
['AsnicarF_2017']

datasets left to download
                                                https  \
1   https://www.dropbox.com/s/n6867hhmo12xnwf/Back...   
2   https://www.dropbox.com/s/ac6cpq50pt5mqj6/Beng...   
3   https://www.dropbox.com/s/uj9cy9914l0cq3x/Brit...   
4   https://www.dropbox.com/s/dtxmz1apdgtf3d1/Cast...   
5   https://www.dropbox.com/s/0469lqd8m779nxe/Chen...   
6   https://www.dropbox.com/s/whsav7tbeefgjww/Chng...   
7   https://www.dropbox.com/s/desbb91o2pq485k/CM_c...   
8   https://www.dropbox.com/s/plett88rfu15v7x/CM_m...   
9   https://www.dropbox.com/s/ly8ncdsw9y3wbln/CM_p...   
10  https://www.dropbox.com/s/qo0zudiqgelkqys/Cost...   
11  https://www.dropbox.com/s/azfhfhfgq2tshv6/Davi...   
12  https://www.dropbox.com/s/35i74mbblc6rivy/Feng...   
13  https://www.dropbox.com/s/3ignf2yinb4lqp7/Ferr...   
14  https://www.dropbox.com/s/42vc198r249t76u/Geve...   
15  https://www.dropbox.com/s/9l9l6meb4fq7zic/Hann...   
16  https://www[...]<jupyter_text>helper scripts big to small for filtering gyrb reads<jupyter_code>def parse_blast(blast_results_file):
    """
    convert raw blast results in a table, selecting only one (the best) hit per metagenomic read
    this blastn table is used by future scripts (prep_df and ) to find the metagenomic fasta, determine the read overlap, 
    """
    #this clause checks to see whether the blastn search of a given metagenomic sample fasta returned any hits
    try:
        blast_results=pd.read_csv(blast_results_file, sep='\t', comment='#',header=None)     
    except:
        print(blast_results_file,'no hits')
        return
    
    qlen_cutoff = .8 #metagenomic reads must have 80% alignment length
    blast_results=pd.read_csv(blast_results_file, sep='\t', comment='#',header=None)
    blast_results.columns=['ASV','sacc','pident','length','qlen','qstart','qend','sstart','send','evalue','sseq']   
    blast_results=blast_results[blast_results['length']/blast_results['qlen'] >= qlen_cutoff]
    blast_results=blast_results.sort_values('evalue')  
    blast_results=blast_results.groupby('sacc').head(1) #get top hit

    sample_name=blast_results_file.split('/')[-1].split('_Bt.txt')[0]
    blast_results['sample_name']=sample_name
    folder=blast_results_file.split('/')[1]
    blast_results['folder']=folder
    outfile = blast_results_file.replace('blastout','parse_blast')
    blast_results.to_csv(outfile,sep='\t',index=False)
    return(blast_results)


def prep_df(df):  
    """
    blastn may return a hit that is not the full length of the ASV. So one can't just take the subject sequence
    This scripts determine whether the blast hit spans the full length of the ASV. If it doesn't, the read is labeled as long or short
    Long hit means that the ASV had gaps in the alignment, Short hit means that the alignment score fell at the ends of the read
    Additionally, hits may be need to reverse complemented
    
    This script determine the metagenomic reads starting position based on where the subject/metagenomic reads aligned with the 
    ASV/query reads. i.e. if ASV start = 3 and the metagenomic reads start is 250000, the adjusted metagenomic start is  250000-3. 
    This is complicated slightly if the read is in the rc orientation, then the ASV start is added to the metagenomic start.
    
    The final step of this script is to return upstream and downsteam positions that can be used to excise the hit from the metagenomic fasta 
    regardless of read orientation.
    """
    
    df['overlap'] = np.where(df['qlen']==df['length'], 'exact', #or if not exact must be long or short
                             np.where(df['qlen']<df['length'], 'long','short'))
    df['orientation'] = np.where(df['send']>df['sstart'], 'in_frame', 'reverse_complement')
    df['sstart_adj'] = np.where(df['orientation']=='in_frame',
                                             df['sstart'] - df['qstart'],
                                             df['sstart'] + df['qstart'] - 1)
    df['send_adj'] = np.where(df['orientation']=='in_frame',
                                             df['sstart_adj'] + df['qlen'],
                                             df['sstart_adj'] - df['qlen'])
    df['upstream'] = df.apply(lambda row:
                                  str(min([int(row['sstart_adj']),int(row['send_adj'])])),
                                  axis=1)
    df['downstream'] = df.apply(lambda row:
                                  str(max([int(row['sstart_adj']),int(row['send_adj'])])),
                                  axis=1)
    return(df)

def fetch_sseq_adj(row):
    """
    This script takes in a row of the blastn hit table after its been modified by prep_df.
    It determines whether blastn hit was full-length or whether is was shorter or longer than the ASV seqs
    It checks to see whether the adjusted metagenomic read start is still in bounds. For instance, if a hit 
    has a metagenomic read start of 10 but a query/ASV start position of 30 this would be out of bounds, indicated
    by a negative upstream value.
    """
    downstream = int(row['downstream'])
    upstream = int(row['upstream'])
    
    if row['length']==row['qlen']: #exact seq match
        return(pd.Series(['NA',row['sseq']]))
    elif upstream < 0:
        return(pd.Series(['NA','outOfBounds']))
    else:
        #go fishing 
        try:
            fasta_path = row['folder']+'/'+row['sample_name']+'.fa'
            seq_dict = SeqIO.to_dict(SeqIO.parse('data/'+fasta_path,'fasta'))
            contig = seq_dict[row['sacc']].seq
            len_cont = len(str(contig))
            if downstream >len_cont: #means contig is too short, hit goes over end
                return(pd.Series([len_cont,'outOfBounds']))
            else:
                sseq_adj = str(seq_dict[row['sacc']].seq)[upstream:downstream]
                if row['orientation']=='reverse_complement':
                    sseq_adj = Seq.reverse_complement(sseq_adj)     
                return(pd.Series([len_cont,sseq_adj]))
        except:
            return('NA','contig_not_found')   
    
def filter_gyrb_reads(folder):     
    """
    the big kahuna. It returns a table with all the gyrb hits for all sample fastas within a folder.
    
    this script puts it all together. Given is a dataset folder that needs to have gyrb reads exacted
    First, all fastas within the folder are listed, then one-by-one fastas are made into blastdbs and searched with the 
    ASV.fasta. Results from on the samples within a folder are parsed (parse_blast function) concatenated into one dataframe.
    To determine which reads are full length and/or reverse complement the dataframe is processed the with prep_df function.
    fetch_sseq_adj then uses the info to identify adjusted reads that can be pulled from the metagenomic fasta or 
    whether these reads are out of bounds. 
    """
    
    fastas = [f for f in os.listdir('data/'+folder) if f.endswith('.fa')]
    print(fastas)
    for fasta in fastas:
        fasta_name = fasta.split('.fa')[0]
        os.system(f'makeblastdb -in data/{folder}/{fasta} -dbtype nucl')
        os.system(f'blastn -query {ASV_fasta} -db data/{folder}/{fasta}  \
                    -outfmt "7 qseqid sseqid pident length qlen qstart qend sstart send evalue sseq" \
                    -out blastout/{folder}/{fasta_name}_Bt.txt')
    
    all_files = glob.glob(f"blastout/{folder}/*.txt")   
    df_from_each_file = [parse_blast(f) for f in all_files]
    df = pd.concat(df_from_each_file, ignore_index=True)   
    
    df = prep_df(df)
    df[['contig_len','sseq_adj']] = df.apply(lambda row: fetch_sseq_adj(row),axis=1)
    df['sseq_same']= df['sseq']==df['sseq_adj']
    df['sseq_adj_len'] = df['sseq_adj'].apply(lambda x: len(str(x)))
    df = df[['ASV','sacc','sample_name','folder','pident','sseq_adj_len','evalue','overlap','orientation','qstart','qend','sstart','send','sstart_adj','send_adj','upstream','downstream','sseq_same','sseq','sseq_adj']]
    print(len(df),'gyrb hits in dataset',folder)
    df.to_csv('sseq_adj/'+folder+'.txt',sep='\t',index=False)

#test filter gyrb
filter_gyrb_reads('AsnicarF_2017')<jupyter_output>['AsnicarF_2017__MV_FEI1_t1Q14.fa', 'AsnicarF_2017__MV_FEI2_t1Q14.fa', 'AsnicarF_2017__MV_FEI3_t1Q14.fa', 'AsnicarF_2017__MV_FEI4_t1Q14.fa', 'AsnicarF_2017__MV_FEI4_t2Q15.fa', 'AsnicarF_2017__MV_FEI5_t1Q14.fa', 'AsnicarF_2017__MV_FEI5_t2Q14.fa', 'AsnicarF_2017__MV_FEI5_t3Q15.fa', 'AsnicarF_2017__MV_FEM1_t1Q14.fa', 'AsnicarF_2017__MV_FEM2_t1Q14.fa', 'AsnicarF_2017__MV_FEM3_t1Q14.fa', 'AsnicarF_2017__MV_FEM4_t1Q14.fa', 'AsnicarF_2017__MV_FEM4_t2Q15.fa', 'AsnicarF_2017__MV_FEM5_t1Q14.fa', 'AsnicarF_2017__MV_FEM5_t2Q14.fa', 'AsnicarF_2017__MV_FEM5_t3Q15.fa', 'AsnicarF_2017__MV_MIM2_t1M14.fa', 'AsnicarF_2017__MV_MIM3_t1M14.fa', 'AsnicarF_2017__MV_MIM4_t2F15.fa', 'AsnicarF_2017__MV_MIM5_t2M14.fa', 'AsnicarF_2017__MV_MIM5_t3F15.fa']
blastout/AsnicarF_2017/AsnicarF_2017__MV_MIM2_t1M14_Bt.txt no hits
blastout/AsnicarF_2017/AsnicarF_2017__MV_MIM3_t1M14_Bt.txt no hits
blastout/AsnicarF_2017/AsnicarF_2017__MV_MIM4_t2F15_Bt.txt no hits
blastout/AsnicarF_2017/AsnicarF_2017__MV_MIM5_t2M14_Bt.txt no h[...]<jupyter_text>Filter metagenomic fastas for gyrb hits<jupyter_code>#list datasets that have already been filtered
datasets_all = list(metagenomic_datasets['folder'])

datasets_complete = [f.split('.')[0] for f in os.listdir('sseq_adj/')]
print('\ndatasets already filtered')
print(datasets_complete)

datasets_to_do = set(datasets_all) - set(datasets_complete)
print('\ndatasets left to filter')
print(datasets_to_do)
#prep for parallel processing
import multiprocessing as mp
print("Number of processors: ", mp.cpu_count())

#filter out gyrb seqs
pool = mp.Pool(mp.cpu_count())  
results = pool.map_async(filter_gyrb_reads, [folder for folder in datasets_to_do])
pool.close()<jupyter_output>Number of processors:  8<jupyter_text>Merge gyrb hits across datasets and add sample metadata<jupyter_code>#merge results from complete datasets
dataset_gyrb_hits =  glob.glob(f"sseq_adj/*.txt")   
df_from_each_file = [pd.read_csv(f,sep='\t') for f in dataset_gyrb_hits]
datasets_merged = pd.concat(df_from_each_file, ignore_index=True)  
datasets_merged.to_csv('final_outputs/gyrb_hits_all_samples.txt',sep='\t',index=False)

#filter out reads that are not full length or have ambig character
datasets_merged = datasets_merged[datasets_merged['sseq_adj']!='outOfBounds']
datasets_merged = datasets_merged[~datasets_merged['sseq_adj'].str.contains('-')]

#print(datasets_merged.head())

#load metadata for samples
metadata = pd.read_csv(metagenomic_samples_file,sep='\t')
#print(metadata.head())
metadata['sample_name'] = metadata['Study']+'__'+metadata['Sample ID']

print("are there any samples that aren't in the metadata file?")
print(set(datasets_merged['sample_name'])-set(metadata['sample_name']))

print("how many samples returned gyrb hits?")
print(len(set(datasets_merged['sample_name'])))

df = datasets_merged.merge(metadata,on='sample_name')
df['Westernized_edit']=df['Westernized'].apply(lambda x: x.replace('Yes','western').replace('No','nonwestern'))
df['sample_group'] = df['folder']+'__'+df['Westernized_edit']
df['count'] = 1 #read detect in that one sample

df_select = df.groupby(['sample_name','sseq_adj'])['count'].sum().reset_index()
seqtab = df_select.pivot(index='sample_name',columns='sseq_adj',values='count').fillna(0)  
seqtab.to_csv('final_outputs/gyrb_metagenomic_seqtab.txt',sep='\t')

print('how many unique gyrb seqs were filtered out of metagenomic samples?')
print(len(seqtab.columns))

print('gives the number of samples each read is found in by percentile')
print(seqtab.sum(axis=0).quantile([.10,.50,.70,.80,.90,.95,.99]))<jupyter_output>are there any samples that aren't in the metadata file?
set()
how many samples returned gyrb hits?
7300
how many unique gyrb seqs were filtered out of metagenomic samples?
7124
gives the number of samples each read is found in by percentile
0.10      1.00
0.50      1.00
0.70      2.00
0.80      3.00
0.90      7.00
0.95     17.00
0.99    143.31
dtype: float64<jupyter_text>Output fasta of unqiue metagenomic readsuse first sample as header<jupyter_code>df_unique = df.groupby('sseq_adj').head(1).reset_index()
with open('final_outputs/gyrb_metagenomic_unique.fasta' ,'w') as f:
    for index,row in df_unique.iterrows():
        header='>'+str(row['sample_name'])+'\n'
        seq=row['sseq_adj']+'\n'
        f.write(header)
        f.write(seq)<jupyter_output>7124
StarcoderdataJupyterStructuredCleanDedup