Nextstrain_Peru tutorial - quipupe/Nextstrain_Peru GitHub Wiki

Tutorial

2021, Mayo 19 - 1226 secuencias en GISAID Todos los archivos mencionados se descargan desde: https://github.com/quipupe/Nextstrain_Peru

Descargar secuencias de GISAID

Las secuencias se encuentran aqui: https://github.com/quipupe/Nextstrain_Peru/blob/main/gisaid_auspice_input_hcov-19_2021_05_19_14.tar

image

Se descargan desde GISAID utilizando la opción Augur format. El archivo es un tar que se descomprime de la siguiente manera tar -xvf nombre_del_archivo y luego se utiliza unxz *.xz para descomprimir los archivos generados. Finalmente el archivo de metadata se debe llamar metadata.tsv y el archivo de secuencias sequences.fasta

image

Modificando los metadatos

La tabla debe ser modificada para que esté lista para ser utilizada por Augur.

Primero, revisar la secciòn division para conocer si hay typos en la escritura de las regiones dentro del país.

image

En Perú, no hay muchos problemas, pero sí en el caso de Chile y Argentina. Donde sí se pueden ver estos problemas de entrada de datos en por ejemplo en submitting_lab

image

Prefiero no tener tildes en los nombres así que escribí el siguiente código en R

#Elegir la carpeta a utilizar
setwd('/home/lgm-upch/Desktop/coronavirus/Mayo12/Peru/tutorial/')

#Cargar las bibliotecas necesarias, instalarlas primero.
library(dplyr)
library(tidyverse)

#Leer el archivo de metadatos
input<- read.table("metadata.tsv", header = T, sep = "\t", stringsAsFactors = FALSE)
attach(input)

#Cambiar los nombres, esto puede servir para cambiar cualquier tipo de dato
output <- input
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Instituto Nacional de Salud Perú"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Enteropatógenos. Instituto Nacional de Salud del Perú"] <- "Laboratorio de Referencia Nacional de Enteropatogenos. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Instituto Nacional de Salud Peru"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Genómica Microbiana, Universidad Peruana Cayetano Heredia"] <- "Laboratorio de Genomica Microbiana, Universidad Peruana Cayetano Heredia"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnologia y Biologia Molecular.Instituto Nacional de Salud.Peru"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Instituto Nacional de Salud."] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Instituto Nacional de Salud.Perú"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnologia y Biologia Molecular. Instituto Nacional de Salud Peru"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Instituto Nacional de Salud Perú."] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnologia y Biologia Molecular.Instituto Nacional de Salud.Perú"] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Infecciones Respiratorias Agudas"] <- "Laboratorio de Infecciones Respiratorias Agudas. Centro Nacional de Salud Publica. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Infecciones Respiratorias Agudas. Centro Nacional de Salud Publica,Instituto Nacional de Salud"] <- "LLaboratorio de Infecciones Respiratorias Agudas. Centro Nacional de Salud Publica. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "Laboratorio de Referencia Nacional de Biotecnología y Biología Molecular. Centro Nacional de Salud Publica. Instituto Nacional de Salud Peru."] <- "Laboratorio de Referencia Nacional de Biotecnologia y Biología Molecular. Instituto Nacional de Salud del Peru"
output$submitting_lab[output$submitting_lab== "LLaboratorio de Infecciones Respiratorias Agudas. Centro Nacional de Salud Publica. Instituto Nacional de Salud del Peru"] <- "Laboratorio de Infecciones Respiratorias Agudas. Centro Nacional de Salud Publica. Instituto Nacional de Salud del Peru"

#Cambiar los valores a factores
output$submitting_lab<-factor(output$submitting_lab)

#Obtener una tabla resumen con las frecuencias por cada laboratorio
summary(output$submitting_lab)

#Guardar los metadatos en un nuevo archivo.
write.table(output, file = './metadata2.tsv', quote = F, row.names = F, sep = "\t")

Al nuevo archivo de metadatos se le debe agregar al final estas dos líneas que corresponden a dos muestras de China que se utilizarán como outgroup

hCoV-19/Wuhan/IPBCAMS-WH-01/2019	betacoronavirus	EPI_ISL_402123	?	2019-12-24	Asia	China	Hubei	Wuhan	Asia	China	Hubei	genome	29899	Human	65	Male	?	B	L	Institute of Pathogen Biology, Chinese Academy of Medical Sciences & Peking Union Medical College	Institute of Pathogen Biology, Chinese Academy of Medical Sciences & Peking Union Medical College	Lili Ren, Jianwei Wang, Qi Jin, Zichun Xiang, Zhiqiang Wu, Chao Wu, Yiwei Liu	https://www.gisaid.org/	?	?	2020-01-11	?
hCoV-19/Wuhan/WH01/2019	betacoronavirus	EPI_ISL_406798	?	2019-12-26	Asia	China	Hubei	Wuhan	Asia	China	Hubei	genome	29866	Human	44	Male	?	B	L	General Hospital of Central Theater Command of People's Liberation Army of China	BGI & Institute of Microbiology, Chinese Academy of Sciences & Shandong First Medical University & Shandong Academy of Medical Sciences & General Hospital of Central Theater Command of People's Liberation Army of China	Weijun Chen, Yuhai Bi, Weifeng Shi and Zhenhong Hu	https://www.gisaid.org/	?	?	2020-01-30	?

Estos dos genomas se encuentran en el archivo china.fasta https://github.com/quipupe/Nextstrain_Peru/blob/main/china.fasta

Para incluir las secuencias en el archivo de secuencias final se utiliza cat sequences.fasta china.fasta > sequences2.fasta.

Crear carpetas y adicionar archivos importantes

mkdir auspice
mkdir config
mkdir data
mkdir results

Los archivos metadata2.tsv y sequences2.fasta van a la carpeta data donde se renombran como metadata.fasta y sequences.fasta

La carpeta config debe tener los archivos coronavirus.gb que es el genoma de referencia, lat_longs.tsv con la información de las latitudes y longitudes decimales de cada division o país dependiendo de la escala del análisis y colors.tsv en el cual se puede afinar los colores identificándolos con un código hexagesimal (HEX). Sin embargo, yo lo dejo vacío para que el programa elija los colores automáticamente.

#lat_longs.tsv
division	Amazonas	-5.0	-78.23333
division	Ancash	-9.52779	-77.5277
division	Arequipa	-16.39889	-71.535
division	Apurimac	-14.0	-73.0
division	Ayacucho	-13.15878	-74.22321
division	Cajamarca	-7.15639	-78.5156
division	Callao	-12.0542	-77.1289
division	Cusco	-13.52264	-71.96734
division	Huancavelica	-12.786389	-74.975555
division	Huanuco	-9.93062	-76.24223
division	Ica	-14.0639	-75.7292
division	Junin	-11.15895	-75.99304
division	La Libertad	-9.63306	-77.7417
division	Lambayeque	-6.70111	-79.90611
division	Lima	-12.04318	-77.02824
division	Loreto	-3.75	-73.2444
division	Moquegua	-17.19832	-70.93567
division	Pasco	-10.66748	-76.25668
division	Piura	-5.0	-80.33333
division	San Martin	-7.0	-76.83333
division	Tacna	-18.006569	-70.246277
division	Ucayali	-8.37915	-74.55387

#colors.tsv
division	Amazonas
division	Ancash
division	Apurimac
division	Arequipa
division	Ayacucho 
division	Cajamarca
division	Callao
division	Cusco
division	Huancavelica
division	Huanuco
division	Ica
division	Junin
division	La Libertad
division	Lambayeque
division	Lima
division	Loreto
division	Moquegua
division	Pasco
division	Piura
division	San Martin
division	Tacna
division	Ucayali

Finalmente, el último archivo en colocarse en la carpeta config se llama auspice_config.json que es el archivo de configuración que permite ver los datos de manera didáctica en Nextstrain.

#auspice_config.json
{
  "title": "Genomic epidemiology of SARS-CoV-2 in Peru",
  "maintainers": [
    {"name": "Laboratorio_de_Genomica_Microbiana_Pedro_Romero_pedro.romero@upch.pe", "url": "[email protected]"}
  ],
  "build_url": "none",
  "colorings": [
    {
      "key": "gt",
      "title": "Genotype",
      "type": "categorical"
    },
    {
      "key": "num_date",
      "title": "Date",
      "type": "continuous"
    },
	{
      "key": "submitting_lab",
      "title": "Institution",
      "type": "categorical"
    },
    {
      "key": "sex",
      "title": "Gender",
      "type": "categorical"
    },
    {
      "key": "division",
      "title": "Region",
      "type": "categorical"
    },
    {
      "key": "pangolin_lineage",
      "title": "Lineage",
      "type":"categorical"
    },
    {
      "key":"GISAID_clade",
      "title":"Clade",
      "type":"categorical"
    }

  ],
  "geo_resolutions": [
    "division"
  ],
  "display_defaults": {
    "map_triplicate": true
  },
  "filters": [
    "division",
	"GISAID_clade",
    "sex",
    "pangolin_lineage",
	"institution"
  ],
  "panels": [
	"tree",
	"map",
	"entropy",
	"frequencies"
]
}

Script para correr el programa Auspice y observar los datos con Augur como en Nextstrain

Primero, instalar ambos programas usando la información desde: https://docs.nextstrain.org/en/latest/install-nextstrain.html#install-nextstrain-with-conda-or-docker

Comandos de Augur https://docs.nextstrain.org/projects/augur/en/stable/index.html

Comandos de Auspice https://docs.nextstrain.org/projects/auspice/en/stable/

conda activate nextstrain
sh nextstrain.sh
conda deactivate
conda activate auspice
auspice view --datasetDir auspice/

#En un buscador web tipear
http://localhost:4000/
#y se observarán los resultados

image

#Script nextstrain.sh
#!/bin/bash
#Number of sequences
grep -c ">" data/sequences.fasta;
#Alignment
augur align --sequences data/sequences.fasta --reference-sequence config/coronavirus.gb --output results/aligned.fasta --nthreads 8;

#Tree reconstruction
augur tree --alignment results/aligned.fasta --output results/tree_raw.nwk --nthreads 8;

#Tree calibration
augur refine --tree results/tree_raw.nwk --alignment results/aligned.fasta --metadata data/metadata.tsv --output-tree results/tree.nwk --output-node-data results/branch_lengths.json --timetree --coalescent opt --clock-rate 0.0008 --clock-std-dev 0.0002 --date-confidence --date-inference marginal --clock-filter-iqd 4 --divergence-units mutations --root hCoV-19/Wuhan/WH01/2019  hCoV-19/Wuhan/IPBCAMS-WH-01/2019 --keep-polytomies;

#Frequencies
augur frequencies --metadata data/metadata.tsv --tree results/tree.nwk --method kde -o auspice/Nextstrain_Peru_tip-frequencies.json --proportion-wide 0.1 --pivot-interval 1 --minimal-frequency 0.1 --minimal-clade-size 1;

#Traits
augur traits --tree results/tree.nwk --metadata data/metadata.tsv --output results/traits.json --columns division --confidence;

#Ancestral reconstruction 
augur ancestral --tree results/tree.nwk --alignment results/aligned.fasta --output-node-data results/nt_muts.json --inference joint; 

#Translation
augur translate --tree results/tree.nwk --ancestral-sequences results/nt_muts.json --reference-sequence config/coronavirus.gb --output results/aa_muts.json;

#Export results
augur export v2 --tree results/tree.nwk --metadata data/metadata.tsv --node-data results/branch_lengths.json results/traits.json results/nt_muts.json results/aa_muts.json --colors config/colors.tsv --lat-longs config/lat_longs.tsv --auspice-config config/auspice_config.json --output auspice/Nextstrain_Peru.json --include-root-sequence;

Para compartir los resultados en Nextstrain, crearse una cuenta Github. Dentro de esta, crear el branch master y dentro de este una carpeta auspice donde irán dos archivos, Nextstrain_Peru.json y Nextstrain_Peru_tip-frequencies.json https://github.com/quipupe/Nextstrain_Peru/tree/master/auspice

En mi caso, mi cuenta es quipupe por lo que para observar los archivos online debo entrar a https://nextstrain.org/community/quipupe/Nextstrain_Peru

Ver que los archivos json generados tienen el mismo nombre que la página donde se verán los resultados.

FIN

#Output file from nextstrain.sh
1228

using mafft to align via:
	mafft --reorder --anysymbol --nomemsave --adjustdirection --thread 8 results/aligned.fasta.to_align.fasta 1> results/aligned.fasta 2> results/aligned.fasta.log 

	Katoh et al, Nucleic Acid Research, vol 30, issue 14
	https://doi.org/10.1093%2Fnar%2Fgkf436

109bp insertion at ref position 0
	AAAAAAAAAAAA: hCoV-19/Peru/LIM-INS-102/2020
	ACAATC: hCoV-19/Peru/LAM-HRL-187/2020
	CGCAGAGACAGAAGAAA: hCoV-19/Peru/LAM-HRL-235/2020, hCoV-19/Peru/LAM-HRL-224/2020
	AAAAAAAAAAAAAAAAAAAAAAAAAAAA: hCoV-19/Peru/ICA-INS-109/2020
2bp insertion at ref position 28251
	TG: hCoV-19/Peru/LIM-INS-813/2020, hCoV-19/Peru/LIM-INS-287/2020, hCoV-19/Peru/LIM-INS-288/2020, hCoV-19/Peru/LIM-INS-814/2020, hCoV-19/Peru/LIM-INS-863/2020, hCoV-19/Peru/LIM-INS-897/2020, hCoV-19/Peru/LIM-INS-756/2020, hCoV-19/Peru/CAL-INS-852/2020, hCoV-19/Peru/LIM-INS-886/2020, hCoV-19/Peru/CAL-INS-894/2020, hCoV-19/Peru/CUS-INS-875/2020, hCoV-19/Peru/LIM-INS-858/2020, hCoV-19/Peru/LIM-INS-873/2020, hCoV-19/Peru/LIM-INS-862/2020, hCoV-19/Peru/LIM-INS-834/2020, hCoV-19/Peru/LIM-INS-905/2021, hCoV-19/Peru/LIM-INS-906/2021, hCoV-19/Peru/LIM-INS-878/2020, hCoV-19/Peru/CAL-INS-821/2020, hCoV-19/Peru/CAL-INS-743/2020, hCoV-19/Peru/un-CDC-2-4070016/2021, hCoV-19/Peru/LIM-UPCH-0360/2020, hCoV-19/Peru/LIM-UPCH-0361/2020, hCoV-19/Peru/un-CDC-2-4070031/2021, hCoV-19/Peru/LIM-UPCH-0359/2020, hCoV-19/Peru/LIM-UPCH-0312/2020, hCoV-19/Peru/LIM-UPCH-0321/2020, hCoV-19/Peru/LIM-UPCH-0318/2020, hCoV-19/Peru/un-CDC-2-4070019/2021, hCoV-19/Peru/un-CDC-2-4069942/2021, hCoV-19/Peru/LIM-UPCH-0311/2020, hCoV-19/Peru/HUC-INS-899/2020, hCoV-19/Peru/LIM-UPCH-0357/2020, hCoV-19/Peru/LIM-UPCH-0334/2020
4bp insertion at ref position 28262
	AACA: hCoV-19/Peru/LIM-INS-175/2021, hCoV-19/Peru/un-CDC-2-4069938/2021
76bp insertion at ref position 29789
	GA: hCoV-19/Peru/JUN-INS-127/2020, hCoV-19/Peru/LAM-INS-171/2020, hCoV-19/Peru/LIM-INS-098/2020
Trimmed gaps in MN908947.3 from the alignment
Building a tree via:
	iqtree2 -ninit 2 -n 2 -me 0.05 -nt 8 -s results/aligned-delim.fasta -m GTR  > results/aligned-delim.iqtree.log
	Nguyen et al: IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies.
	Mol. Biol. Evol., 32:268-274. https://doi.org/10.1093/molbev/msu300

augur refine is using TreeTime version 0.8.1

15.45	TreeTime.reroot: with method or node: ['hCoV-19/Wuhan/WH01/2019',
     	'hCoV-19/Wuhan/IPBCAMS-WH-01/2019']

16.21	TreeTime.reroot: with method or node: ['hCoV-19/Wuhan/WH01/2019',
     	'hCoV-19/Wuhan/IPBCAMS-WH-01/2019']
pruning leaf  MN908947.3

16.94	WARNING: Previous versions of TreeTime (<0.7.0) RECONSTRUCTED sequences of
     	tips at positions with AMBIGUOUS bases. This resulted in unexpected
     	behavior is some cases and is no longer done by default. If you want to
     	replace those ambiguous sites with their most likely state, rerun with
     	`reconstruct_tip_states=True` or `--reconstruct-tip-states`.

36.94	TreeTime.reroot: with method or node: ['hCoV-19/Wuhan/WH01/2019',
     	'hCoV-19/Wuhan/IPBCAMS-WH-01/2019']

47.39	###TreeTime.run: INITIAL ROUND

80.37	TreeTime.reroot: with method or node: ['hCoV-19/Wuhan/WH01/2019',
     	'hCoV-19/Wuhan/IPBCAMS-WH-01/2019']

81.26	###TreeTime.run: rerunning timetree after rerooting

117.90	###TreeTime.run: ITERATION 1 out of 2 iterations

174.74	###TreeTime.run: ITERATION 2 out of 2 iterations

380.68	###TreeTime.run: FINAL ROUND - confidence estimation via marginal
      	reconstruction

Inferred a time resolved phylogeny using TreeTime:
	Sagulenko et al. TreeTime: Maximum-likelihood phylodynamic analysis
	Virus Evolution, vol 4, https://academic.oup.com/ve/article/4/1/vex042/4794731

updated tree written to results/tree.nwk
node attributes written to results/branch_lengths.json

augur traits is using TreeTime version 0.8.1
Assigned discrete traits to 1228 out of 1228 taxa.

NOTE: previous versions (<0.7.0) of this command made a 'short-branch
length assumption. TreeTime now optimizes the overall rate numerically
and thus allows for long branches along which multiple changes
accumulated. This is expected to affect estimates of the overall rate
while leaving the relative rates mostly unchanged.

Inferred ancestral states of discrete character using TreeTime:
	Sagulenko et al. TreeTime: Maximum-likelihood phylodynamic analysis
	Virus Evolution, vol 4, https://academic.oup.com/ve/article/4/1/vex042/4794731

results written to results/traits.json

augur ancestral is using TreeTime version 0.8.1

Inferred ancestral sequence states using TreeTime:
	Sagulenko et al. TreeTime: Maximum-likelihood phylodynamic analysis
	Virus Evolution, vol 4, https://academic.oup.com/ve/article/4/1/vex042/4794731

ancestral mutations written to results/nt_muts.json

Read in 12 features from reference sequence file
amino acid mutations written to results/aa_muts.json

Validating schema of 'results/aa_muts.json'...
Validating config file config/auspice_config.json against the JSON schema
Validating schema of 'config/auspice_config.json'...
WARNING: Color map file contains invalid line: division	Amazonas

WARNING: Color map file contains invalid line: division	Ancash

WARNING: Color map file contains invalid line: division	Apurimac

WARNING: Color map file contains invalid line: division	Arequipa

WARNING: Color map file contains invalid line: division	Ayacucho 

WARNING: Color map file contains invalid line: division	Cajamarca

WARNING: Color map file contains invalid line: division	Callao

WARNING: Color map file contains invalid line: division	Cusco

WARNING: Color map file contains invalid line: division	Huancavelica

WARNING: Color map file contains invalid line: division	Huanuco

WARNING: Color map file contains invalid line: division	Ica

WARNING: Color map file contains invalid line: division	Junin

WARNING: Color map file contains invalid line: division	La Libertad

WARNING: Color map file contains invalid line: division	Lambayeque

WARNING: Color map file contains invalid line: division	Lima

WARNING: Color map file contains invalid line: division	Loreto

WARNING: Color map file contains invalid line: division	Moquegua

WARNING: Color map file contains invalid line: division	Pasco

WARNING: Color map file contains invalid line: division	Piura

WARNING: Color map file contains invalid line: division	San Martin

WARNING: Color map file contains invalid line: division	Tacna

WARNING: Color map file contains invalid line: division	Ucayali

WARNING: division->Hubei did not have an associated lat/long value (matching performed in lower case). Auspice won't be able to display this location.

Validating produced JSON
Validating schema of 'auspice/Nextstrain_Peru.json'...
Validating that the JSON is internally consistent...
	WARNING:  "Hubei", a value of the geographic resolution "division", appears in the tree but not in the metadata.
		This will cause transmissions & demes involving this location not to be displayed in Auspice
	WARNING:  The filter "institution" does not appear as a property on any tree nodes.
Validation of 'auspice/Nextstrain_Peru.json' succeeded, but there were warnings you may want to resolve.