README.md

# Workflow Metagenomics
# Tim Reska 





![Image](https://ascgitlab.helmholtz-muenchen.de/ttreska57/PhD_Notes/-/blob/main/Copy%20of%20computational_workflow_overviewv0.2.png)

# Applied software and the corresponding console commands
## basecalling 
-guppy basecaller

The first step for each sequencing run is to basecall your data, i.e., produce fastq sequencing files from the initial fast5 files. We will use the standard basecaller is guppy developed by Oxford Nanopore.

```shell
$ guppy_basecaller -i /input_directory_fast5 -s /guppy_output -c <configuration_filge.cfg> [options]
```
The command above will call guppy on the input fast5 directory option (-i), write the output to the directory given with option -s in fastq format (option -o). The configuration file comprises the used nanopore device and the corresponding kit


## Data Analysis and Quality Control
-Nanoplot

After the basecalling step we want an overview of the data: what is the number of reads, average length, quality etc.

```shell
$ NanoPlot --fastq /guppy_output -o nanoplot/ --title Passed_reads --loglength
```

The options --fastq and -o set the input and output for the command. We also want the length of the sequences shown log-transformed (--loglength) and also set a title for the plots (--title). The above command will create several plots and summary files including a *.html file that can be opened in your browser.

-fastqc 


## Adapter Removal
-Porechop

Porechop is a tool for finding and removing adapters from Oxford Nanopore reads. Adapters on the ends of reads are trimmed off, and when a read has an adapter in its middle, it is treated as chimeric and chopped into separate reads

```shell
$ porechop -i <fastq_output_guppy.fastq> -o <porechopped.fastq> 
```
The above command will use the default values of porechop to search for adapters in the input fastq file, trim the reads and write them to file porechopped.fastq

## Filtering
-NanoFilt

NanoFilt is a python script that can be used to filter reads based on length and quality as well as trimming parts of the sequence

```shell
$ NanoFilt -l 200 -q 7 < <porechopped.fastq> > <nanofiltered.fastq>
```
NanoFilt does not provide options for input or output files. Therefore we will use the two operators > and <

With the option -l 200 we remove all sequences shorter than 200 nucleotides and with -q we filter all reads with a phred quality score lower than 7

## Assembly
-Flye

Flye is a de novo assembler for long-reads. Flye also produces a polished consensus sequence for the assembly which significantly reduces the error rate.

```shell
$ flye --meta --nano-hq <nanofiltered.fastq> -o /output_dir --genome-size 1m -i 1
```
The above command uses the high-quality polished nanopore reads and assembles the genome. The estimated --genome-size is an obligatory input. The --meta option is for metagenomic approaches. The -i option is for polishing iterations.

## Polishing
-minimap2

Minimap2 is a versatile sequence alignment program that aligns DNA against a reference or draft genome. 

```shell
$ minimap2 -ax map-ont -t 8 <flye_assembly.fa> <nanofiltered.fastq> > <aligned_reads.sam>
```

The above command uses the raw nanopore reads and maps them to the draft genome generated by flye. The option -ax map-ont sets the default parameters for nanopore reads.

Minimap2 is also used to find overlaps between sequences.

```shell
$ minimap2 -ax map-ont -t 8 <raw_ont_reads.fa> <raw_ont_reads.fasta> > <overlaps.paf>
```
-Racon

Consensus assemblies try to reduce error rates by choosing the most likely sequence of a given assembly and a set of raw reads. Although this does not incorporate the raw signal information of the flow cell to correct individual reads it can significantly improve the quality of an assembly. To improve a draft assembly with racon map the reads that should be used for error correction against the assembly. Use minimap to map the trimmed reads from prac3/nanofilt against the miniasm assembly and subsequently use the filtered reads and the mapping to build the consensus assembly

```shell
$ racon racon [options ...] <correctedreads.fasta> <overlaps.sam|paf> <flyeassembly.fasta>
```
-Medaka

Medaka is a tool to create consensus sequences and variant calls from nanopore sequencing data. This task is performed using neural networks applied a pileup of individual sequencing reads against a draft assembly

```shell
$ medaka_consensus –d flye_assembly.fasta -i filtered.fastq \ -o medaka_output –m R941_min_high_g303 –b 50
```

The above command runs medaka on the flye assembly using the filtered reads. The –m model specifies the basecalling model used by guppy to produce the reads. Medaka will write its results to the folder medaka_output.

## Sorting
-Samtools

```shell
$ samtools view -bS <aligned_reads.sam> > <aligned_reads.bam> [options]
$ samtools sort <aligned_sample.bam> -o <sample_reads_sorted.bam> [options]
```

## Binning
-MetaBat2

MetaBAT2 takes a metagenome assembly and the reads that produced the assembly and organizes the contigs into putative genomes, called "bins".

```shell
$ runMetaBat.sh <options> <flye_assembly.fasta> <sample_reads_sorted.bam>
```
The above command runs medaka and the following files are the Output Object:The BinnedContig Object represents the directory of binned contigs created by MetaBAT2. This object can be used for downstream analysis. Output Bin Summary Report:The number of bins produced, the number of contigs that were binned and the total number of contigs in the assembly. Downloadable files: The enitre output of the MetaBAT2 run may be downloaded as a zip file. This zip file also contains a table of read-depth coverage per contig ("*.depth.txt")

-MaxBin2

Let’s run a MaxBin binning on the flye assembly. First, we need to generate an abundance file from the mappes reads:

```shell
pileup.sh in=megahit.sam out=cov.txt | awk '{print $1"\t"$5}' cov.txt | grep -v '^#' > abundance.txt
qsub -cwd -pe multislot 24 -N maxbin -l mtc=1 -b y \ run_MaxBin.pl -thread 24 -contig final.contigs.fa -out maxbin -abund abundance.tx
```

MaxBin2 clusters metagenomic contigs (assembled contiguous genome fragments) into different "bins", each of which corresponds to a putative genome. It uses nucleotide composition information, source strain abundance (measured by depth-of-coverage by aligning the reads to the contigs), and phylogenetic marker genes to perform binning through an Expectation-Maximization (EM) algorithm.

-Vamb
First we have to combine all our contigs to one fasta-list

```shell
concatenate.py /path/to/catalogue.fna.gz /path/to/assemblies/sample1/contigs.fasta
/path/to/assemblies/sample2/contigs.fasta
```
Second we use our aligned and sorted reads together with the fasta catalogue to generate the corresponding bins.
```shell
vamb --outdir path/to/outdir --fasta /path/to/catalogue.fna.gz --bamfiles /path/to/bam/*.bam -o C --minfasta 200000
```

Vamb is a metagenomic binner which feeds sequence composition information from a contig catalogue and co-abundance information from BAM files into a variational autoencoder and clusters the latent representation. It performs excellently with multiple samples, and pretty good on single-sample data. Vamb is implemented purely in Python (with a little bit of Cython) and can be used both from command line and from within a Python interpreter.

## Bin Refining
-DAS Tool

General command: 

```shell
DAS_Tool [options] -i <contig2bin> -c <contigs_fasta> -o <outputbasename>
```

Example command:

```shell
DAS_Tool  -i sample_data/sample.vamb_contigs2bin.tsv,\
sample_data/sample.maxbin2_contigs2bin.tsv,\
sample_data/sample.metabat_contigs2bin.tsv,\
-l vamb,maxbin,metabat \
-c sample_data/sample.contigs.fa \
-o sample_output/DASToolRun
```


DAS Tool is an automated method that integrates the results of a flexible number of binning algorithms to calculate an optimized, non-redundant set of bins from a single assembly


## Bin coverage
-CoverM

CoverM calculates the coverage of a set of reads on a set of contigs

```shell
coverm contig --coupled read1.fastq.gz read2.fastq.gz --reference assembly.fna
```

CoverM calculates coverage of genomes/MAGs coverm genome (help) or individual contigs coverm contig (help). Calculating coverage by read mapping, its input can either be BAM files sorted by reference, or raw reads and reference genomes in various formats

## Bin QC
-CheckM

CheckM works on a directory of genome bins in FASTA format. By default, CheckM assumes genomes consist of contigs/scaffolds in nucleotide space and that the files to process end with the extension fna. You can specify a different extension with the –x flag. CheckM calls genes internally using prodigal, taking care to identify genes with recoded stop codons

```shell
checkm lineage_wf -t 8 -x fa /home/bins /home/checkm
```
Or, to process files of called genes in amino acid space which have the extension faa, use:

```shell
checkm lineage_wf --genes -t 8 -x faa <bin folder> <output folder>
```

CheckM provides a set of tools for assessing the quality of genomes recovered from isolates, single cells, or metagenomes. It provides robust estimates of genome completeness and contamination by using collocated sets of genes that are ubiquitous and single-copy within a phylogenetic lineage. CheckM also provides tools for identifying genome bins that are likely candidates for merging based on marker set compatibility, similarity in genomic characteristics, and proximity within a reference genome tree.

## Bin Classification
-GDTB-Tk

GTDB-Tk is a software toolkit for assigning objective taxonomic classifications to bacterial and archaeal genomes based on the Genome Database Taxonomy 


Step 1: Gene calling (identify)
```shell
gtdbtk identify --genome_dir /gtdbtk/genomes --out_dir /gtdbtk/identify --extension gz 
```

Step 2: Aligning genomes: The align step will align all identified markers, determine the most likely domain and output a concatenated MSA
```shell
gtdbtk align --identify_dir /gtdbtk/identify --out_dir /gtdbtk/align
``` 

Step 3: Classifying genomes (classify)
The classify step will place the genomes into a reference tree, then determine their most likely classification.
```shell
!gtdbtk classify --genome_dir /gtdbtk/genomes --align_dir /gtdbtk/align --out_dir /gtdbtk/classify -x gz 
```

## Bin Clustering
-dRep

dRep is a python program which performs rapid pair-wise comparison of genome sets. One of it’s major purposes is for genome de-replication, but it can do a lot more.

It has in general three commands: compare, dereplicate, and check dependencies. For creating a set of non redundant genomes we use dRep dereplicate.

```shell
dRep dereplicate complete_only -g *.fna --S_algorithm gANI
```

De-replication is the process of identifying sets of genomes that are the “same” in a list of genomes, and removing all but the “best” genome from each redundant set. How similar genomes need to be to be considered “same”, how to determine which genome is “best”, and other important decisions are discussed in Important Concepts

## Consensus error calculation
-QUAST

QUAST evaluates genome assemblies. QUAST works both with and without a reference genome. The tool accepts multiple assemblies, thus is suitable for comparison. 

```shell
quast.py test_data/contigs_1.fasta \
           test_data/contigs_2.fasta \
        -r test_data/reference.fasta.gz \
        -g test_data/genes.txt \
        -1 test_data/reads1.fastq.gz -2 test_data/reads2.fastq.gz \
        -o quast_test_output
```

QUAST stands for QUality ASsessment Tool. It evaluates genome/metagenome assemblies by computing various metrics. The current QUAST toolkit includes the general QUAST tool for genome assemblies, MetaQUAST, the extension for metagenomic datasets, QUAST-LG, the extension for large genomes (e.g., mammalians), and Icarus, the interactive visualizer for these tools.


## Protein truncation interference
-IDEEL


## Protein Sequence Prediction
-Prodigal

Fast, reliable protein-coding gene prediction for prokaryotic genomes.

```shell
prodigal -i my.genome.fna -o my.genes -a my.proteins.faa
prodigal -i my.metagenome.fna -o my.genes -a my.proteins.faa -p meta
```


## Homopolymer Counting
-Counterr

Counterr is a light-weight command line tool that computes errors in sequencing data by comparing the reads to a reference genome.

```shell
counterr -bam file.bam -genome assembly.fa -output_dir output -verbose
```

Counterr takes as input an alignment of the reads (bam + bai) and the corresponding reference (fasta) and outputs summaries (figures/tables) of errors in the reads. The tool computes both context independent and context dependent error statistics of the reads. The latter can be particularly helpful for uncovering systematic errors latent in the sequencing technology (due to hardware, basecaller, etc.).


## Antimicrobial Resistance
-EPI2ME ARMA?