Análisis Genomas: E. coli

Contenido

• Introducción
  • Descarga de archivos
• Anotación de genomas
  • Prokka
  • Bakta
  • AMRfinder
  • Virulence Factor DataBase

Introducción

Vamos a trabajar con ensambles de distintas cepas de E. coli hospedero Bovino.

NOTA ⚠️👁️‍🗨️: Para las figuras del curso se utilizaron todos los 110 genomas con hospedero Bovino. Para este tutorial generaremos un subconjunto de 15, ya incluidos genomas de referencia para patotipos EHEC y EPEC.

Para iniciar vamos a descargar los datos disponibles que existan usando datasets de NCBI, lo puedes instalar aquí.

Con el siguiente código vamos a descargar los genomas depoitados en NCBI de E. coli

mkdir 1_data
cd 1_data

datasets \
summary genome taxon '562' \
--assembly-level complete\
--assembly-source refseq \
--as-json-lines > Escherichia_coli_562.json

dataformat tsv genome \
--inputfile Escherichia_coli_562.json > Escherichia_coli_562.tsv

Este conjunto de metadatos Escherichia_coli_562.tsv nos servirá para seleccionar los ensambles de genoma que necesitamos, este archivo cuenta con 58533 filas y 174 columnas.

Por ello es necesario utilizar la plataforma o método que más prefieras, y quedarnos con un conjunto más pequeño de genomas para descargar.

En este tutorial haremos un ejemplo de parseo de datos con R.

raw_tsv <- readr::read_tsv("0_data/Escherichia_coli_562.tsv")
dim(raw_tsv)


# Change colnames
colnames(raw_tsv) <- gsub(" ", ".", colnames(raw_tsv))

cols <- c("Assembly.BioSample.Accession", "Assembly.BioSample.Attribute.Name",
          "Assembly.BioSample.Attribute.Value", "Assembly.BioSample.Collected.by",
          "Assembly.BioSample.Collection.date", "Assembly.BioSample.Description.Comment",
          "Assembly.BioSample.Description.Organism.Name", "Assembly.BioSample.Description.Organism.Taxonomic.ID",
          "Assembly.BioSample.Description.Title", "Assembly.BioSample.Geographic.location",
          "Assembly.BioSample.Host", "Assembly.BioSample.Host.disease",
          "Assembly.BioSample.Isolation.source", "Assembly.BioSample.Last.updated",
          "Assembly.BioSample.Latitude./.Longitude", "Assembly.BioSample.Serovar",
          "Assembly.BioSample.Owner.Name", "Assembly.BioSample.Strain"
          
)



Refseq_for_taxon<- raw_tsv |>
  dplyr::select('Assembly.Accession',"Assembly.Level",
         'Organism.Name', starts_with('ANI'),starts_with('CheckM'), 
         cols,
         "Assembly.BioProject.Accession") |>
  dplyr::filter(ANI.Best.ANI.match.Organism == "Escherichia coli") |>
  dplyr::filter(grepl("GCF",Assembly.Accession)) |>
  dplyr::distinct()

colnames(Refseq_for_taxon) <- stringr::str_replace(colnames(Refseq_for_taxon), 
                                                   pattern = "Assembly.BioSample.",
                                                   replacement = "")

library(dplyr)


# Create the column 'Assembly.Biosample.Attribute' by combining the values and eliminating duplicates.
ecoli_df <- Refseq_for_taxon |>
  dplyr::group_by(Assembly.Accession) %>%
  summarise(
    Attribute = paste0(
      Attribute.Name, " : ", 
      Attribute.Value, 
      collapse = "; "
    ),
    .groups = "drop"
  ) %>%
  left_join(
    Refseq_for_taxon %>% 
      select(-Attribute.Name, -Attribute.Value) %>% 
      distinct(Assembly.Accession, .keep_all = TRUE),
    by = "Assembly.Accession"
  )


ecoli_subset <- ecoli_df %>% 
  select(where(~ n_distinct(.) > 1)) %>% 
  select(where(~ !all(is.na(.)))) %>% 
  filter(
    !is.na(CheckM.completeness),
    !is.na(Collection.date),
    !is.na(Host)
  )


ecoli_subset <- ecoli_subset %>% 
  mutate(find = CheckM.completeness >= 90 & CheckM.contamination < 5 & ANI.Best.ANI.match.ANI >= 96)


ecoli_subset <- ecoli_subset %>% 
  select(where(~ n_distinct(.) > 1)) %>% 
  select(where(~ !all(is.na(.)))) %>% 
  filter(
    !is.na(CheckM.completeness),
    !is.na(Collection.date),
    !is.na(Host)
  )


ecoli_subset$Host <- tolower(ecoli_subset$Host)

table(ecoli_subset$Host)

ecoli_subset2 <- ecoli_subset[ecoli_subset$Host == "bovine",]

ecoli_subset2$Geographic.location <- tolower(ecoli_subset2$Geographic.location)

ecoli_subset2[stringr::str_detect(ecoli_subset2$Geographic.location, pattern = "china"),
              "Geographic.location"] <- "china"

ecoli_subset2[stringr::str_detect(ecoli_subset2$Geographic.location, pattern = "usa"),
              "Geographic.location"] <- "usa"

ecoli_subset2 <- ecoli_subset2[ecoli_subset2$Geographic.location != "missing",]

# Quiero quedarme al azar con 4 de china y 4 de usa, y conservar a francia y suiza

ecoli_bovine <- ecoli_subset2 %>% 
  group_by(Geographic.location) %>% 
  slice_sample(n = 4) %>% 
  ungroup() %>% 
  select(Assembly.Accession, Strain, Geographic.location,
         Host, Owner.Name, Assembly.BioProject.Accession)

readr::write_tsv(ecoli_bovine, "0_data/Ecoli_bovineHost.tsv")

table(ecoli_subset2$Geographic.location)

En este tutorial este será el archivo con el que trabajaremos Ecoli_bovineHost.tsv:

Assembly.Accession	Strain	Geographic.location	Source	Owner.Name	Assembly.BioProject.Accession	Patotype	Phylogroup
GCF_029625135.1	E47	china	bovine	Chinese Academy of Agricultural Sciences	PRJNA949005	NA	NA
GCF_025908275.1	TL-14	china	bovine	Inner Mongolia Minzu University	PRJNA890057	NA	NA
GCF_029319035.1	E94	china	bovine	Chinese Academy of Agricultural Sciences	PRJNA942474	NA	NA
GCF_029227895.1	ZYB62	china	bovine	Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Science	PRJNA727665	NA	NA
GCF_009931435.1	EC42405	france	bovine	Anses	PRJNA556131	NA	NA
GCF_004771235.1	PF9285	switzerland	bovine	University of Bern	PRJNA530748	NA	NA
GCF_049267065.1	F3_13_1_3	usa	bovine	United State Department of Agriculture	PRJNA728709	NA	NA
GCF_013168215.1	7409	usa	bovine	United State Department of Agriculture	PRJNA528413	NA	NA
GCF_013167395.1	NE122	usa	bovine	United State Department of Agriculture	PRJNA528413	NA	NA
GCF_048571305.1	7_12_30A	usa	bovine	United State Department of Agriculture	PRJNA728709	NA	NA
GCF_000025165.1	CB9615	Germany	reference	Nankai University	PRJNA224116	EPEC	E
GCF_000026545.1	E2348/69	UK	reference	The Wellcome Trust Sanger Institute	PRJNA224116	EPEC	B2
GCF_000026265.1	IAI1	-	reference	EBI	PRJNA33373	Comensal	B1
GCF_000005845.2	K-12 MG1655	usa	reference	Univ. Wisconsin	PRJNA225	Comensal	A
GCF_000010745.1	12009	japan	reference	University of Tokyo	PRJDA32511	EHEC	B1
GCF_029962285.1	EDL933	china	reference	East China University of Science and Technology	PRJNA963085	EHEC	E

Descarga de archivos

Para descargar los archivos vamos a utilizar el siguiente código:

for g in $(cut -f1 Ecoli_bovineHost.tsv | awk 'NR > 1'); do
        if [ -f $g.zip ]; then
            echo $g'.zip exists';
        else
            datasets download genome accession $g \ 
            --include genome --filename $g'.zip';
            sleep 2;
        fi
    done

Descomprimir los archivos

for g in $(cut -f1 Ecoli_bovineHost.tsv | awk 'NR > 1'); do
        unzip $g.zip -d $g;
        mv $g/ncbi_dataset/data/$g/*.fna $g.fna;
        rm -r $g/;
done

rm *.zip

Con esto ya contamos con un sub conjunto de datos listo para trabajar.

Anotación de genomas

Prokka

Prokka es un anotador de genomas para organismos procariontes, lo puedes descargar de aquí.

genomes=$( ls $HOME/1_data/*.fna )

for g in $genomes; do
 prokka \
 $g \
 --cpus 6 \
 --outdir prokka/${id%.fna} \
 --genus Escherichia
done
 

Bakta

Bakta es uno de los anotadores de genomas “más” completos ya que usa distintas bases, utiliza a DIAMOND de alineador. Lo puedes descargar aquí. Es necesario descargar la base de datos para poder ejecutarlo, para este ejemplo se guardó en el directorio $HOME/db

Precaución ⚠️: Necesitas 8GB de RAM, disponibilidad de tiempo y memoria en tu máquina, con más de 10 genomas podría demorar según la capacidad de nucleos que tengas. Para este tutorial solo te mostraremos como es el código pero no se utilizará este paso para los siguientes análisis.

source $HOME/bin/bakta-env/bin/activate

genomes=$(ls $HOME/1_data)

experiment_file="$HOME/Experiment_Design.tsv"


for g in $genomes; do
    id=$(basename "$g" | awk -F'/' '{print $NF}' | sed 's/.fna$//')

    strain=$(awk -F'\t' -v id="$id" '$1 == id {print $3}' "$experiment_file")

    bakta \
    --db $HOME/db \
    --output $HOME/2_annotation/bakta/${id} \
    --genus Escherichia \
    --strain "$strain" \
    --threads 6 \
    --force $HOME/1_data/${g}
done

AMRfinder

workdir=$HOME/2_annotation/prokka
output_dir=$HOME/2_annotation/AMR

for dir in $workdir/*; do
  if [ -d "$dir" ]; then
      sample=$(basename "$dir")
      protein=$dir/$sample.faa
      dna=$dir/$sample.fna
      gff=$dir/$sample.gff3

  temp_gff="${output_dir}/amrfinder_${sample}.gff"
        perl -pe '/^##FASTA/ && exit; s/(\W)Name=/$1OldName=/i; s/ID=([^;]+)/ID=$1;Name=$1/' $gff > $temp_gff

      amrfinder \
      -p $protein \
      -n $dna \
      -g $temp_gff \
      --plus --organism Escherichia \
      -o $output_dir/${sample}_amrfinder.tsv \
      --threads 6
done

Virulence Factor DataBase