7 Mapeo y Binning

M.C. Diana Oaxaca y Dra. Mirna Vazquez Rosas Landa

02 de agosto de 2022

7.1 Secuenciación

El ADN ambiental se secuencia utilizando alguna plataforma de secuenciacion o varias:

7.2 Control de calidad

Del proceso de secuenciacion se recuperan lecturas de aproximadamente 300 bps. Las secuencias las recortamos para eliminar artefactos generados durante el proceso de secuenciación. Aquí algunas herramientas utilizadas:

7.3 Ensamble

Las lecturas se ensamblan para armar contigs y scaffolds. Existen diferentes programas para ensamblar metagenomas, aquí algunos de ellos:

7.3.1 Discutamos

https://docs.google.com/document/d/1iiw-q-90nATg-RNTd9nU8L1XE5xoDRC6j19JC1GiPIk/edit?usp=sharing

7.4 Mapeo

Profundidad: La profundidad de cada contig se calcula mapeando las lecturas al ensamble. Este paso permite evaluar la calidad del ensable y es necesario para hacer la reconstrucción de genomas ya que, como veremos más adelante, es uno de los parámetros que toman en cuenta los “bineadores”.

Vamos a mapear utilizando la herramienta BBMap del programa BBtools. Y samtools.

¡Manos a la obra!

Conéctate al servidor:

ssh USER@hpc-matematicas-z.fciencias.unam.mx

Crea tu carpeta y una liga simbólica a los datos:

mkdir -p 01.Mapeo/{data,results}
cd 01.Mapeo/
ln -s /home/diana/samples/htn/data/htn.fasta data/  
ln -s /home/mirna/03.Mapeo/01.Trimm_reads/htn*.fastq  data/

Primero, activa tu ambiente de conda.

conda activate bbmap_env

¡Ahora sí! Explora las opciones de bbmap, y vamos a hacer nuestro primer mapeo.

bbmap.sh ref=data/htn.fasta in1=data/htn_1-short.fastq in2=data/htn_2-short.fastq out=results/htn.sam kfilter=22 subfilter=15 maxindel=80

cd results
samtools view -bShu htn.sam | samtools sort -@ 94 -o htn_sorted.bam
samtools index htn_sorted.bam

7.4.1 Discutamos

https://docs.google.com/document/d/1iiw-q-90nATg-RNTd9nU8L1XE5xoDRC6j19JC1GiPIk/edit?usp=sharing

7.5 Binning

Utilizaremos varios programas para hacer la reconstrucción de los genomas y haremos una comparación de estos.

NOTA: Cada programa tiene una ayuda y un manual de usuario, es importante revisarlo y conocer cada parámetro que se ejecute. En terminal se puede consultar el manual con el comando man y también se puede consultar la ayuda con -h o --help, por ejemplo fastqc -h.

La presente práctica sólo es una representación del flujo de trabajo, sin embargo, no sustituye los manuales de cada programa y el flujo puede variar dependiendo del tipo de datos y pregunta de investigación.

7.5.1 MaxBin

Crea tu espacio de trabajo y una liga símbólica hacia los datos que se usarán:

mkdir -p 02.MaxBin/{data,results}
cd 02.MaxBin/
ln -s /home/diana/samples/htn/data/htn.fasta data/  
ln -s /home/diana/samples/htn/data/htn-depth.txt  data/

Okay, ahora activa tu ambiente.

conda activate maxbin_env

Explora las opciones y ahora sí, a calcular bins.

run_MaxBin.pl -contig data/htn.fasta -out results/maxbin -abund data/htn-depth.txt -max_iteration 2

7.5.2 MetaBat

Okay vamos a utilizar otro progama. Crea tus ligas simbólicas :)

mkdir -p 03.Metabat/{data,results}
cd 03.Metabat/
ln -s /home/diana/samples/htn/data/htn.fasta data/  
ln -s /home/diana/samples/htn/data/htn_sorted.bam  data/

Para MetaBat lo primero que tenemos que hacer es crear un archivo de profundidad utilizando el script jgi_summarize_bam_contig_depths.

Entonces, primero activamos el ambiente.

conda activate metabat_env

Como cualquier otro programa jgi_summarize_bam_contig_depths tiene opciones, podemos revisarlas.

jgi_summarize_bam_contig_depths --outputDepth data/htn-depth.txt data/htn_sorted.bam

Okay… exploremos el archivo con head

head data/htn-depth.txt

Para metabat sólo necesitamos dos archivos principales:

El ensamble
El archivo de profundidad

metabat -i data/htn.fasta -a data/htn-depth.txt -o results/bins -t 4 --minCVSum 0 --saveCls -d -v --minCV 0.1 -m 2000

7.5.3 CONCOCT

Okay, vamos a utilizar otro programa. Crea tus ligas simbólicas :)

mkdir -p 04.Concoct/{data,results}
cd 04.Concoct/
ln -s /home/diana/samples/htn/data/htn.fasta data/  
ln -s /home/diana/samples/htn/data/htn_sorted.bam  data/

Primero, activemos el ambiente

conda activate concoct_env

Primero, los contigs se tienen que partir en pedazos más pequeños

cut_up_fasta.py data/htn.fasta -c 10000 -o 0 --merge_last -b results/SplitAssembly-htn.bed > results/htn.fasta-split10K.fa

Para crear la tabla de cobertura se necesita primero indexar el archivo bam

samtools index data/htn_sorted.bam

concoct_coverage_table.py results/SplitAssembly-htn.bed data/htn_sorted.bam > results/concoct_coverage_table_htn.tsv

¡Ahora sí! A correr concoct.

Normalmente correríamos 500 iteraciones, pero esta vez sólo haremos una.

concoct --coverage_file results/concoct_coverage_table_htn.tsv --composition_file results/htn.fasta-split10K.fa --clusters 400 --kmer_length 4 --threads 4 --length_threshold 3000 --basename concot --seed 4 --iterations 1

Combinar contigs

merge_cutup_clustering.py concot_clustering_gt3000.csv > results/merged-htn-gt3000.csv

Extraer bins como fasta individualmente

mkdir results/bins-concot
extract_fasta_bins.py  data/htn.fasta results/merged-htn-gt3000.csv --output_path results/bins-concot

7.5.4 Discutamos

https://docs.google.com/document/d/1iiw-q-90nATg-RNTd9nU8L1XE5xoDRC6j19JC1GiPIk/edit?usp=sharing

7.6 Refinamiento

7.6.1 DASTool

Preparing input files.

Fasta_to_Scaffolds2Bin.sh -i /home/mirna/05.Concoct/bins-concot -e fa > htn_concot.scaffolds2bin.tsv

Fasta_to_Scaffolds2Bin.sh -i /home/mirna/04.Metabat2 -e fa > htn.scaffolds2bin.tsv

PATH=/home/programs:$PATH
/home/programs/DAS_Tool-1.1.2/DAS_Tool -i htn_maxbin.contigs2bin.tsv,htn_metabat.scaffolds2bin.tsv,htn_concoct.scaffolds2bin.tsv -l maxbin,metabat,concoct -c data/htn.fasta -o results/htn_bins --debug -t 4  --search_engine diamond --write_bins 1

7.6.2 CheckM

Muy bien, crea un nuevo directorio y entra en él.

mkdir 06.CheckM
cd 06.CheckM/

Ahora activemos el ambiente.

conda activate checkm_env

checkm  lineage_wf -t 4 -x fa /home/mirna/05.DAS_tool/results/htn_bins_DASTool_bins DAStools-log_htn  -f CheckM-DAS_Tool_bins.txt

Vamos a explorar la salida de checkM

Primero me puedes decir ¿Cuántas lineas tiene tu archivo?

Okay… ahora vamos a remover esas lineas feas.

sed -e '1,3d' CheckM-DAS_Tool_bins.txt | sed -e '37d' >CheckM-DAS_Tool_bins_mod.txt

library(tidyverse)
# CheckM -------------------------------------------------------------------####
checkm<-read.table("CheckM-DAS_Tool_bins_mod.txt", sep = "", header = F, na.strings ="", stringsAsFactors= F)
# Extracting good quality bins Megahit ------------------------------------####
colnames(checkm)<-c("Bin_Id", "Marker", "lineage", "Number_of_genomes", 
                         "Number_of_markers", "Number_of_marker_sets", 
                         "0", "1", "2", "3", "4", "5", "Completeness", 
                         "Contamination", "Strain_heterogeneity")  

good_bins<-checkm %>%
  select(Bin_Id, Marker, Completeness, Contamination) %>%
  filter(Completeness >= 50.00) %>%
  filter(Contamination <= 10.00)

Okay… quizá podamos recuperar algunos más.

medium_bins<-checkm %>%
  select(Bin_Id, Marker, Completeness, Contamination) %>%
  filter(Completeness >= 50.00) %>%
  filter(Contamination <= 20.00)

Muy bien, vamos a extraer esos bins.

bins<-medium_bins$Bin_Id

write.table(bins, "lista_medium_bins", quote = F, row.names = F, col.names = F)

7.6.3 Mis Bins

mkdir  -p 08.Bins/{Genoma,Proteoma}
cd 08.Bins

sed 's#bin#cp /home/mirna/05.DAS_tool/results/htn_bins_DASTool_bins/bin#g' lista_medium_bins | sed 's#$#.fa .#g' > copy_bins.sh

Ahora un ejercicio.

grep ">" *.fa

¿Cuál es el problema?

change_bin_name<-function(ruta, ambiente){
ruta_original<-getwd()
setwd(ruta)
filez <- list.files()
newname<-paste0(ambiente, "_", filez)
file.rename(from=filez, to=newname)
filez <- list.files()
file.rename(from=filez, to=sub(pattern="\\.", replacement="_", filez))
setwd(ruta_original)
}

change_bin_name("/home/mirna/07.Bins/Genoma", "htn")

library(phylotools)
library(tidyverse)

add_names_to_seqs <- function(nombre_del_archivo){
  filenames <- unlist(strsplit(nombre_del_archivo, "/"))
  filenames <- filenames[[grep("fa", filenames)]]
  divide <- unlist(strsplit(filenames, "\\."))
  bin_name <- divide[1]
  termination <- divide[2]
  old_name <- get.fasta.name(nombre_del_archivo)
  new_name <- paste0( bin_name, "-scaffold-", old_name) 
  ref2 <- data.frame(old_name, new_name)
  out_file <- paste0(bin_name, "_renamed", ".", termination)
  rename.fasta(infile = nombre_del_archivo, ref_table = ref2, outfile = out_file)
}

files <- list.files(".")
files <- paste0("/home/mirna/07.Bins/Genoma/", files)

map(files, add_names_to_seqs)

Veamos si funcionó

grep ">" *.fa

Muy bien, pongamos eso en una nueva carpeta y esperemos lo mejor jaja. No es cierto, sí movámoslo a otra carpeta, pero quitemos el renamed.

change_bin_name<-function(ruta){
ruta_original<-getwd()
setwd(ruta)
filez <- list.files()
file.rename(from=filez, to=sub(pattern="_renamed", replacement="", filez))
setwd(ruta_original)
}

change_bin_name("/home/mirna/07.Bins/Genoma/01.Bins_named")

¿Creen que puedan optimizar esos scripts? ¡Discute en tu equipo si tienes una mejor idea!