more features to align interleaved files (with parallelization)

bin/download_bacmet.txt

@@ -0,0 +1,13 @@
+How to download HTTP URL for the BacMet protein database.
+
+$ wget http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta
+
+IMPORTANT: databases are updated recurrently, so always check the availability of the http
+
+Download HTTPS URL for the bacmet.py script, necessary for using DIAMOND against the BacMet database with makefile.diamond and, eventually, parsing such output with makefile.bacmet.
+
+$ wget https://gist.github.com/08ab264f0ea127cc242d99bc3844c743.git 
+
+IMPORTANT: have bacmet.py inside one of the directories listed in the $PATH variable.
+
+N.B.: the python script "bacmet.py" requires python 3.6*-3.8*

bin/download_deeparg.txt

@@ -0,0 +1,29 @@
+IMPORTANT: how to download the deeparg database
+
+1) $ git clone https://bitbucket.org/gusphdproj/deeparg-ss.git
+
+then, go to the "database" directory, choose the version (v1 or v2)
+and link the file "features.fasta" into the DIAMOND work directory
+
+2) install deeparg, e.g. via conda
+
+$ conda create -n deeparg_env python=2.7.18
+$ conda activate deeparg_env
+$ pip install deeparg==1.0.2
+$ deeparg download_data -o /path/to/local/directory/
+then look for the "features.fasta" file as above.
+
+More about the deeparg installation on bitbucket:
+https://bitbucket.org/gusphdproj/deeparg-ss/src/master/
+More about the database versions in this thread:
+https://bench.cs.vt.edu/argminer/#/forum/selected_question;id=1582104387832
+
+Download HTTPS URL for the deeparg.py script, necessary for using DIAMOND against the DeepARG database with makefile.diamond, and eventually parsing such output with makefile.deeparg
+
+$ wget https://gist.github.com/5f3b236bd36acf7a5a7fcf3cff99a56d.git 
+
+IMPORTANT: after running mirror_deeparg.py, have deeparg.py inside one of the directories listed in the $PATH variable.
+N.B.: the python script "deeparg.py" requires python 3.6*-3.8*
+
+
+

lib/make/makefile.bacmet

@@ -0,0 +1,65 @@
+# makefile.bacmet
+#
+# Library makefile to make running DIAMOND simpler against the bacmet database.
+#
+# Include the full path of this file in your Makefile and making databases and
+# running the aligner can be done via make. See documentation in in
+# the biomakefile Git repository: doc/makefile.md and documentation for
+# individual targets below.
+#
+# Author: diego.brambilla@tim.it
+
+SHELL := /bin/bash
+
+# *** Parameters ***
+
+# Override in your Makefile by setting a parameter *after* the row that
+# includes this file, see documentation in doc/makefile.md.
+
+# bacmet.py specific macro
+# DIAMOND output associate some scores to each match, 
+# including the aminoacidic (a.a.) alignement length.
+# Defines the desired minimum a.a. alignemnt length threshold that bacmet.py will use, default is 37.5
+ALENGTH = 37.5
+
+
+# *** Internal ***
+
+# MAKECALL_BACMET is a macro that defines what will be output to the .makecall
+# file, the file that records versions, file stamps, parameters etc.
+#
+# *Don't redefine!*
+MAKECALL_DATE     = echo "`date +"%Y-%m-%d %H:%M:%S"`" > $@.makecall
+MAKECALL_INFILES     = echo "	Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall
+MAKECALL_BACMET = $(MAKECALL_DATE); $(MAKECALL_INFILES)
+
+
+# *** BacMet targets ***
+
+# parses DIAMOND m8 matrix to get unique tab-separated hits
+sort_best_m8: $(subst .bacmet.m8,.bacmet.sort.best.m8, $(wildcard *.bacmet.m8))
+
+# filters out DIAMOND m8 matrix to select hits based on alignment length cutoff (input)
+
+bacmet_alength_filter: $(subst .bacmet.m8,.bacmet.filter.sort.best.m8, $(wildcard *.bacmet.m8))
+
+# parses DIMAOND m8 matrix to get an annotated gene list
+sort_bacmet_genelist_tsv: $(subst .bacmet.m8,.genelist.tsv, $(wildcard *.bacmet..m8))
+
+
+# *** Make rules ***
+
+%.bacmet.sort.best.m8: %.bacmet.m8
+	$(MAKECALL_BACMET)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--"
+	sort -u -k1,1 $< | sed 's/|/\t/g' > $@
+
+%.bacmet.filter.sort.best.m8: %.bacmet.sort.best.m8
+	$(MAKECALL_BACMET)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--"
+	python bacmet.py $< $(ALENGTH) $@
+
+%.genelist.tsv: %.m8
+	$(MAKECALL_BACMET)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--"
+	cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@

lib/make/makefile.blast

@@ -19,10 +19,19 @@ SHELL := /bin/bash
 NCBI_BLAST_DB_PATH  = YOU NEED TO SET THE NCBI_BLAST_DB_PATH VARIABLE (see makefile.blast)
 BLAST_PROTEIN_DB    = YOU NEED TO SET THE BLAST_PROTEIN_DB VARIABLE (see makefile.blast)
 BLAST_NUCLEOTIDE_DB = YOU NEED TO SET THE BLAST_NUCLEOTIDE_DB VARIABLE (see makefile.blast)
+BLAST_PROTEIN_DB =  YOU NEED TO SET THE BLAST_PROTEIN_DB VARIABLE (see makefile.blast)
+#Note: you can set identity thresholds with BLASTn but not with BLASTx or BLASTp!
+BLASTX_OPTS = -max_target_seqs 1 -evalue 1e-7 -outfmt 6
+# Override in your local Makefile to custom BLASTn settings and thresholds
+BLASTN_OPTS = -max_target_seqs 1 -evalue 0.01 -perc_identity 95 -outfmt 6
+# Number of threads per process used by GNU Parallel.
+# IMPORTANT! Redefine on your local Makefile according to available resources,
+# or all_parallel_blastn6 will not work
+N_TASKS = 
 
 # *** Internal ***
 
-# MAKECALL_BLAST_BLAST is a macro that defines what will be output to the .makecall
+# MAKECALL_BLAST is a macro that defines what will be output to the .makecall
 # file, the file that records versions, file stamps, parameters etc.
 #
 # *Don't redefine!*
@@ -31,37 +40,76 @@ MAKECALL_BLAST_NCBI_DBS    = echo "	Called with path to NCBI dbs: $(NCBI_BLAST_D
 MAKECALL_BLAST_INFILES     = echo "	Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall
 MAKECALL_BLAST             = $(MAKECALL_BLAST_VERSION); $(MAKECALL_BLAST_NCBI_DBS); $(MAKECALL_BLAST_INFILES)
 
+# BLASTn-specific macros
+#
+# *Don't redefine!*
+MAKECALL_BLASTN_VERSION     = echo "$$(date +"%Y%m%d %H:%M:%S"): $@ was made with $$(blastn -version|grep 'blastn:')" > $@.makecall
+MAKECALL_BLASTN_PARAMS    = echo "	Called with parameters: $(BLASTN_OPTS)" >> $@.makecall
+
+MAKECALL_BLASTN_DB  = echo "  BLAST database: $(BLAST_NUCLEOTIDE_DB), timestamp: $(ls -lL $(BLAST_NUCLEOTIDE_DB).blastdb)" >> $@.makecall
+MAKECALL_BLAST_INFILES     = echo "	Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall
+MAKECALL_BLASTN             = $(MAKECALL_BLASTN_VERSION); $(MAKECALL_BLASTN_PARAMS); $(MAKECALL_BLASTN_DB); $(MAKECALL_BLAST_INFILES)
+
+# BLASTx-specific macros
+#
+# *Don't redefine!*
+MAKECALL_BLASTX_VERSION     = echo "$$(date +"%Y%m%d %H:%M:%S"): $@ was made with $$(blastx -version|grep 'blastx:')" > $@.makecall
+MAKECALL_BLASTX_PARAMS    = echo "	Called with parameters: $(BLASTX_OPTS)" >> $@.makecall
+
+MAKECALL_BLASTX_DB  = echo "  BLAST database: $(BLAST_PROTEIN_DB), timestamp: $(ls -lL $(BLAST_PROTEIN_DB).blastdb)" >> $@.makecall
+MAKECALL_BLAST_INFILES     = echo "	Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall
+MAKECALL_BLASTX             = $(MAKECALL_BLASTX_VERSION); $(MAKECALL_BLASTX_PARAMS); $(MAKECALL_BLASTX_DB); $(MAKECALL_BLAST_INFILES)
+
+
 # *** Targets ***
 
-# UNTESTED and assumes you've also included makefile.misc that can convert a fastq.gz to an fna:
-all_mega.blastn6s: $(subst .fna,.mega.blastn6.tsv.gz,$(wildcard *.fastq.gz))
+# UNTESTED and assumes you've previously converted fastq.gz files to .fna through makefile.misc:
+all_mega.blastn6s: $(subst .fna,.mega.blastn6.tsv.gz,$(wildcard *.fna))
+
+# Running BLASTn on uncompressed FASTA files with GNU Parallel
+# Returns a .m8 tabular file (BLAST default output format, BLOSUM62 matrix)
+all_parallel_blastn6: $(subst .fna,blastn.m8,$(wildcard *.fna))
+
+%blastn.m8: %.fna
+	$(MAKECALL_BLASTN)
+	cat $< | parallel -j $(N_TASKS) --block 100k --recstart '>' --pipe blastn $(BLASTN_OPTS) -db $(BLAST_NUCLEOTIDE_DB) -query - > $@
+	echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
+
+# Running BLASTx on uncompressed FASTA files with GNU Parallel
+# Returns a .m8 tabular file (BLAST default output format, BLOSUM62 matrix)
+all_parallel_blastx6: $(subst .fna,.blastx.m8,$(wildcard *.fna))
+
+%.blastx.m8: %.fna
+	$(MAKECALL_BLASTX)
+	cat $< | parallel -j $(N_TASKS) --block 100k --recstart '>' --pipe blastx $(BLASTX_OPTS) -db $(BLAST_PROTEIN_DB) -query - > $@
+	echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
 
 %.refseq_protein.blasttab.gz: %.fna
 	$(MAKECALL_BLAST)
-	blastx -db $(NCBI_BLAST_DB_PATH)/refseq_protein -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
+	blastx -db $(NCBI_BLAST_DB_PATH)/refseq_protein -query $< -outfmt 6 $(BLASTX_OPTS) | gzip -c > $@
 	echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
 
 %.refseq_protein.blasttab.gz: %.faa
 	$(MAKECALL_BLAST)
-	blastp -db $(NCBI_BLAST_DB_PATH)/refseq_protein -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
+	blastp -db $(NCBI_BLAST_DB_PATH)/refseq_protein -query $< -outfmt 6 $(BLASTX_OPTS) | gzip -c > $@
 	echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
 
 %.nt.mega.blastn6.tsv.gz: %.fna
 	$(MAKECALL_BLAST)
-	blastn -task megablast -db $(NCBI_BLAST_DB_PATH)/nt -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
+	blastn -task megablast -db $(NCBI_BLAST_DB_PATH)/nt -query $< -outfmt 6 $(BLASTN_OPTS) | gzip -c > $@
 	echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
 
 %.blastp6.tsv.gz: %.faa
-	blastp -db $(BLAST_PROTEIN_DB) -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
+	blastp -db $(BLAST_PROTEIN_DB) -query $< -outfmt 6 $(BLASTX_OPTS) | gzip -c > $@
 
 %.blastx.blastp6.tsv.gz: %.fna
-	blastx -db $(BLAST_PROTEIN_DB) -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
+	blastx -db $(BLAST_PROTEIN_DB) -query $< -outfmt 6 $(BLASTX_OPTS) | gzip -c > $@
 
 %.mega.blastout.gz: %.fna
-	blastn -task megablast -db $(BLAST_NUCLEOTIDE_DB) -query $< $(BLAST_OPTS) | gzip -c > $@
+	blastn -task megablast -db $(BLAST_NUCLEOTIDE_DB) -query $< $(BLASTN_OPTS) | gzip -c > $@
 
 %.mega.blastout: %.fna
-	blastn -task megablast -db $(BLAST_NUCLEOTIDE_DB) -query $< $(BLAST_OPTS) > $@
+	blastn -task megablast -db $(BLAST_NUCLEOTIDE_DB) -query $< $(BLASTN_OPTS) > $@
 
 %.mega.blastn6.tsv.gz: %.fna
 	blastn -task megablast -db $(BLAST_NUCLEOTIDE_DB) -query $< -outfmt 6 $(BLAST_OPTS) | gzip -c > $@
@@ -87,3 +135,8 @@ all_mega.blastn6s: $(subst .fna,.mega.blastn6.tsv.gz,$(wildcard *.fastq.gz))
 %.blastdb: %.fna.gz
 	gunzip -c $< | makeblastdb -out $(basename $@) -title $(basename $@) -dbtype nucl
 	touch $@
+# Create .md5 files out of compressed BLAST tablar files (.m8.gz)
+all_m8s.gz.md5s: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz))
+
+%.m8.gz.md5: %.m8.gz
+	md5sum $< > $@

lib/make/makefile.bowtie2

@@ -50,6 +50,7 @@ all_pe_sams: $(subst .r1.fastq.gz,.pe.bowtie2.sam,$(wildcard *.r1.fastq.gz))
 
 all_pe_sorted_bams: $(subst .r1.fastq.gz,.pe.bowtie2.sorted.bam,$(wildcard *.r1.fastq.gz))
 
+
 all_pe_sams.gz: $(subst .r1.fastq.gz,.pe.bowtie2.sam.gz,$(wildcard *.r1.fastq.gz))
 
 all_se_sams: $(subst .r1.fastq.gz,.se.bowtie2.sam,$(wildcard *.r1.fastq.gz))
@@ -58,7 +59,8 @@ all_se_sorted_bams: $(subst .r1.fastq.gz,.se.bowtie2.sorted.bam,$(wildcard *.r1.
 
 all_se_sams.gz: $(subst .fastq.gz,.se.bowtie2.sam.gz,$(wildcard *.fastq.gz))
 
-all_interleaved_sorted_bams: $(subst .int.fastq.gz,.int.bowtie2.sorted.bam,$(wildcard *.int.fastq.gz))
+all_interleaved_sorted_bams: $(subst .int.fastq.gz,.int.bowtie2.sorted.bam,$(wildcard *.int.fastq.gz)) $(subst .int.fna.gz,.int.bowtie2.sorted.bam,$(wildcard *.int.fna.gz))
+
 
 %.pe.bowtie2.sam: %.r1.fastq.gz %.r2.fastq.gz
 	@$(MAKECALL_BOWTIE2_ALIGN)
@@ -95,6 +97,11 @@ all_interleaved_sorted_bams: $(subst .int.fastq.gz,.int.bowtie2.sorted.bam,$(wil
 	bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) --interleaved $(word 1,$^) 2> $@.err | samtools view -Sb | samtools sort > $@
 	@echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
 
+%.int.bowtie2.sorted.bam: %.int.fna.gz
+	@$(MAKECALL_BOWTIE2_ALIGN)
+	bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) --interleaved $(word 1,$^) 2> $@.err | samtools view -Sb | samtools sort > $@
+	@echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall
+
 # SILVA mapping targets
 all_silvaparcs: $(subst .r1.fastq.gz,.pe.bowtie2-silvaparc.sam.gz,$(wildcard *.r1.fastq.gz))
 

lib/make/makefile.deeparg

@@ -0,0 +1,66 @@
+# makefile.deeparg
+#
+# Library makefile to make running DIAMOND simpler against the deeparg database.
+#
+# Include the full path of this file in your Makefile and making databases and
+# running the aligner can be done via make. See documentation in in
+# the biomakefile Git repository: doc/makefile.md and documentation for
+# individual targets below.
+#
+# Author: diego.brambilla@tim.it
+
+SHELL := /bin/bash
+
+# *** Parameters ***
+
+# Override in your Makefile by setting a parameter *after* the row that
+# includes this file, see documentation in doc/makefile.md.
+
+#deeparg.py specific macro
+# DIAMOND output associate some scores to each match, 
+# including the aminoacidic (a.a.) alignement length.
+# Defines the desired minimum a.a. alignemnt length threshold that deeparg.py will use, default is 37.5
+ALENGTH = 37.5
+
+
+# *** Internal ***
+
+# MAKECALL_DEEPARG is a macro that defines what will be output to the .makecall
+# file, the file that records versions, file stamps, parameters etc.
+#
+# *Don't redefine!*
+MAKECALL_DATE     = echo "`date +"%Y-%m-%d %H:%M:%S"`" > $@.makecall
+MAKECALL_INFILES     = echo "	Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall
+MAKECALL_DEEPARG = $(MAKECALL_DATE); $(MAKECALL_INFILES)
+
+
+# *** DeepARG targets ***
+
+# parses DIAMOND m8 matrix to get unique tab-separated hits
+# N.B.: need to use DIAMOND with the -f 6 flag
+sort_best_m8: $(subst .argdb.m8,.argdb.sort.best.m8, $(wildcard *.argdb.m8))
+
+#filters out DIAMOND m8 matrix to select hits based on alignment length cutoff (input)
+deeparg_alength_filter: $(subst .argdb.m8,.argdb.filter.sort.best.m8, $(wildcard *.argdb.m8))
+
+# parses DIMAOND m8 matrixto get an annotated gene list
+# N.B.: need to set OUTFMT=6 
+sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb.m8))
+
+
+# *** Make rules ***
+
+%.argdb.sort.best.m8: %.argdb.m8
+	$(MAKECALL_DEEPARG)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--"
+	sort -u -k1,1 $< | sed 's/|/\t/g' > $@
+
+%.argdb.filter.sort.best.m8: %.argdb.sort.best.m8
+	$(MAKECALL_DEEPARG)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--"
+	python deeparg.py $< $(ALENGTH) $@
+
+%.genelist.tsv: %.m8
+	$(MAKECALL_DEEPARG)
+	echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--"
+	cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@