<jupyter_start><jupyter_text>Alignment and tree building using Nextstrain[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/imperial-genomics-facility/viral-genome-notebook-image/master?urlpath=lab%2Ftree%2Fexamples%2Fcoronavirus_analysis_alignment_and_tree_building_using_Nextstrain.ipynb)* __Notebook version__: `v0.0.1`* __Created by:__ [Dr. <NAME>](https://www.linkedin.com/in/hiren-ghosh-phd-6181853a), `Imperial BRC Genomics Facility`* __Maintained by:__ `Imperial BRC Genomics Facility`* __Docker image:__ [imperialgenomicsfacility/viral-genome-analysis-notebooks](https://hub.docker.com/r/imperialgenomicsfacility/viral-genome-analysis-notebooks)* __Github repository:__ [imperial-genomics-facility/viral-genome-notebook-image](https://github.com/imperial-genomics-facility/viral-genome-notebook-image)* __Created on:__ `2020-April-21 14:42`* __Contact us:__ [Imperial BRC Genomics Facility](https://www.imperial.ac.uk/medicine/research-and-impact/facilities/genomics-facility/contact/)* __License:__ [Apache License 2.0](https://github.com/imperial-genomics-facility/scanpy-notebook-image/blob/master/LICENSE) Introduction[Nextstrain](https://nextstrain.org/docs/getting-started/introduction) is an open-source project to harness the scientific and public health potential of pathogen genome data. Real time tracking of pathogen evolution data and visualization is provided by [nextstrain.org](https://nextstrain.org/). For example, the real time tracking of genomic data for COVID-19 pandemic can be found here: [nextstrain.org/ncov/global](https://nextstrain.org/ncov/global). This notebook will focus on using the bioinformatics toolkits provided by [Nextstrain](https://nextstrain.org/docs/getting-started/introduction) and setting up a docker based local [Nextstrain](https://nextstrain.org/docs/getting-started/introduction) server for output visualization.  Configure notebook for run<jupyter_code>CPU_THREADS = 1<jupyter_output><empty_output><jupyter_text>Load required python libraries<jupyter_code>import os,requests,json,time
import pandas as pd<jupyter_output><empty_output><jupyter_text>Fetch genome sequences<jupyter_code>def fetch_genome_fasta_from_ncbi(refseq_id,output_path='.',file_format='fasta'):
  '''
  A function for fetching the genome fasta sequences from NCBI
  
  :param refseq_id: NCBI genome id
  :param output_path: Path to dump genome files, default '.'
  :param file_format: Output fileformat, default fasta, supported formats are 'fasta' and 'gb'
  :return: output_file
  '''
  try:
    url = \
      'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id={0}&rettype={1}'.\
        format(refseq_id,file_format)
    r = requests.get(url)
    if r.status_code != 200:
      raise ValueError(
              'Failed to download data from NCBI for id {0}, error code: {1}'.\
                format(refseq_id,r.status_code))
    fasta_data = r.content.decode('utf-8')
    output_file = \
      os.path.join(
        os.path.abspath(output_path),
        '{0}.{1}'.format(refseq_id,file_format))
    with open(output_file,'w') as fp:
      fp.write(fasta_data)
    print('Downloaded genome seq for {0}'.format(refseq_id))
    return output_file
  except Exception as e:
    print('Failed to download data for {0} from NCBI, error: {1}'.format(refseq_id,e))
genome_list = [
    'MN988668.1',
    'MN938384.1',
    'MN975262.1',
    'MN985325.1',
    'MN988713.1',
    'MN994467.1',
    'MN994468.1',
    'MN997409.1'
]
fasta_list = []
for id in genome_list:
  fasta_list.\
    append(
      fetch_genome_fasta_from_ncbi(
        refseq_id=id,
        output_path='/tmp'))
  time.sleep(2)
fasta_list
reference_gb = \
  fetch_genome_fasta_from_ncbi(
    'NC_045512.2',
    file_format='gb',
    output_path='/tmp')
reference_gb
MERGED_FASTA = '/tmp/corona_merged.fasta'
!rm -f $MERGED_FASTA
for f in fasta_list:
  !cat $f >> $MERGED_FASTA

## count the genome sequences in the merged fasta file
!grep '>' $MERGED_FASTA|wc -l<jupyter_output>8<jupyter_text>Multiple sequence alignment and tree building<jupyter_code>ALIGNED_FASTA = '/tmp/corona_merged.afa'
!augur align \
  --sequences $MERGED_FASTA \
  --reference-sequence $reference_gb \
  --nthreads $CPU_THREADS \
  --output $ALIGNED_FASTA \
  --fill-gaps
RAW_TREE = '/tmp/corona_raw.tree'
!augur tree \
  --alignment $ALIGNED_FASTA \
  --method fasttree \
  --output $RAW_TREE \
  --nthreads $CPU_THREADS<jupyter_output>Cannot specify model unless using IQTree. Model specification ignored.
Building a tree via:
	FastTreeMP -nosupport -nt /tmp/corona_merged.afa  1> /tmp/corona_raw.tree 2> /tmp/corona_raw.tree.log
	Price et al: FastTree 2 - Approximately Maximum-Likelihood Trees for Large Alignments.
	PLoS ONE 5(3): e9490. https://doi.org/10.1371/journal.pone.0009490


Building original tree took 1.8078341484069824 seconds<jupyter_text>Metadata configs for Nextstrain<jupyter_code>metadata = [{
  'strain':'NC_045512.2',
  'virus':'SARS-CoV-2',
  'accession':'NC_045512',
  'date':'2020-03-30',
  'region':'China',
  'country':'China',
  'division':'China',
  'city':'Wuhan',
  'db':'genbank',
  'segment':'genome',
  'authors':'Baranov,P.V., Henderson,C.M., Anderson,C.B., Gesteland,R.F.,Atkins,J.F. and Howard,M.T.',
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2',
  'title':'Programmed ribosomal frameshifting in decoding the SARS-CoV genome',
  'journal':'Virology 332 (2), 498-510 (2005)',
  'paper_url':'https://www.ncbi.nlm.nih.gov/pubmed/15680415'
  },
  {
  'strain':'MN988668.1',
  'virus':'SARS-CoV-2',
  'accession':'MN988668',
  'date':'2020-03-11',
  'region':'China',
  'country':'China',
  'division':'China',
  'city':'Wuhan',
  'db':'genbank',
  'segment':'genome',
  'authors':'Chen,L., Liu,W., Zhang,Q., Xu,K., Ye,G., Wu,W., Sun,Z., Liu,F.,Wu,K., Mei,Y., Zhang,W., Chen,Y., Li,Y., Shi,M., Lan,K. and Liu,Y.',
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN988668.1',
  'title':'RNA based mNGS approach identifies a novel human coronavirus from two individual pneumonia cases in 2019 Wuhan outbreak',
  'journal':'Emerg Microbes Infect (2019) In press',
  'paper_url':''
  },
  {
  'strain':'MN938384.1',
  'virus':'SARS-CoV-2',
  'accession':'MN938384',
  'date':'2020-02-11',
  'region':'China',
  'country':'China',
  'division':'China',
  'city':'Shenzhen',
  'db':'genbank',
  'segment':'genome',
  'authors':"""Chan,J.F.-W., Yuan,S., Kok,K.H., To,K.K.-W., Chu,H., Yang,J.,
            Xing,F., <NAME>., Poon,R.W.-S., Tsai,H.W., Lo,S.K.-F.,
            Chan,K.H., Poon,V.K.-M., Chan,W.M., Ip,J.D., Cai,J.P.,
            Cheng,V.C.-C., Chen,H., Hui,C.K.-M. and Yuen,K.Y.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN938384.1',
  'title':'A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster',
  'journal':'Lancet (2020) In press',
  'paper_url':''
  },
  {
  'strain':'MN975262.1',
  'virus':'SARS-CoV-2',
  'accession':'MN975262',
  'date':'2020-02-11',
  'region':'China',
  'country':'China',
  'division':'China',
  'city':'Shenzhen',
  'db':'genbank',
  'segment':'genome',
  'authors':"""Chan,J.F.-W., Yuan,S., <NAME>,K.K.-W., Chu,H., Yang,J.,
            Xing,F., Liu,J., Yip,C.C.-Y., Poon,R.W.-S., Tsai,H.W., Lo,S.K.-F.,
            Chan,K.H., Poon,V.K.-M., Chan,W.M., Ip,J.D., Cai,J.P.,
            Cheng,V.C.-C., Chen,H., Hui,C.K.-M. and Yuen,K.Y.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN975262.1',
  'title':"""A familial cluster of pneumonia associated with the 2019 novel
            coronavirus indicating person-to-person transmission: a study of a
            family cluster""",
  'journal':'Lancet (2020) In press',
  'paper_url':''
  },
  {
  'strain':'MN985325.1',
  'virus':'SARS-CoV-2',
  'accession':'MN985325',
  'date':'2020-03-27',
  'region':'USA',
  'country':'USA',
  'division':'USA',
  'city':'Atlanta',
  'db':'genbank',
  'segment':'genome',
  'authors':"""<NAME>., Lu,X., Kamili,S., Sakthivel,S.K., Murray,J.,
            Queen,K., Tao,Y., Paden,C.R., Zhang,J., Li,Y., Uehara,A., Wang,H.,
            Goldsmith,C., Bullock,H.A., Wang,L., Whitaker,B., Lynch,B.,
            Gautam,R., Schindewolf,C., Lokugamage,K.G., Scharton,D.,
            Plante,J.A., Mirchandani,D., Widen,S.G., Narayanan,K., Makino,S.,
            Ksiazek,T.G., Plante,K.S., Weaver,S.C., <NAME>.,
            Menachery,V.D. and Thornburg,N.J.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN985325.1',
  'title':"""Severe Acute Respiratory Syndrome Coronavirus 2 from Patient with
            2019 Novel Coronavirus Disease, United States""",
  'journal':'Emerging Infect. Dis. 26 (6) (2020) In press',
  'paper_url':''
  },
  {
  'strain':'MN988713.1',
  'virus':'SARS-CoV-2',
  'accession':'MN988713',
  'date':'2020-02-11',
  'region':'USA',
  'country':'USA',
  'division':'USA',
  'city':'Illinois',
  'db':'genbank',
  'segment':'genome',
  'authors':"""Tao,Y., Queen,K., Paden,C.R., Zhang,J., Li,Y., Uehara,A., Lu,X.,
            Lynch,B., Sakthivel,S.K.K., Whitaker,B.L., Kamili,S., Wang,L.,
            Murray,J.R., Gerber,S.I., Lindstrom,S. and Tong,S.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN988713.1',
  'title':'Novel betacoronavirus, Illinois',
  'journal':'',
  'paper_url':''
  },
  {
  'strain':'MN994467.1',
  'virus':'SARS-CoV-2',
  'accession':'MN994467',
  'date':'2020-02-11',
  'region':'USA',
  'country':'USA',
  'division':'USA',
  'city':'California',
  'db':'genbank',
  'segment':'genome',
  'authors':"""<NAME>., Li,Y., Paden,C.R., Zhang,J., Lu,X.,
            Lynch,B., Sakthivel,S.K.K., Whitaker,B.L., Kamili,S., Wang,L.,
            Murray,J.R., Gerber,S.I., Lindstrom,S. and Tong,S.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN994467.1',
  'title':'nCoV-2019 California case 1',
  'journal':'',
  'paper_url':''
  },
  {
  'strain':'MN994468.1',
  'virus':'SARS-CoV-2',
  'accession':'MN994468',
  'date':'2020-02-11',
  'region':'USA',
  'country':'USA',
  'division':'USA',
  'city':'California',
  'db':'genbank',
  'segment':'genome',
  'authors':"""Uehara,A., Queen,K., Tao,Y., Li,Y., Paden,C.R., Zhang,J., Lu,X.,
            Lynch,B., Sakthivel,S.K.K., Whitaker,B.L., Kamili,S., Wang,L.,
            Murray,J.R., Gerber,S.I., Lindstrom,S. and Tong,S.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN994468.1',
  'title':'nCoV-2019 California case 2',
  'journal':'',
  'paper_url':''
  },
  {
  'strain':'MN997409.1',
  'virus':'SARS-CoV-2',
  'accession':'MN997409',
  'date':'2020-02-11',
  'region':'USA',
  'country':'USA',
  'division':'USA',
  'city':'Arizona',
  'db':'genbank',
  'segment':'genome',
  'authors':"""Tao,Y., Paden,C.R., Queen,K., Uehara,A., Li,Y., Zhang,J., Lu,X.,
            Lynch,B., Sakthivel,S.K.K., Whitaker,B.L., Kamili,S., Wang,L.,
            Murray,J.R., Gerber,S.I., Lindstrom,S. and Tong,S.""",
  'url':'https://www.ncbi.nlm.nih.gov/nuccore/MN997409.1',
  'title':'nCoV-2019 sequence from Arizona case',
  'journal':'',
  'paper_url':''
  }
]

## load metadata and write tsv files
METADATA_FILE = '/tmp/metadata.tsv'

df = pd.DataFrame(metadata)
df.to_csv(METADATA_FILE,sep='\t',index=False)
REFINED_TREE = '/tmp/corona_refined.tree'
BRANCH_LENGTHS_JSON = '/tmp/corona_branch_lengths.json'
!augur refine \
  --tree $RAW_TREE \
  --alignment $ALIGNED_FASTA \
  --metadata $METADATA_FILE \
  --output-tree $REFINED_TREE \
  --output-node-data $BRANCH_LENGTHS_JSON \
  --timetree \
  --coalescent opt \
  --date-confidence \
  --date-inference marginal \
  --clock-filter-iqd 4
## Annotate the Phylogeny
### Reconstruct Ancestral Traits
TRAITS_JSON = '/tmp/corona_traits.json'
!augur traits \
  --tree $REFINED_TREE \
  --metadata $METADATA_FILE \
  --output $TRAITS_JSON \
  --columns region country \
  --confidence
### Infer Ancestral Sequences
NT_MUTS_JSON = '/tmp/corona_nt_muts.json'
!augur ancestral \
  --tree $REFINED_TREE \
  --alignment $ALIGNED_FASTA \
  --output-node-data $NT_MUTS_JSON \
  --inference joint
### Identify Amino-Acid Mutations
AA_MUTS_JSON = '/tmp/corona_aa_muts.json'
!augur translate \
  --tree $REFINED_TREE \
  --ancestral-sequences $NT_MUTS_JSON \
  --reference-sequence $reference_gb \
  --output $AA_MUTS_JSON
%%file /tmp/lat_longs.tsv
country,china,31.0740838,107.7450643
country,usa,42.6584011,-80.2716978
region,china,31.0740838,107.7450643
region,usa,42.6584011,-80.2716978
%%file /tmp/colors.tsv
country,china,#511EA8
country,usa,#cd5700
region,china,#511EA8
region,usa,#cd5700
!sed -i 's|,|\t|g' /tmp/lat_longs.tsv
!sed -i 's|,|\t|g' /tmp/colors.tsv
%%file /tmp/auspice_config.json
{
  "title": "Covid-19 5 samples",
  "maintainers": [
    {"name": "<NAME>", "url": "your url"}
  ],
  "build_url": "your build url",
  "colorings": [
    {
      "key": "gt",
      "title": "Genotype",
      "type": "categorical"
    },
    {
      "key": "num_date",
      "title": "Date",
      "type": "continuous"
    },
    {
      "key": "author",
      "title": "Author",
      "type": "categorical"
    },
    {
      "key": "country",
      "title": "Country",
      "type": "categorical"
    },
    {
      "key": "region",
      "title": "Region",
      "type": "categorical"
    }
  ],
  "geo_resolutions": [
    "country",
    "region"
  ],
  "panels": [
     "tree",
     "map"
  ],
  "display_defaults": {
    "map_triplicate": true
  },
  "filters": [
    "country",
    "region",
    "author"
  ]
}
### Export the Results
NEXTSTRAIN_SERVER_JSON = 'corona.json'
!augur export v2 \
  --tree $REFINED_TREE \
  --metadata $METADATA_FILE \
  --node-data $BRANCH_LENGTHS_JSON \
              $TRAITS_JSON \
              $NT_MUTS_JSON \
              $AA_MUTS_JSON \
  --colors /tmp/colors.tsv \
  --lat-longs /tmp/lat_longs.tsv \
  --auspice-config /tmp/auspice_config.json \
  --output $NEXTSTRAIN_SERVER_JSON<jupyter_output>Validating schema of '/tmp/corona_aa_muts.json'...
Validating config file /tmp/auspice_config.json against the JSON schema
Validating schema of '/tmp/auspice_config.json'...
Validating produced JSON
Validating schema of 'corona.json'...
Validating that the JSON is internally consistent...
Validation of 'corona.json' succeeded.
{'fraction_non_alphanumeric': 0.19017063883796362, 'fraction_numerical': 0.040473840078973346, 'mean_word_length': 3.647116644823067, 'pattern_counts': {'":': 26, '<': 25, '<?xml version=': 0, '>': 30, 'https://': 24, 'lorem ipsum': 0, 'www.': 12, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '446759', 'license': 'Apache-2.0', 'n_tokens_mistral': 6215, 'n_tokens_neox': 5446, 'n_words': 970, 'path': 'examples/coronavirus_analysis_alignment_and_tree_building_using_Nextstrain.ipynb', 'repo_name': 'imperial-genomics-facility/viral-genome-notebook-image', 'size': '14256'}