MetaAMRSpotter: Automated workflow with shell scripting for uncovering hidden AMR hotspots from metagenomes

Vidya Niranjan, Chandrashekar K, Anagha S Setlur, M Purushotham Rao, S Pooja

Published: 2024-05-21 DOI: 10.17504/protocols.io.e6nvw1jyzlmk/v1

Disclaimer

This protocol can be run on Linux & Ubuntu systems with enough RAM and memory (for databases and tools) to enable appropriate data run and generation.

Abstract

This protocol employs a novel, open-source automated pipeline scripted entirely in shell for analyzing metagenomic data from various samples. Designed to streamline the workflow, the pipeline integrates functionalities for pathogen identification, antimicrobial resistance (AMR) gene detection, and listing the probable antibiotics to which the genes are resistant. This user-friendly pipeline eliminates the need for manual tools installation and configuration, simplifying the analysis process. It directly analyzes raw sequencing reads, if there is presence of appropriate reference genomes and runs through the pipeline for each sample. This protocol runs nine tools together, with just one input given at the start of the program. Demonstrated using publicly available data on both a desktop Linux system and a high-performance computing cluster, the pipeline acknowledges potential variations arising from different software tools and versions, providing users the flexibility to modify them as needed. This approach offers a robust solution for metagenomic data analysis from varied samples, facilitating efficient and accurate detection and uncovering of hidden AMR hotspots.

Keywords: AMR gene prediction, metagenomics, automated pipeline, shell scripting

Before start

All necessary tools and databases must be downloaded and installed.
Make sure the reference file for the respective selected organism is indexed and placed in the reference folder.

Steps

RETRIEVAL OF METAGENOME SAMPLES

The metagenomic samples of various organisms can be retrieved from NCBI SRA. Respective organisms' reference genomes can also be downloaded from NCBI Ref-Seq and indexed.

The following command can be used for indexing the reference:

#Indexing 
bowtie2-build <reference.fasta> <index_name>

DIRECTORY SPECIFICATION AND UNZIPPING FILES

Specify the directory of the file and unzip all .gz files in the directory. Check if the file is found first and if yes, then proceed to the next step.

Code provided below:

Note

# Specify thedirectory where your .gz files are located # Unzip all .gzfiles in the directoryfor file in *.gz;do # Check if the file exists and is a .gz file if [ -e "

file..." gunzip "

file" fidone

RUNNING FASTQC

This tool describes the quality of the raw sequence data which is a result of high through-put sequencing techniques. The tool measures length distribution, GC content and level of duplications. Quality score for the sequence which has the potency to have low-quality regions will be detected and the tool also analyzes the adapter sequence and overexpressed k-mers which could lead to errors.

The following code was used to run FastQC for selected genomes.

Note

# Check if the "fastqc" directory exists,and create it if notif [ ! -d "fastqc" ]; then mkdir fastqcfi # Create an array to store unique prefixesprefixes=() # Loop through all FASTQ and FASTQ.gz files in thecurrent directoryfor input_file in *.fastq; do # Extract theprefix from the input file name prefix=

input_file" | sed -E's/_[12].(fastq|fastq.gz)//') # Check if theprefix is already in the array if [[ ! "

prefix " ]]; then prefixes+=("

You can't use 'macro parameter character #' in math mode prefix")</span><span> </span><span> # Find theinput files for this prefix</span><span> input_r1="

{prefix}_1.fastq" input_r2="

You can't use 'macro parameter character #' in math mode{prefix}_2.fastq"</span><span> </span><span>    # Defineoutput file names based on the prefix</span><span>   output_r1_paired="result/

{prefix}/trimmomatic/1_paired.fastq" output_r2_paired="result/

{prefix}/trimmomatic/1_unpaired.fastq" output_r2_unpaired="result/

{prefix}/fastqc" mkdir -p"result/

prefix" # Run fastqc on the input files and save the resultsin the "fastqc" directory fastqc "

{prefix}/fastqc/" -t 16 fastqc"

{prefix}/fastqc/" -t 16 echo "FastQC analysis completed for$prefix."

RUNNING TRIMMOMATIC

This tool is designed to pre-process the next-generation sequencing data. Trimming will enhance the quality of the file. The tool supports both single-end read and paired-end reads data. The tool eliminates low-quality reads which optimizes ensuring all the high-quality data are retained.

The below code runs Trimmomatic.

Note

# Define Trimmomatic command with desired parameters java -jartrimmomatic-0.39.jar PE -phred33 "

input_r2""

output_r1_unpaired""

output_r2_unpaired" ILLUMINACLIP:TruSeq3-PE.fa:2:30:10LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 # Check theexit status of the Trimmomatic command if [

prefix." else echo"Trimmomatic encountered an error for $prefix." fi

ALIGNMENT USING BOWTIE2

Bowtie aligns sequences against the references and it supports gapped, local and paired-end alignment. It generates genome index using a technique called Burrows-Wheeler Transform (BWA) via similar algorithms such as Needleman-Wunch and Smith Waterman algorithms. This tool will optimize the sequence read by alignment process.

The provided code runs Bowtie2 tool.

Note

# Define the indexfile index="reference/indexed_file_name" output_dir="result/

output_dir" input_file_1="result/

{prefix}/trimmomatic/2_paired.fastq" # Check if input files were found if [ -z "

input_file_2" ]; then echo "Input files not found in thecurrent directory." exit 1 fi # Run bowtie2 bowtie2 -x "

input_file_1" -2 "

output_dir/unmapped.fastq" # Check the exitstatus of bowtie2 if [ $? -eq 0 ]; then echo "bowtie2 alignment completedsuccessfully." else echo "bowtie2 encountered anerror."Fi

SPADES ASSEMBLY

SPAdes is a genome assembly tool which works on de Bruijn Graph algorithm that reconstructs the entire genomes sequence by reading the fragments. It provides simplified graphs for the user. The tool measures the distance between k-mers and adjusts the scores to accurate distances. The contig file generated has valuable information and are high-quality assemblies that are optimized and analyzed for sequenced data.

The below code runs the Spades assembly.

Note

echo "SPAdes assembly with error correction" forward="result/

{prefix}/mapping/unmapped.2.fastq" output_path="result/

You can't use 'macro parameter character #' in math mode {prefix}/assembly/spades/"</span><span> </span><span> # Run SPAdes</span><span> spades.py--meta --pe1-1 "

forward" --pe1-2 "

output_path" --phred-offset 33 -t 8 -m 16 --only-error-correction # Rename andcompress output files output_dir="

output_dir/unmapped.1.00.0_0.cor.fastq.gz" output_file_2="

{prefix}/assembly/spades/corrected/unmapped.1.00.0_0.cor.fastq.gz" reverse="result/

{prefix}/assembly/spades/contigs" # Run SPAdes spades.py--meta --pe1-1 "

reverse" -o"

You can't use 'macro parameter character #' in math mode output_path" -t 16 --only-assembler</span><span> </span><span> # Rename andcompress output files</span><span> output_dir="

output_path/corrected/contigs" output_file="$output_dir/contigs.fasta" echo "SPAdes assembly is done"

RUNNING QUAST

The tool abbreviation stands for Quality Assessment Tool to analyze the genome assembled. The tool compares the sequence by either comparing with the available reference genome to identify the gaps in the contigs or performs de novo comparison without the reference genome and predicts the assembly quality. This tool optimizes the assembled file and predicts the low gene-coverage and provides possible results with tables and graphs.

PROFILING USING METAPHLAN

Metaphlan is the most diversely used computational tool to perform microbial profiling. It mainly focuses on metagenomic shot-gun sequencing data. The database has pre-defined markers specific to the microbial community and the sequencing data aligns against the database. The assigned reads are taxonomical labels to the given samples. It provides insights in composition and diversity of microbial populations. It is essential to determine the microbes in agriculture, health-disease, pathology and food production.

The codes for running Quast and Metaphlan are provided below:

Note

QUAST & METAPHLAN CODE:# Rename and compress output files output_dir="

output_dir/contigs.fasta" echo"SPAdes assembly is done" echo "QUAST" input="result/

{prefix}/quast" # Run quast quast.py"

output_path" --min-contig 100 echo"QUAST is done" echo"Metaphlan" output_dir="result/

output_dir" output_file_2="result/

{prefix}/metaphlan/profiled.txt" input_r1_pattern="result/

{prefix}/trimmomatic/2_paired.fastq" # Check ifinput files were found metaphlan"

input_r2_pattern" --bowtie2out"

output_file_1" echo"Metaphlan is done"

IDENTIFICATION OF AMR GENES

ABRICATE AMR

Abricate is a computational tool to identify antimicrobial resistance genes and virulence genes. Microbial genome is considered as the input. The tool uses database which contains AMR genes and is specific to sequences associated to resistance to antibiotics. Followed by comparison of sequence against the reference to determine the high similarity to sequence in AMR gene database.

10.

ABRICATE PLASMID FINDER

This tool is used to identify the plasmids in the bacterial genome that could have adapted antimicrobial resistance genes, this serves as reference to identify the similarity. The report generated has names corresponding to the matched plasmid and the alignment coverage scores.

11.

ABRICATE VIRULENCE FACTOR

This is used to find similarity against virulence factor using a pre-built database. The report generated consist of specific virulence factors' names, the alignment coverage and the virulence factor associated. The virulence factors indicate the risk of pathogenicity specific to the bacterium and this helps to understand the crucial need of developing potential therapies.

Codes for running abricate and detection of AMR genes:

Note

OVERALL ABRICATE CODE TO BE RUN:echo"abricate" input="result/

{prefix}/abricate" mkdir -p "

{prefix}/abricate/amr.txt" output_pf="result/

{prefix}/abricate/vf.txt" log_amr="result/

{prefix}/abricate/pf.log" log_vf="result/

input" > "

log_amr" echo "Abricate Plasmidfinder" abricate --threads 16 --mincov 60 --dbplasmidfinder "

output_pf" 2>"

input" > "

log_vf" echo "abricate is done" echo "abricate summary" input_amr="result/

{prefix}/abricate/pf.txt" input_vf="result/

{prefix}/abricate/amr_summary.txt" output_pf="result/

{prefix}/abricate/vf_summary.txt" echo "Abricate AMR Summary" abricate --summary "

output_amr" echo "Abricate PlasmidfinderSummary" abricate --summary "

output_pf" echo "Abricate VirulencefactorSummary" abricate --summary "

output_vf" echo "abricate summary is done"

STITCHING THESE CODES TOGETHER - FORMULATING WHOLE PIPELINE

12.

These individual codes for each tool were stitched together thereby automating the entire protocol for easy use. All users who would like to use this protocol may choose to stitch the code to run the workflow.

EXPECTED OUTCOMES - HUMAN, POULTRY AND GOAT DEMO

13.

This shell scripted workflow has been run for human genomes, poultry and goat to identify the AMR genes and the possible antibiotics they are resistant to. Expected outcomes are provided below.

A	B	C	D	E	F	G	H	I	J	K	L	M	N	O
#FILE	SEQUENCE	START	END	STRAND	GENE	COVERAGE	COVERAGE_MAP	GAPS	%COVERAGE	%IDENTITY	DATABASE	ACCESSION	PRODUCT	RESISTANCE
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_12457_length_561_cov_2.610672	2	402	-	dfrA40	1-401/513	============...	0/0	78.17	83.04	ncbi	NG_148594.1	trimethoprim-resistant dihydrofolate reductase DfrA40	TRIMETHOPRIM
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_1391_length_1533_cov_4.531800	367	1155	-	aadA27	1-789/789	===============	0/0	100	98.61	ncbi	NG_054660.1	ANT(3'')-II family aminoglycoside nucleotidyltransferase AadA27	SPECTINOMYCIN;STREPTOMYCIN
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_1775_length_1342_cov_4.503497	399	1235	+	aph(6)-Id	1-837/837	===============	0/0	100	99.88	ncbi	NG_047464.1	aminoglycoside O-phosphotransferase APH(6)-Id	STREPTOMYCIN
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_2203_length_1199_cov_3.133741	55	879	-	blaOXA-1044	1-825/825	===============	0/0	100	96.48	ncbi	NG_079234.1	OXA-211 family carbapenem-hydrolyzing class D beta-lactamase OXA-1044	BETA-LACTAM
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_391_length_3677_cov_10.529266	2203	3390	+	tet(39)	1-1188/1188	===============	0/0	100	100	ncbi	NG_048137.1	tetracycline efflux MFS transporter Tet(39)	TETRACYCLINE
result/ERR4083685/assembly/spades/contigs/contigs.fasta	NODE_46_length_16755_cov_30.079760	384	894	+	dfrA44	1-511/513	===============	0/0	99.61	86.5	ncbi	NG_073446.2	trimethoprim-resistant dihydrofolate reductase DfrA44	TRIMETHOPRIM

AMR genes detected in Goat metagenome sample

A	B	C	D	E	F	G	H	I	J	K	L	M	N	O
#FILE	SEQUENCE	START	END	STRAND	GENE	COVERAGE	COVERAGE_MAP	GAPS	%COVERAGE	%IDENTITY	DATABASE	ACCESSION	PRODUCT	RESISTANCE
result/ERR5295139/assembly/spades/contigs/contigs.fasta	NODE_10225_length_1216_cov_3.423773	211	1077	-	aadE	1-867/867	===============	0/0	100	99.89	ncbi	NG_047378.1	aminoglycoside 6-adenylyltransferase AadE	STREPTOMYCIN
result/ERR5295139/assembly/spades/contigs/contigs.fasta	NODE_1051_length_7106_cov_4.404482	1337	1978	+	catD	1-639/639	========/======	04-May	99.84	99.07	ncbi	NG_047622.1	type A-11 chloramphenicol O-acetyltransferase CatD	CHLORAMPHENICOL
result/ERR5295139/assembly/spades/contigs/contigs.fasta	NODE_10859_length_1170_cov_2.605381	71	613	+	sat4	1-543/543	===============	0/0	100	100	ncbi	NG_048072.1	streptothricin N-acetyltransferase Sat4	STREPTOTHRICIN
result/ERR5295139/assembly/spades/contigs/contigs.fasta	NODE_11304_length_1141_cov_5.878453	92	955	+	aadS	1-864/864	===============	0/0	100	99.65	ncbi	NG_047380.1	aminoglycoside 6-adenylyltransferase AadS	STREPTOMYCIN

Human AMR genes from human metagenome sample

A	B	C	D	E	F	G	H	I	J	K	L	M	N	O
#FILE	SEQUENCE	START	END	STRAND	GENE	COVERAGE	COVERAGE_MAP	GAPS	%COVERAGE	%IDENTITY	DATABASE	ACCESSION	PRODUCT	RESISTANCE
result/SRR6323357/assembly/spades/contigs/contigs.fasta	NODE_10852_length_623_cov_1.850352	1	623	-	aph(6)-Id	136-758/837	..============.	0/0	74.43	100	ncbi	NG_047465.1	aminoglycoside O-phosphotransferase APH(6)-Id	STREPTOMYCIN
result/SRR6323357/assembly/spades/contigs/contigs.fasta	NODE_13674_length_514_cov_2.904139	137	477	-	lnu(C)	156-495/495	....====/======	01-Jan	68.69	98.83	ncbi	NG_047924.1	lincosamide nucleotidyltransferase Lnu(C)	LINCOSAMIDE
result/SRR6323357/assembly/spades/contigs/contigs.fasta	NODE_14335_length_498_cov_1.162528	155	498	+	aac(6')_E64	1-344/435	============...	0/0	79.08	100	ncbi	NG_242168.1	aminoglycoside 6'-N-acetyltransferase AAC(6')-E64	AMIKACIN;KANAMYCIN;TOBRAMYCIN
result/SRR6323357/assembly/spades/contigs/contigs.fasta	NODE_15256_length_476_cov_2.216152	1	476	-	catA9	9-484/651	============...	0/0	73.12	99.79	ncbi	NG_047564.1	type A-9 chloramphenicol O-acetyltransferase	CHLORAMPHENICOL

Detection of poultry AMR genes from poultry metagenome sample

The reference genomes must be taken according to the genome in question being studied.

Citation

1. As noted above, the AMR genes can be detected for various samples, with the antibiotics they are resistant to. 2. This workflow may be applied to diverse range of samples to uncover any hidden AMR hotspots.

CONCLUSION

14.

This study thus introduces an open-source pipeline for streamlined and quick analysis of metagenomic data from various samples. Scripted entirely in shell, it integrates pathogen identification, AMR gene detection, and antibiotic resistance prediction. The pipeline directly analyzes raw reads whose quality checks have been completed priorly, eliminating manual tool setup and simplifying workflows. Demonstrated on diverse samples, it offers flexibility for customization and facilitates efficient AMR gene detection. Thus, this workflow may be applied to diverse range of samples to uncover any hidden AMR hotspots.

ACKNOWLEDGEMENTS

15.

We would like to thank Dr. Akshatha Prasanna, Assistant Professor, Department of Biotechnology, Dayananda Sagar College of Engineering for her inspiring work that led us to this study. The authors are also extremely grateful to Mr. Akshay Uttarkar, Research Scholar at RV College of Engineering, for providing all his valuable inputs.

Special thanks to our research interns Vasupradha SH, Shreya Vinod and Rajnee Joel for helping the authors run the protocol for different samples.

MetaAMRSpotter: Automated workflow with shell scripting for uncovering hidden AMR hotspots from metagenomes

Disclaimer

Abstract

Before start

Steps

RETRIEVAL OF METAGENOME SAMPLES

DIRECTORY SPECIFICATION AND UNZIPPING FILES

RUNNING FASTQC

RUNNING TRIMMOMATIC

ALIGNMENT USING BOWTIE2

SPADES ASSEMBLY

RUNNING QUAST

PROFILING USING METAPHLAN

IDENTIFICATION OF AMR GENES

STITCHING THESE CODES TOGETHER - FORMULATING WHOLE PIPELINE

EXPECTED OUTCOMES - HUMAN, POULTRY AND GOAT DEMO

CONCLUSION

ACKNOWLEDGEMENTS

推荐阅读