<jupyter_start><jupyter_text>Blast markers versus dipterian genomes: Anopheles gambiae, Aedes aegyptii, Culex pipiens, Drosophila melanogaster.Results:- 591 of 663 markers with unique location in agam genome left for analysisMinor results:- add unaligned species to metadata- reorder and rename columns  Libraries and parameters<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pybedtools
from Bio import AlignIO
from Bio.Blast.Applications import NcbiblastnCommandline

DATA_DIR = "../../../data/"
# blast db from Anopheles gambiae, Aedes segypti, Culex pipiens, Drosophila melanogaster genomes
COMB_DB = DATA_DIR + 'genome/AgAaCpDm' 
# out
AMPL_FILE = "data/20180619_comb_ampl.maf"
META_FILE = "data/20180619_comb_ampl_data.csv"
BLAST_FILE = "data/20180706_blast_4sp.xml"
UPD_META_FILE = "data/20180706_comb_ampl_data.xlsx"
UPD_AMPL_FILE = "data/20180706_comb_ampl.maf"<jupyter_output><empty_output><jupyter_text>BLAST run<jupyter_code># extract agam sequences for markers removing alignment gaps
# headers are 0-based numbers of markers
agam_fa = 'data/temp_ampl.fasta'

with open('data/temp_ampl.fasta', 'w') as o:
    i = 0
    ampl_seq = AlignIO.parse(AMPL_FILE, "maf")
    for a in ampl_seq:
        o.write(">{}\n".format(i))
        o.write(str(a[0].seq).replace('-',''))
        o.write('\n')
        i += 1
    
!head {agam_fa}
help(NcbiblastnCommandline)
# generate blast command
cline = NcbiblastnCommandline(cmd='blastn', out=BLAST_FILE, outfmt=5, query=agam_fa, 
                              db=COMB_DB, evalue=0.001)
print(cline)
# run blast command. Return values are stdout and stderr
cline()<jupyter_output><empty_output><jupyter_text>BLAST parse<jupyter_code># parse hits with NCBIXML
# http://biopython.org/DIST/docs/tutorial/Tutorial.html#fig:blastrecord
from Bio.Blast import NCBIXML
blast_records = list(NCBIXML.parse(open(BLAST_FILE)))
# http://biopython.org/DIST/docs/api/Bio.Blast.Record.HSP-class.html
test_hsp = blast_records[0].alignments[1].hsps[0] 
print(str(test_hsp))
print(test_hsp.identities, test_hsp.align_length)
# extract nececcary blast data
from collections import defaultdict
bd = dict()
species = ('AgamP3', 'CpipJ2', 'AaegL5', 'BDGP6')
for record in blast_records:
    q = record.query
    bd[q] = defaultdict(list)
    for aln in record.alignments:
        #extract species from hit definition
        (s,c) = [l.split(':')[1:3] for l in aln.hit_def.split(' ') if len(l.split(':')) > 2][0]
        if s not in species:
            raise ValueError('Unknown genome ' + s)
#         if s not in bd[q].keys():
#             bd[q][s] = defaultdict(list)
        for hsp in aln.hsps:
            if (hsp.num_alignments is not None):
                raise ValueError('Many alignmed fragments per HSP:\n' + str(hsp))
            bd[q][s + '_length'].append(hsp.align_length)
            bd[q][s + '_identity'].append(hsp.identities / hsp.align_length)
            bd[q][s + '_q_span'].append('{}-{}'.format(hsp.query_start, hsp.query_end))
            bd[q][s + '_s_span'].append('{}:{}-{}'.format(c, hsp.sbjct_start, hsp.sbjct_end))
            bd[q][s + '_expect'].append(hsp.expect)
            i+=1
        bd[q][s + '_hits'] = len(bd[q][s + '_length'])
    for (k,v) in bd[q].items():
        if type(bd[q][k]) is list:
            bd[q][k] = ';'.join([str(s) for s in v])
            
bd['129']
# example of duplicated marker
# transform keys to int
bd =  {int(k): v for k, v in bd.items()}
bd[1]<jupyter_output><empty_output><jupyter_text>Combine metadata<jupyter_code># load metadata
init_meta = pd.read_csv(META_FILE, index_col=0)
init_meta
# combine metadata
comb_meta = pd.concat([init_meta, pd.DataFrame(bd).T], axis=1)
comb_meta.shape
comb_meta.columns<jupyter_output><empty_output><jupyter_text>Add forgotten data from alignment - unaln_species<jupyter_code># list of all species
all_sp = []
for aln in AlignIO.parse(AMPL_FILE, "maf"):
    if len(aln) == 21:
        all_sp = [seq.id.split('.')[0] for seq in aln]
    break
all_sp
# unaligned species per amplicon
unaln = dict()
i = 0
for aln in AlignIO.parse(AMPL_FILE, "maf"):
    aln_sp = [seq.id.split('.')[0] for seq in aln]
    unaln[i] = ';'.join([sp for sp in all_sp if sp not in aln_sp])
    i+=1
unaln[662]
# add to metadata
comb_meta['unaligned_species'] = pd.Series(unaln)
comb_meta.head()<jupyter_output><empty_output><jupyter_text>Remove markers duplicated in agam<jupyter_code># 72 markers were removed as duplicates in agam
comb_meta = comb_meta[comb_meta['AgamP3_hits'] == 1]
comb_meta.shape
# write filtered alignment 
i = 0
with open(UPD_AMPL_FILE, "w") as o:
    for aln in AlignIO.parse(AMPL_FILE, "maf"):
        if i in comb_meta.index:
            AlignIO.write(aln, o, "maf")
        i += 1
! grep -c 'AgamP3' {UPD_AMPL_FILE}
# reset index for compatibility with filtered alignment
comb_meta.reset_index(drop=True, inplace=True)
comb_meta.index<jupyter_output><empty_output><jupyter_text>Reorder columns and save<jupyter_code># reorder and rename columns
comb_meta = comb_meta[['chr', 'start', 'end', 'seqid', # agam coordinates
            'aligned_species', 'unaligned_species', 'aligned_len', 'indels', 'snvs', # amplicon stats
            'target_start', 'target_end', 'target_aligned_len', 'target_indels', 'target_snvs', # target stats
            'id_lineages', 'informativeness', 'unid_species', # identification power
            'exon', 'gene', 'intron', 'mRNA', 'repeat', 'utr', # agam annotations
            'AgamP3_hits', 'AgamP3_length', 'AgamP3_identity', 
            'AgamP3_q_span', 'AgamP3_s_span', 'AgamP3_expect', # agam blast
            'AaegL5_hits', 'AaegL5_length', 'AaegL5_identity',
            'AaegL5_q_span', 'AaegL5_s_span', 'AaegL5_expect', # aedes blast
            'CpipJ2_hits', 'CpipJ2_length', 'CpipJ2_identity', 
            'CpipJ2_q_span', 'CpipJ2_s_span', 'CpipJ2_expect', # culex blast
            'BDGP6_hits', 'BDGP6_length', 'BDGP6_identity', 
            'BDGP6_q_span', 'BDGP6_s_span','BDGP6_expect' # drosophila blast
            #'id lineage proportion' # duplicate to informativeness, remove
            ]]
comb_meta = comb_meta.rename(columns={"aligned_species": "num_aligned_species", 
                              "indels": "total_indels",
                              "snvs": "total_snvs",
                              "id_lineages":"num_id_lineages",
                              "informativeness": "prop_id_lineages"
                             })
comb_meta.shape
comb_meta.to_excel(UPD_META_FILE)<jupyter_output><empty_output><jupyter_text>Sandbox<jupyter_code># parse hits with SearchIO - it does not extract identities/matches
from Bio import SearchIO
blast_qresult = list(SearchIO.parse(BLAST_FILE, "blast-xml"))
for hit in blast_qresult[129]: 
    for hsp in hit:
        print(hsp)
#hierarchy results->hit->hsp->hsp_framgent<jupyter_output>/Users/am60/miniconda3/envs/biipy/lib/python3.5/site-packages/Bio/SearchIO/__init__.py:211: BiopythonExperimentalWarning: Bio.SearchIO is an experimental submodule which may undergo significant changes prior to its future official release.
  BiopythonExperimentalWarning)
/Users/am60/miniconda3/envs/biipy/lib/python3.5/site-packages/Bio/SearchIO/BlastIO/blast_xml.py:336: BiopythonParserWarning: Renaming hit ID '2R' to a BLAST-generated ID 'gnl|BL_ORD_ID|2318' since the ID was already matched by your query '129'. Your BLAST database may contain duplicate entries.
  (hit.id, hit.blast_id, query_id), BiopythonParserWarning)
/Users/am60/miniconda3/envs/biipy/lib/python3.5/site-packages/Bio/SearchIO/BlastIO/blast_xml.py:336: BiopythonParserWarning: Renaming hit ID '2R' to a BLAST-generated ID 'gnl|BL_ORD_ID|2318' since the ID was already matched by your query '278'. Your BLAST database may contain duplicate entries.
  (hit.id, hit.blast_id, query_id), BiopythonParserWarning)
/Users/am60/mini[...]<jupyter_text>This is an example of 'bad' marker - multiple hits in three of four genomes. Other informative measures - % query coverage and % identity for each hit.Warnings generated during parsing result from same chromosome names in different genomes: here it was 2R in Agam and Dmel. Need to take into account the additional field.<jupyter_code># look at one of the introns
meta = pd.read_csv(META_FILE)
meta[meta['intron'] == 'Yes']
# example intron
for hit in blast_qresult[33]:
    for hsp in hit:
        print(hsp)
# two hits in Aaeg correspond to flanking exons (according to query range), other species don't have that
for hit in blast_qresult[33]:
    #print(hit)
    for hsp in hit:
            print(hsp.midline)
        
print<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12245841035120147, 'fraction_numerical': 0.018946395563770795, 'mean_word_length': 3.8799323562570462, 'pattern_counts': {'":': 5, '<': 36, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '319005', 'license': 'MIT', 'n_tokens_mistral': 3240, 'n_tokens_neox': 3045, 'n_words': 723, 'path': 'work/1_panel_design/20180706_blast_diptera.ipynb', 'repo_name': 'malariagen/ano-spp', 'size': '8656'}