<jupyter_start><jupyter_text>For filtered dataset - sample filters on, genotype filters off:- generate amplicon (gene) trees for pseudohaploid samples- generate species tree reconciling gene treesManual steps: none<jupyter_code>%run common.ipynb<jupyter_output><empty_output><jupyter_text>Install and test astral<jupyter_code># astral not packaged in conda, install from github
# output jar
astral = '../../../bin/Astral/astral.5.6.3.jar'
# install
if not os.path.isfile(astral):
    ! wget https://github.com/smirarab/ASTRAL/raw/master/Astral.5.6.3.zip -P {BIN_DIR}
    ! unzip {os.path.join(BIN_DIR, 'Astral.5.6.3.zip')} -d {BIN_DIR}
    ! rm {os.path.join(BIN_DIR, '*.zip*')}
# test installation
! java -jar {astral} -i {os.path.dirname(astral)}/test_data/song_primates.424.gene.tre<jupyter_output>================== ASTRAL ===================== 

This is ASTRAL version 5.6.3
Gene trees are treated as unrooted
424 trees read from ../../../bin/Astral/test_data/song_primates.424.gene.tre
All output trees will be *arbitrarily* rooted at Marmoset

======== Running the main analysis
Number of taxa: 14 (14 species)
Taxa: [Marmoset, Orangutan, Human, Chimpanzee, Gorilla, Macaque, Galago, Mouse_Lemur, Tree_Shrew, Rat, Tarsier, Rabbit, Horse, Sloth]
Taxon occupancy: {Human=424, Rat=424, Tarsier=424, Galago=424, Rabbit=424, Macaque=424, Sloth=424, Marmoset=424, Tree_Shrew=424, Chimpanzee=424, Mouse_Lemur=424, Horse=424, Orangutan=424, Gorilla=424}
Number of gene trees: 424
0 trees have missing taxa
Calculating quartet distance matrix (for completion of X)
Species tree distances calculated ...
Building set of clusters (X) from gene trees 
------------------------------
gradient0: 339
Number of Clusters after addition by distance: 339
calculating extra bipartitions to be added at level 1 ..[...]<jupyter_text>Import and filter sequencing data<jupyter_code>seq_data = pd.read_csv(CLUSTERING, dtype={'target':str})
seq_data.columns
# filter sequencing data 
hq_seq_data = seq_data[(seq_data.species != 'unknown')] # &
#                        (~seq_data.outlier_genotype) & 
#                        (~seq_data.split_alleles)]
display(seq_data.shape, hq_seq_data.shape)<jupyter_output><empty_output><jupyter_text>Make pseudohaploid samplesGenerate two pseudohaploid sample IDs for each of the diploid sample. For homozygous samples, use same sequence for both alleles. For multiple allelic sequences, use first two sequences as these appear in `hq_seq_data`.<jupyter_code>def get_first_seq(df):
    return df.consensus.iloc[0]
first_seq = hq_seq_data.groupby(by=['s_Sample','target']).apply(get_first_seq).reset_index()
first_seq['species'] = hq_seq_data.groupby(by=['s_Sample','target'])['species'].max().reset_index(drop=True)
first_seq['hap_Sample'] = first_seq.s_Sample + '-1'
first_seq.head()
def get_second_seq(df):
    if df.shape[0] > 1:
        return df.consensus.iloc[1]
    else:
        return df.consensus.iloc[0]
second_seq = hq_seq_data.groupby(by=['s_Sample','target']).apply(get_second_seq).reset_index()
second_seq['species'] = hq_seq_data.groupby(by=['s_Sample','target'])['species'].max().reset_index(drop=True)
second_seq['hap_Sample'] = second_seq.s_Sample + '-2'
second_seq.head()
display(second_seq.shape)
comb_seq = pd.concat([first_seq, second_seq])
comb_seq.shape<jupyter_output><empty_output><jupyter_text>Species-sample mapping<jupyter_code># sequence-species mapping file
with open(MAPPING, 'w') as o:
    for _, r in comb_seq.groupby('species')['hap_Sample'] \
                      .unique() \
                      .reset_index() \
                      .iterrows():
        o.write('{}:{}\n'.format(r[0], ','.join(r[1])))
! head -1 {MAPPING}<jupyter_output>Anopheles_aconitus:VBS00053-1,VBS00055-1,VBS00053-2,VBS00055-2<jupyter_text>Alignment<jupyter_code># seequence export and alignment    
for ampl in AMPLS:
    sys.stdout.write('\rProcessing {}'.format(ampl))
    # subset amplicon data
    ampl_data = comb_seq[comb_seq.target == ampl]
    with open('temp.fa', 'w') as o:
        for (i, row) in ampl_data.iterrows():
            o.write('>{}\n{}\n'.format(row.hap_Sample,
                                       row[0]))
    ! mafft temp.fa > {ALN_SP_TREE.format(ampl)} 2> /dev/null
    ! rm temp.fa
sys.stdout.write('\nDone!\n')<jupyter_output>Processing 61
Done!<jupyter_text>Gene trees<jupyter_code># write all trees into single file
! rm -f {GENE_TREES}
for ampl in AMPLS:
    sys.stdout.write('\rProcessing {}'.format(ampl))
    ! FastTree -nt {ALN_SP_TREE.format(ampl)} >> {GENE_TREES} 2> /dev/null
! wc -l {GENE_TREES}<jupyter_output>Processing 61      62 data/5_gene_trees.nwk<jupyter_text>Species tree<jupyter_code>! java -jar {astral} -i {GENE_TREES} -o {SPECIES_TREE} -a {MAPPING}<jupyter_output>================== ASTRAL ===================== 

This is ASTRAL version 5.6.3
Gene trees are treated as unrooted
62 trees read from data/5_gene_trees.nwk
All output trees will be *arbitrarily* rooted at Anopheles_atroparvus

======== Running the main analysis
Number of taxa: 306 (56 species)
Taxa: {Anopheles_sinensis=[anopheles-sinensis-sinensisscaffoldsasins2-1, anopheles-sinensis-sinensisscaffoldsasins2-2, anopheles-sinensis-chinascaffoldsasinc2-1, anopheles-sinensis-chinascaffoldsasinc2-2], Anopheles_oryzalimnetes=[Aory-SP141-2, Aory-SP141-1], Anopheles_vagus=[VBS00158-1, VBS00158-2, VBS00157-2, VBS00157-1, VBS00156-2, VBS00156-1, VBS00154-1, VBS00154-2], Anopheles_farauti=[anopheles-farauti-far1scaffoldsafarf2-1, anopheles-farauti-far1scaffoldsafarf2-2], Anopheles_maculatus_A=[VBS00108-2, VBS00108-1, VBS00107-2, VBS00107-1], Anopheles_maculatus_B=[anopheles-maculatus-maculatus3scaffoldsamacm1-1, anopheles-maculatus-maculatus3scaffoldsamacm1-2], Anopheles_tenebrosus=[Aten-185-2, A[...]<jupyter_text>Sample tree<jupyter_code>with open(SAMPLE_MAPPING, 'w') as o:
    for _, r in comb_seq.groupby('s_Sample')['hap_Sample'] \
                      .unique() \
                      .reset_index() \
                      .iterrows():
        
        o.write('{}:{}\n'.format(r[0], ','.join(r[1])))
! head -1 {SAMPLE_MAPPING}
! java -jar {astral} -i {GENE_TREES} -o {SAMPLE_TREE} -a {SAMPLE_MAPPING}<jupyter_output>================== ASTRAL ===================== 

This is ASTRAL version 5.6.3
Gene trees are treated as unrooted
62 trees read from data/5_gene_trees.nwk
All output trees will be *arbitrarily* rooted at VBS00112

======== Running the main analysis
Number of taxa: 306 (153 species)
Taxa: {Aimp-M0001=[Aimp-M0001-1, Aimp-M0001-2], Avin-B0009=[Avin-B0009-1, Avin-B0009-2], Aten-191=[Aten-191-2, Aten-191-1], anopheles-farauti-far1scaffoldsafarf2=[anopheles-farauti-far1scaffoldsafarf2-1, anopheles-farauti-far1scaffoldsafarf2-2], Athe-6-2=[Athe-6-2-2, Athe-6-2-1], Abro-22=[Abro-22-2, Abro-22-1], Acol-570=[Acol-570-1, Acol-570-2], Athe-6-1=[Athe-6-1-2, Athe-6-1-1], Abro-21=[Abro-21-2, Abro-21-1], Anils-7=[Anils-7-1, Anils-7-2], Avin-M0012=[Avin-M0012-2, Avin-M0012-1], anopheles-gambiae-pimperenascaffoldsagams1=[anopheles-gambiae-pimperenascaffoldsagams1-1, anopheles-gambiae-pimperenascaffoldsagams1-2], anopheles-funestus-fumozchromosomesafunf3=[anopheles-funestus-fumozchromosomesafunf3-1, ano[...]
{'fraction_non_alphanumeric': 0.15938803894297635, 'fraction_numerical': 0.0521557719054242, 'mean_word_length': 4.609204368174727, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 38, 'https://': 1, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '126674', 'license': 'MIT', 'n_tokens_mistral': 3016, 'n_tokens_neox': 2714, 'n_words': 512, 'path': 'work/7_species_id/5_species_tree.ipynb', 'repo_name': 'malariagen/ano-spp', 'size': '7190'}