From 0f0176995144555b6f0b2bc9ad1c10de4e2bb67d Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Mon, 16 Nov 2020 11:28:28 +0100 Subject: [PATCH 01/12] more recipies for interleaved files with Bowtie2 and DIAMOND --- bin/bacmet.py | 32 ++++++ bin/deeparg.py | 32 ++++++ lib/make/makefile.bowtie2 | 16 ++- lib/make/makefile.diamond | 181 +++++++++++------------------- lib/make/makefile.diamond_bacmet | 98 ++++++++++++++++ lib/make/makefile.diamond_deeparg | 114 +++++++++++++++++++ 6 files changed, 358 insertions(+), 115 deletions(-) create mode 100755 bin/bacmet.py create mode 100755 bin/deeparg.py create mode 100644 lib/make/makefile.diamond_bacmet create mode 100644 lib/make/makefile.diamond_deeparg diff --git a/bin/bacmet.py b/bin/bacmet.py new file mode 100755 index 0000000..dc4defe --- /dev/null +++ b/bin/bacmet.py @@ -0,0 +1,32 @@ +#!/usr/bin/python + +""" +Usage: + +join.arg.py [arg1] [num] [arg2] + +arg1: input tabular file with annotation +num: minimum a.a. length threshold +arg2: output file name + +""" +import io +import sys +import string +import pandas as pd +#import click + + +if len(sys.argv) == 0: + print (__doc__) + exit(0) + +# Change the header names to suit your input tabular file structure +arg=pd.read_csv(sys.argv[1], sep='\t', index_col=0, header=None, names=["read_ID", "db_id", "mrg", "features", "secondary_id", "phenotype", "pident", "alignmentlength", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qcovhsp"]) +argf=arg.loc[(arg['alignmentlength']>=float(sys.argv[2]))] +# +# +argf.to_csv(sys.argv[3],sep="\t", index=True) +# + + diff --git a/bin/deeparg.py b/bin/deeparg.py new file mode 100755 index 0000000..9225a17 --- /dev/null +++ b/bin/deeparg.py @@ -0,0 +1,32 @@ +#!/usr/bin/python + +""" +Usage: + +join.arg.py [arg1] [num] [arg2] + +arg1: input tabular file with annotation +num: minimum a.a. length threshold +arg2: output file name + +""" +import io +import sys +import string +import pandas as pd +#import click + + +if len(sys.argv) == 0: + print (__doc__) + exit(0) + +# Change the header names to suit your input tabular file structure +arg=pd.read_csv(sys.argv[1], sep='\t', index_col=0, header=None, names=["read_ID", "db_id", "features", "db", "drug", "arg","pident", "alignmentlength", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qcovhsp"]) +argf=arg.loc[(arg['alignmentlength']>=float(sys.argv[2]))] +# +# +argf.to_csv(sys.argv[3],sep="\t", index=True) +# + + diff --git a/lib/make/makefile.bowtie2 b/lib/make/makefile.bowtie2 index c68c250..eaf222c 100644 --- a/lib/make/makefile.bowtie2 +++ b/lib/make/makefile.bowtie2 @@ -48,7 +48,8 @@ MAKECALL_BOWTIE2_BUILD = $(MAKECALL_BOWTIE2_VERSION); $(MAKECALL_BOWTIE 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_sorted_bams: $(subst .r1.fastq.gz,.pe.bowtie2.sorted.bam,$(wildcard *.r1.fastq.gz)) $(subst .r1.fastq,.pe.bowtie2.sorted.bam,$(wildcard *.r1.fastq)) + all_pe_sams.gz: $(subst .r1.fastq.gz,.pe.bowtie2.sam.gz,$(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.fasta,.int.bowtie2.sorted.bam,$(wildcard *.int.fasta)) + %.pe.bowtie2.sam: %.r1.fastq.gz %.r2.fastq.gz @$(MAKECALL_BOWTIE2_ALIGN) @@ -70,6 +72,11 @@ all_interleaved_sorted_bams: $(subst .int.fastq.gz,.int.bowtie2.sorted.bam,$(wil bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $@.err | samtools view -Sb | samtools sort > $@ @echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall +%.pe.bowtie2.sorted.bam: %.r1.fastq %.r2.fastq + @$(MAKECALL_BOWTIE2_ALIGN) + bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $@.err | samtools view -Sb | samtools sort > $@ + @echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall + %.pe.bowtie2.sam.gz: %.r1.fastq.gz %.r2.fastq.gz @$(MAKECALL_BOWTIE2_ALIGN) bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $(basename $@).err | gzip -c > $@ @@ -95,6 +102,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.fasta + @$(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)) diff --git a/lib/make/makefile.diamond b/lib/make/makefile.diamond index 4c8bf1b..cee891c 100644 --- a/lib/make/makefile.diamond +++ b/lib/make/makefile.diamond @@ -1,56 +1,52 @@ -# makefile.diamond -# -# Library makefile to make running the Diamond aligner simpler. -# -# 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: erik.rikard.daniel@gmail.com +include /marconi_work/mBI20_AmbCo/diego/biomakefiles/lib/make/makefile.diamond +include /marconi_work/mBI20_AmbCo/diego/biomakefiles/lib/make/makefile.misc -SHELL := /bin/bash +# DIAMOND-specific macros +MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) +MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/arg (`ls -lL $(DIAMOND_DB_PATH)/arg.dmnd`)" >> $@.makecall +DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) +OUTFMT= 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovhsp +DIAMOND_DB_PATH = arg_db +DIAMOND_DB = /marconi_scratch/userexternal/afranzet/icme9/arg_db/arg +#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +ALENGTH = 37.5 -# *** Parameters *** -# Override in your Makefile by setting a parameter *after* the row that -# includes this file, see documentation in doc/makefile.md. +# From: https://github.com/erikrikarddaniel/biomakefiles/blob/master/doc/running_diamond_megan.md +# There doesn't seem to be a way of producing one meganized daa file from two input daa files. +# It's therefore preferable to run Diamond on an interleaved fastq file (can be gzipped) containing both reads from each pair. +# To interleave sequences, look in the makefile.misc file. -# DIAMOND_DB_PATH sets the directory where to look for generic database files -# like the protein RefSeq. The default is the current directory. -DIAMOND_DB_PATH = . -# Specify this when you want to point to a specific, non-generic, database. -# Full path, without ".dmnd". -DIAMOND_DB = +# *** Misc Targets *** -# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment -# -# See the Diamond documentation available at GitHub: -# https://github.com/bbuchfink/diamond. -# -# Look particularly at options for sensitivity, scoring and memory and -# performance. -# -# To set options, override this by inserting a row similar to the one below but -# with whatever options you'd like to add in your Makefile *after* you include -# this file. -DIAMOND_ALIGN_OPTS = +# Compress with gzip a .fastq file +#%.fastq.gz: %.fastq +# gzip -c $< > $@ + +%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz + python interleave_fastq.py $^ | gzip -c > $@ -# *** Internal *** +%.intlv.fastq: %.R1.fastq %.R2.fastq + python interleave_fastq.py $^ > $@ -# MAKECALL_DIAMOND is a macro that defines what will be output to the .makecall -# file, the file that records versions, file stamps, parameters etc. +# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). # -# *Don't redefine!* -MAKECALL_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall -MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall -MAKECALL_DB_RFP = echo " Database: $(DIAMOND_DB_PATH)/refseq_protein (`ls -lL $(DIAMOND_DB_PATH)/refseq_protein.dmnd`)" >> $@.makecall -MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall -MAKECALL_DIAMOND_RFP = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_RFP); $(MAKECALL_INFILES) -MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) +# The target uses a python script -- interleave_fastq.py -- cloned from this +# gist: +# +# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e +# *** Warning! *** +# The aforementioned script works for bot .gzip and uncompressed files +# but it has been written in PYTHON 2!!! +# +# Dependencies will automatically be unzipped if not present in unzipped format +# thanks to the pattern above. Unzipped files will be deleted after +# interleaving. +interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) -# *** Targets *** +# The same, but not zipping +interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) # Running the Diamond aligner # @@ -58,86 +54,45 @@ MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) # Non-generic databases (requires the DIAMOND_DB to be set to the full path, # except the .dmnd suffix, of the database file). -# RefSeq targets +# *** DIAMOND targets *** -%.daa: %.fna.gz - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) +%.argdb.m8: %.fastq.gz + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +#same pattern for a different target (.fna files) -%.daa: %.fna - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) +%.argdb.m8: %.fna + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -# Run diamond against the database in DIAMOND_DB for all fastq.gz files in the directory. -fastq.gz2daas: $(subst .fastq.gz,.daa,$(wildcard *.fastq.gz)) +%.argdb.sort.best.m8: %.argdb.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + sort -u -k1,1 $< > $@ -%.daa: %.fastq.gz - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +%.argdb.filter.sort.best.m8: %.argdb.sort.best.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" + python deeparg.py $< $(ALENGTH) $@ -# Nucleotides in fastq.gz file against the NCBI RefSeq protein database -fastq.gz2refseq_protein.daas: $(subst .fastq.gz,.refseq_protein.daa,$(wildcard *.fastq.gz)) +%.genelist.tsv: %.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ -%.refseq_protein.daa: %.fastq.gz - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Protein sequences in faa.gz file against the NCBI RefSeq protein database -faa2refseq_protein.daas: $(subst .faa.gz,.refseq_protein.daa,$(wildcard *.faa.gz)) $(subst .faa,.refseq_protein.daa,$(wildcard *.faa)) - -%.refseq_protein.daa: %.faa.gz - $(MAKECALL_DIAMOND_RFP) - diamond blastp -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +# *** Make rules *** +# compress all .fastq files in the work directory. +compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) -%.refseq_protein.daa: %.faa - $(MAKECALL_DIAMOND_RFP) - diamond blastp -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Nucleotides in fna.gz file against the NCBI RefSeq protein database -fna2refseq_protein.daas: $(subst .fna.gz,.refseq_protein.daa,$(wildcard *.fna.gz)) $(subst .fna,.refseq_protein.daa,$(wildcard *.fna)) - -%.refseq_protein.daa: %.fna.gz - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.refseq_protein.daa: %.fna - $(MAKECALL_DIAMOND_RFP) - diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Reformating the daa output -all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) - -%.blasttab.tsv.gz: %.daa - diamond view --daa $< -f 6 | gzip -c > $@ - -# Creating the database - -# From protein sequences in a file ending with .faa -# (Currently, Diamond only works with protein sequence databases.) -%.dmnd: %.faa.gz - $(MAKECALL_MAKEDB) - gunzip -c $< | diamond makedb -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.dmnd: %.faa - $(MAKECALL_MAKEDB) - diamond makedb --in $< -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +# parses DIAMOND m8 matrix to get unique tab-separated hits +sort_best_m8: $(subst .argdb.m8,.argdb.sort.best.m8, $(wildcard *.argdb.m8)) -all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) +#filters out DIAMOND m8 matrix to select hits based on alignment length cutoff (input) +# N.B.: depends on the python script "deeparg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) +argdb_alength_filter: $(subst .argdb.m8,.argdb.filter.sort.best.m8, $(wildcard *.argdb.m8)) -%.daa.md5: %.daa - md5sum $< > $@ +# parses DIMAOND m8 matrix to get an annotated gene list +sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb..m8)) -all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) +# Nucleotides in fastq file against the NCBI RefSeq protein database +fastq2arg_db.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) -%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz - md5sum $< > $@ diff --git a/lib/make/makefile.diamond_bacmet b/lib/make/makefile.diamond_bacmet new file mode 100644 index 0000000..725468c --- /dev/null +++ b/lib/make/makefile.diamond_bacmet @@ -0,0 +1,98 @@ +# makefile.diamond_bacmet +# # +# # Library makefile to make running DIAMOND simpler against the BacMet database. +# # +# # Include the full path of this file in your Makefile ... +# # +# # Author: diego.brambilla@lnu.se + +include makefile.diamond +include makefile.misc + +# Set this to suffixes corresponding to the types of files you wish to fetch with some +# of the t$(DIAMOND_DB)ets. +ACCEPTED_SUFFIXES = + +#In case you are not using ftp, you can override this macro and paste here the https URL. +HTTP = http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta +#IMPORTANT: databases are updated recurrently, so always check the availability of the http + +# DIAMOND-BacMet-specific macros +MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) +MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/$(DIAMOND_DB).dmnd`)" >> $@.makecall +DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) +OUTFMT= 6 + +# DIAMOND_DB_PATH sets the directory where to look for generic database files +# like the protein RefSeq. The default is the current directory. +DIAMOND_DB_PATH = . + +# Specify this when you want to point to a specific, non-generic, database. +# Full path, without ".dmnd". +DIAMOND_DB = BacMet2_EXP_database + +#bacmet.py specific macro: define minimum a.a. alignemnt length threshold +ALENGTH = 37.5 + +# *** DIAMOND BacMet targets *** + +# Fetch the BacMet protein database, and unpack the tar files. +mirror_bacmet: + wget $(HTTP) + +# Creating the database +# From protein sequences in a file ending with .fasta +# (Currently, Diamond only works with protein sequence databases.) +create_db:$(subst .fasta,.dmnd,$(wildcard *.fasta)) + +# Nucleotides in fastq/fna files against the BacMet protein database +f2bacmet_db.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fna,.bacmet.m8,$(wildcard *.fna)) + +# 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) +# N.B.: depends on the python script "bacmetg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) +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 *** + +%.dmnd: %.fasta + $(MAKECALL_MAKEDB) + diamond makedb --in $< -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.bacmet.m8: %.fastq.gz + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.bacmet.m8: %.fna + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.bacmet.sort.best.m8: %.bacmet.m8 + 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 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" + python ../../bin/bacmet.py $< $(ALENGTH) $@ + +%.genelist.tsv: %.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ + +#Create md5sum unique hashes for blast tabular files + +all_m8s.gz.md5s: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) + +%.m8.gz.md5: %.m8.gz + md5sum $< > $@ + diff --git a/lib/make/makefile.diamond_deeparg b/lib/make/makefile.diamond_deeparg new file mode 100644 index 0000000..062bf41 --- /dev/null +++ b/lib/make/makefile.diamond_deeparg @@ -0,0 +1,114 @@ +# makefile.diamond_deeparg +# # +# # Library makefile to make running DIAMOND simpler against the deeparg database. +# # +# # Include the full path of this file in your Makefile ... +# # +# # Author: diego.brambilla@lnu.se + +include makefile.diamond +include make/makefile.misc + +# Set this to suffixes corresponding to the types of files you wish to fetch with some +# of the t$(DIAMOND_DB)ets. +ACCEPTED_SUFFIXES = + +# DIAMOND-specific macros +MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) +MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/arg (`ls -lL $(DIAMOND_DB_PATH)/arg.dmnd`)" >> $@.makecall +DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) +OUTFMT= 6 + +# DIAMOND_DB_PATH sets the directory where to look for generic database files +# like the protein RefSeq. The default is the current directory. +DIAMOND_DB_PATH = arg_db + +# Specify this when you want to point to a specific, non-generic, database. +# Full path, without ".dmnd". +DIAMOND_DB = deeparg + +#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +ALENGTH = 37.5 + +# *** misc targets *** + +# compress all .fastq files in the work directory. +compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) + +# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). +# +# The target uses a python script -- interleave_fastq.py -- cloned from this +# gist: +# +# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e +# *** Warning! *** +# The aforementioned script works for bot .gzip and uncompressed files +# but it has been written in PYTHON 2!!! +# +# Dependencies will automatically be unzipped if not present in unzipped format +# thanks to the pattern above. Unzipped files will be deleted after +# interleaving. +interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) + +# The same, but not zipping +interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) + +# *** DIAMOND deeparg targets *** + +# Running the DIAMOND aligner + +# parses DIAMOND m8 matrix to get unique tab-separated hits +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) +# N.B.: depends on the python script "deeparg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) +deeparg_alength_filter: $(subst .argdb.m8,.argdb.filter.sort.best.m8, $(wildcard *.argdb.m8)) + +# parses DIMAOND m8 matrix to get an annotated gene list +sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb..m8)) + +# Nucleotides in fastq file against the NCBI RefSeq protein database +fastq2arg_db.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) + +# *** Make rules *** + +# Compress with gzip a .fastq file +%.fastq.gz: %.fastq + gzip -c $< > $@ + +%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz + python interleave_fastq.py $^ | gzip -c > $@ + +%.intlv.fastq: %.R1.fastq %.R2.fastq + python interleave_fastq.py $^ > $@ + +%.argdb.m8: %.fastq.gz + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.argdb.m8: %.fna + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.argdb.sort.best.m8: %.argdb.m8 + 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 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" + python ../../bin/deeparg.py $< $(ALENGTH) $@ + +%.genelist.tsv: %.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ + +#Create md5sum unique hashes for blast tabular files + +all_m8s.gz.md5s: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) + +%.m8.gz.md5: %.m8.gz + md5sum $< > $@ + From 484050e5dbb689ddaeafc4bc59b157099cf477db Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Mon, 16 Nov 2020 12:03:27 +0100 Subject: [PATCH 02/12] more recipies for parallel BLASTn & BLASTx --- lib/make/makefile.blast | 57 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/lib/make/makefile.blast b/lib/make/makefile.blast index 477c6a0..c3ce422 100644 --- a/lib/make/makefile.blast +++ b/lib/make/makefile.blast @@ -19,6 +19,17 @@ 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_NUCLEOTIDE_DB_PATH = . +BLAST_PROTEIN_DB_PATH = . +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 *** @@ -31,11 +42,50 @@ 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_PATH)/$(BLAST_NUCLEOTIDE_DB), timestamp: $(ls -lL $(BLAST_NUCLEOTIDE_DB_PATH)/$(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_PATH)/$(BLAST_PROTEIN_DB), timestamp: $(ls -lL $(BLAST_PROTEIN_DB_PATH)/$(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)) +# 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_PATH)/$(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_PATH)/$(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 > $@ @@ -77,4 +127,9 @@ 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 $@ + 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 $< > $@ From b825119b5194911904b082932974d4da8c09ea8c Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Mon, 16 Nov 2020 12:16:12 +0100 Subject: [PATCH 03/12] adding GNU Parallel recomendations in makefile.last --- lib/make/makefile.last | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/make/makefile.last b/lib/make/makefile.last index 1d83ad3..67348f8 100644 --- a/lib/make/makefile.last +++ b/lib/make/makefile.last @@ -37,7 +37,9 @@ LASTAL_OPTS = -e200 # more than a few percent of a core. The problems does not seem to be due to too # much IO, since I can start *many* parallel processes without cpu usage per # process going down. -PARALLEL_OPTS = -P 0 +# IMPORTANT! Do not change -P from the default (0), it will cause misbehaviours +# instead, set -j (numer of CPUS per parallel process) +PARALLEL_OPTS = -j 1 # If you use something else than the sanger flavour of fastq you should reset # this one, see the LAST documentation: http://last.cbrc.jp/. Not needed for From dc95d894ea866c2e4657d038cd65cbea99276f93 Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Tue, 17 Nov 2020 11:12:23 +0100 Subject: [PATCH 04/12] restoring makefile.diamond to version from commit 5ff8ec --- lib/make/makefile.diamond | 181 ++++++++++++++++++++++++-------------- 1 file changed, 113 insertions(+), 68 deletions(-) diff --git a/lib/make/makefile.diamond b/lib/make/makefile.diamond index cee891c..4c8bf1b 100644 --- a/lib/make/makefile.diamond +++ b/lib/make/makefile.diamond @@ -1,52 +1,56 @@ -include /marconi_work/mBI20_AmbCo/diego/biomakefiles/lib/make/makefile.diamond -include /marconi_work/mBI20_AmbCo/diego/biomakefiles/lib/make/makefile.misc - -# DIAMOND-specific macros -MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) -MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/arg (`ls -lL $(DIAMOND_DB_PATH)/arg.dmnd`)" >> $@.makecall -DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) -OUTFMT= 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovhsp -DIAMOND_DB_PATH = arg_db -DIAMOND_DB = /marconi_scratch/userexternal/afranzet/icme9/arg_db/arg -#deeparg.py specific macro: define minimum a.a. alignemnt length threshold -ALENGTH = 37.5 - - -# From: https://github.com/erikrikarddaniel/biomakefiles/blob/master/doc/running_diamond_megan.md -# There doesn't seem to be a way of producing one meganized daa file from two input daa files. -# It's therefore preferable to run Diamond on an interleaved fastq file (can be gzipped) containing both reads from each pair. -# To interleave sequences, look in the makefile.misc file. +# makefile.diamond +# +# Library makefile to make running the Diamond aligner simpler. +# +# 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: erik.rikard.daniel@gmail.com +SHELL := /bin/bash -# *** Misc Targets *** +# *** Parameters *** -# Compress with gzip a .fastq file -#%.fastq.gz: %.fastq -# gzip -c $< > $@ +# Override in your Makefile by setting a parameter *after* the row that +# includes this file, see documentation in doc/makefile.md. -%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz - python interleave_fastq.py $^ | gzip -c > $@ +# DIAMOND_DB_PATH sets the directory where to look for generic database files +# like the protein RefSeq. The default is the current directory. +DIAMOND_DB_PATH = . -%.intlv.fastq: %.R1.fastq %.R2.fastq - python interleave_fastq.py $^ > $@ +# Specify this when you want to point to a specific, non-generic, database. +# Full path, without ".dmnd". +DIAMOND_DB = -# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). +# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment +# +# See the Diamond documentation available at GitHub: +# https://github.com/bbuchfink/diamond. # -# The target uses a python script -- interleave_fastq.py -- cloned from this -# gist: +# Look particularly at options for sensitivity, scoring and memory and +# performance. # -# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e -# *** Warning! *** -# The aforementioned script works for bot .gzip and uncompressed files -# but it has been written in PYTHON 2!!! +# To set options, override this by inserting a row similar to the one below but +# with whatever options you'd like to add in your Makefile *after* you include +# this file. +DIAMOND_ALIGN_OPTS = + +# *** Internal *** + +# MAKECALL_DIAMOND is a macro that defines what will be output to the .makecall +# file, the file that records versions, file stamps, parameters etc. # -# Dependencies will automatically be unzipped if not present in unzipped format -# thanks to the pattern above. Unzipped files will be deleted after -# interleaving. -interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) +# *Don't redefine!* +MAKECALL_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall +MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_DB_RFP = echo " Database: $(DIAMOND_DB_PATH)/refseq_protein (`ls -lL $(DIAMOND_DB_PATH)/refseq_protein.dmnd`)" >> $@.makecall +MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall +MAKECALL_DIAMOND_RFP = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_RFP); $(MAKECALL_INFILES) +MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) -# The same, but not zipping -interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) +# *** Targets *** # Running the Diamond aligner # @@ -54,45 +58,86 @@ interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) # Non-generic databases (requires the DIAMOND_DB to be set to the full path, # except the .dmnd suffix, of the database file). -# *** DIAMOND targets *** +# RefSeq targets -%.argdb.m8: %.fastq.gz - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) +%.daa: %.fna.gz + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -#same pattern for a different target (.fna files) -%.argdb.m8: %.fna - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) +%.daa: %.fna + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -%.argdb.sort.best.m8: %.argdb.m8 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" - sort -u -k1,1 $< > $@ +# Run diamond against the database in DIAMOND_DB for all fastq.gz files in the directory. +fastq.gz2daas: $(subst .fastq.gz,.daa,$(wildcard *.fastq.gz)) -%.argdb.filter.sort.best.m8: %.argdb.sort.best.m8 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" - python deeparg.py $< $(ALENGTH) $@ +%.daa: %.fastq.gz + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -%.genelist.tsv: %.m8 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" - cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ +# Nucleotides in fastq.gz file against the NCBI RefSeq protein database +fastq.gz2refseq_protein.daas: $(subst .fastq.gz,.refseq_protein.daa,$(wildcard *.fastq.gz)) -# *** Make rules *** -# compress all .fastq files in the work directory. -compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) +%.refseq_protein.daa: %.fastq.gz + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Protein sequences in faa.gz file against the NCBI RefSeq protein database +faa2refseq_protein.daas: $(subst .faa.gz,.refseq_protein.daa,$(wildcard *.faa.gz)) $(subst .faa,.refseq_protein.daa,$(wildcard *.faa)) + +%.refseq_protein.daa: %.faa.gz + $(MAKECALL_DIAMOND_RFP) + diamond blastp -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -# parses DIAMOND m8 matrix to get unique tab-separated hits -sort_best_m8: $(subst .argdb.m8,.argdb.sort.best.m8, $(wildcard *.argdb.m8)) +%.refseq_protein.daa: %.faa + $(MAKECALL_DIAMOND_RFP) + diamond blastp -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Nucleotides in fna.gz file against the NCBI RefSeq protein database +fna2refseq_protein.daas: $(subst .fna.gz,.refseq_protein.daa,$(wildcard *.fna.gz)) $(subst .fna,.refseq_protein.daa,$(wildcard *.fna)) + +%.refseq_protein.daa: %.fna.gz + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.refseq_protein.daa: %.fna + $(MAKECALL_DIAMOND_RFP) + diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Reformating the daa output +all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) + +%.blasttab.tsv.gz: %.daa + diamond view --daa $< -f 6 | gzip -c > $@ + +# Creating the database + +# From protein sequences in a file ending with .faa +# (Currently, Diamond only works with protein sequence databases.) +%.dmnd: %.faa.gz + $(MAKECALL_MAKEDB) + gunzip -c $< | diamond makedb -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.dmnd: %.faa + $(MAKECALL_MAKEDB) + diamond makedb --in $< -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -#filters out DIAMOND m8 matrix to select hits based on alignment length cutoff (input) -# N.B.: depends on the python script "deeparg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) -argdb_alength_filter: $(subst .argdb.m8,.argdb.filter.sort.best.m8, $(wildcard *.argdb.m8)) +all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) -# parses DIMAOND m8 matrix to get an annotated gene list -sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb..m8)) +%.daa.md5: %.daa + md5sum $< > $@ -# Nucleotides in fastq file against the NCBI RefSeq protein database -fastq2arg_db.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) +all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) +%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz + md5sum $< > $@ From fc3d74093d6f21a8c787fc8a4168f9d27faa058b Mon Sep 17 00:00:00 2001 From: DiegoBrambilla <36502951+DiegoBrambilla@users.noreply.github.com> Date: Tue, 17 Nov 2020 11:14:39 +0100 Subject: [PATCH 05/12] parsing makefile.blast --- lib/make/makefile.blast | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/make/makefile.blast b/lib/make/makefile.blast index c3ce422..dd30cc4 100644 --- a/lib/make/makefile.blast +++ b/lib/make/makefile.blast @@ -127,7 +127,7 @@ all_parallel_blastx6: $(subst .fna,.blastx.m8,$(wildcard *.fna)) %.blastdb: %.fna.gz gunzip -c $< | makeblastdb -out $(basename $@) -title $(basename $@) -dbtype nucl - touch $ + 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)) From 17b4776b46c0bf8747db2fc8a9f9b228d58274ed Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Wed, 8 Dec 2021 15:05:45 +0100 Subject: [PATCH 06/12] improving documentation and syntax for makefile.deeparg, makefile.bacmet and makefile.trimgalore --- lib/make/makefile.bacmet | 209 +++++++++++++++++++++++++++++++++ lib/make/makefile.deeparg | 217 +++++++++++++++++++++++++++++++++++ lib/make/makefile.trimgalore | 2 +- 3 files changed, 427 insertions(+), 1 deletion(-) create mode 100644 lib/make/makefile.bacmet create mode 100644 lib/make/makefile.deeparg diff --git a/lib/make/makefile.bacmet b/lib/make/makefile.bacmet new file mode 100644 index 0000000..9c0303d --- /dev/null +++ b/lib/make/makefile.bacmet @@ -0,0 +1,209 @@ +# 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. + +# DIAMOND_DB_PATH sets the directory where to look for generic database files +# like the protein RefSeq. The default is the current directory. +DIAMOND_DB_PATH = . + +# Specify this when you want to point to a specific, non-generic, database. +# Set DIAMOND_DB in a local Makefile to the +# full path of the database file, without the ".dmnd" suffix. +DIAMOND_DB = + +#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +ALENGTH = 37.5 + +# Set this to suffixes corresponding to the types of files you wish to fetch with some +# of the t$(DIAMOND_DB)ets. +ACCEPTED_SUFFIXES = + +# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment +# +# See the Diamond documentation available at GitHub: +# https://github.com/bbuchfink/diamond. +# +# Look particularly at options for sensitivity, scoring and memory and +# performance. +# +# To set options, override this by inserting a row similar to the one below but +# with whatever options you'd like to add in your Makefile *after* you include +# this file. +DIAMOND_ALIGN_OPTS = + +# Output format of the DIAMOND results. The default value refers to a tabular file. +# For more format options see https://github.com/bbuchfink/diamond. +OUTFMT= 6 + +# Download HTTP URL for the BacMet protein database +# In case you are not using ftp, you can override this macro and paste here the https URL. +HTTP = 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 the downstream analyses +HTTPS = https://gist.github.com/08ab264f0ea127cc242d99bc3844c743.git +# IMPORTANT: bacmet.py needs to be placed inside the directory that has the DIAMOND output +# so run mirror_bacmet.py into the directory that contains .bacmet.sort.best.m8 files +# N.B.: the python script "bacmetg.py" requires python 3.6*-3.8* + +# *** Internal *** + +# MAKECALL_DIAMOND_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_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall +MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall +MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/$(DIAMOND_DB).dmnd`)" >> $@.makecall +MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) +MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) + +# *** misc targets *** + +# compress all .fastq files in the work directory. +compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) + +# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). +# +# The target uses a python script -- interleave_fastq.py -- cloned from this +# gist: +# +# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e +# *** Warning! *** +# The aforementioned script works for bot .gzip and uncompressed files +# but it has been written in PYTHON 2!!! +# +# Dependencies will automatically be unzipped if not present in unzipped format +# thanks to the pattern above. Unzipped files will be deleted after +# interleaving. +interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) + +# The same, but not zipping +interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) + +# *** DIAMOND BacMet targets *** + +# Fetch the BacMet protein database, and unpack the tar files. +mirror_bacmet: + wget $(HTTP) + +# Fetch the bacmet.py script +# Uses git as a prerequisite, if not available you can manually download the file via wget +mirror_bacmet.py: + git clone $(HTTPS) + +# Creating the database +# From protein sequences in a file ending with .faa +# (Currently, Diamond only works with protein sequence databases.) +create_bacmetdb:$(subst .faa,.dmnd,$(wildcard *.faa)) + +# Nucleotides or proteins in fastq/fna/faa files against the BacMet protein database +f2bacmetdb.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fna,.bacmet.m8,$(wildcard *.fna)) $(subst .faa,.bacmet.m8,$(wildcard *.faa)) + +# 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 *** + +# Compress with gzip a .fastq file +%.fastq.gz: %.fastq + gzip -c $< > $@ + +%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz + python interleave_fastq.py $^ | gzip -c > $@ + +%.intlv.fastq: %.R1.fastq %.R2.fastq + python interleave_fastq.py $^ > $@ + +#Create the DIAMOND database +# From protein sequences in a file ending with .faa +# (Currently, Diamond only works with protein sequence databases.) +%.dmnd: %.faa.gz + $(MAKECALL_MAKEDB) + gunzip -c $< | diamond makedb -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.dmnd: %.faa + $(MAKECALL_MAKEDB) + diamond makedb --in $< -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) +%.bacmet.m8: %.fastq.gz + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.bacmet.m8: %.fna + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.faa files) +%.bacmet.m8: %.faa + $(MAKECALL_DIAMOND_BACMET) + diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.bacmet.sort.best.m8: %.bacmet.m8 + 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 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" + python bacmet.py $< $(ALENGTH) $@ + +%.genelist.tsv: %.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ + +# Reformating the daa output +all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) + +%.blasttab.tsv.gz: %.daa + diamond view --daa $< -f 6 | gzip -c > $@ + +#Create md5sum unique hashes for DIAMOND .daa file extension +all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) + +%.daa.md5: %.daa + md5sum $< > $@ + +#Create md5sum unique hashes for blast tabular files + +all_m8s.md5sums: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) $(subst .m8,.m8.md5,$(wildcard *.m8)) + +%.m8.gz.md5: %.m8.gz + md5sum $< > $@ + +%.m8.md5: %.m8 + md5sum $< > $@ + +all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) + +%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz + md5sum $< > $@ diff --git a/lib/make/makefile.deeparg b/lib/make/makefile.deeparg new file mode 100644 index 0000000..cf4b996 --- /dev/null +++ b/lib/make/makefile.deeparg @@ -0,0 +1,217 @@ +# 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. + +# DIAMOND_DB_PATH sets the directory where to look for generic database files +# like the protein RefSeq. The default is the current directory. +DIAMOND_DB_PATH = . + +# Specify this when you want to point to a specific, non-generic, database. +# Set DIAMOND_DB in a local Makefile to the +# full path of the database file, without the ".dmnd" suffix. +DIAMOND_DB = + +#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +ALENGTH = 37.5 + +# Set this to suffixes corresponding to the types of files you wish to fetch with some +# of the t$(DIAMOND_DB)ets. +ACCEPTED_SUFFIXES = + +# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment +# +# See the Diamond documentation available at GitHub: +# https://github.com/bbuchfink/diamond. +# +# Look particularly at options for sensitivity, scoring and memory and +# performance. +# +# To set options, override this by inserting a row similar to the one below but +# with whatever options you'd like to add in your Makefile *after* you include +# this file. +DIAMOND_ALIGN_OPTS = + +# Output format of the DIAMOND results. The default value refers to a tabular file. +# For more format options see https://github.com/bbuchfink/diamond. +OUTFMT= 6 + +# 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 +# $ source 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 the downstream analyses +HTTPS = https://gist.github.com/5f3b236bd36acf7a5a7fcf3cff99a56d.git +# IMPORTANT: deeparg.py needs to be placed inside the directory that has the DIAMOND output +# so run mirror_deeparg.py into the directory that contains .argdb.sort.best.m8 files +# N.B.: the python script "deeparg.py" requires python 3.6*-3.8* + +# *** Internal *** + +# MAKECALL_DIAMOND_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_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall +MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall +MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/*.dmnd`)" >> $@.makecall +MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) +MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) + +# *** misc targets *** + +# compress all .fastq files in the work directory. +compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) + +# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). +# +# The target uses a python script -- interleave_fastq.py -- cloned from this +# gist: +# +# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e +# *** Warning! *** +# The aforementioned script works for bot .gzip and uncompressed files +# but it has been written in PYTHON 2!!! +# +# Dependencies will automatically be unzipped if not present in unzipped format +# thanks to the pattern above. Unzipped files will be deleted after +# interleaving. +interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) + +# The same, but not zipping +interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) + +# *** DIAMOND DeepARG targets *** + +# Fetch the deeparg.py script +# Uses git as a prerequisite, if not available you can manually download the file via wget +mirror_deeparg.py: + git clone $(HTTPS) + +# Creating the database +# From protein sequences in a file ending with .faa +# (Currently, Diamond only works with protein sequence databases.) +create_deepargdb:$(subst .faa.gz,.dmnd, $(wildcard *.faa.gz)) $(subst .faa,.dmnd, $(wildcard *.faa)) + +# Nucleotides or proteins in fastq/fna/faa files against the DeepARG protein database +f2deepargdb.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) $(subst .faa,.argdb.m8,$(wildcard *.faa)) + +# parses DIAMOND m8 matrix to get unique tab-separated hits +# N.B.: need to set OUTFMT=6 +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 *** + +# Compress with gzip a .fastq file +%.fastq.gz: %.fastq + gzip -c $< > $@ + +%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz + python interleave_fastq.py $^ | gzip -c > $@ + +%.intlv.fastq: %.R1.fastq %.R2.fastq + python interleave_fastq.py $^ > $@ + +#Create the DIAMOND database +# From protein sequences in a file ending with .faa +# (Currently, Diamond only works with protein sequence databases.) +%.dmnd: %.faa.gz + $(MAKECALL_MAKEDB) + gunzip -c $< | diamond makedb -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.dmnd: %.faa + $(MAKECALL_MAKEDB) + diamond makedb --in $< -d $(basename $@) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) +%.argdb.m8: %.fastq.gz + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.argdb.m8: %.fna + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Pattern for a BLASTp against a protein target (.faa files) +%.argdb.m8: %.faa + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +%.argdb.sort.best.m8: %.argdb.m8 + 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 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" + python deeparg.py $< $(ALENGTH) $@ + +%.genelist.tsv: %.m8 + echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" + cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ + +# Reformating the daa output +all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) + +%.blasttab.tsv.gz: %.daa + diamond view --daa $< -f 6 | gzip -c > $@ + +#Create md5sum unique hashes for DIAMOND .daa file extension +all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) + +%.daa.md5: %.daa + md5sum $< > $@ + +#Create md5sum unique hashes for blast tabular files + +all_m8s.md5sums: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) $(subst .m8,.m8.md5,$(wildcard *.m8)) + +%.m8.gz.md5: %.m8.gz + md5sum $< > $@ + +%.m8.md5: %.m8 + md5sum $< > $@ + +all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) + +%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz + md5sum $< > $@ diff --git a/lib/make/makefile.trimgalore b/lib/make/makefile.trimgalore index a4797b0..e6acb02 100644 --- a/lib/make/makefile.trimgalore +++ b/lib/make/makefile.trimgalore @@ -4,7 +4,7 @@ # # Include the full path of this file in your Makefile ... # -# Author: diego.brambilla@lnu.se +# Author: diego.brambilla@tim.it SHELL := /bin/bash From 8243fd8c5d34576b7c53338f964f7ec40592ce6b Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Wed, 8 Dec 2021 15:44:02 +0100 Subject: [PATCH 07/12] merging with erikrikarddaniel/biomakefiles master, removingredundant files --- bin/bacmet.py | 32 --------- bin/deeparg.py | 32 --------- lib/make/makefile.diamond_bacmet | 98 ------------------------- lib/make/makefile.diamond_deeparg | 114 ------------------------------ 4 files changed, 276 deletions(-) delete mode 100755 bin/bacmet.py delete mode 100755 bin/deeparg.py delete mode 100644 lib/make/makefile.diamond_bacmet delete mode 100644 lib/make/makefile.diamond_deeparg diff --git a/bin/bacmet.py b/bin/bacmet.py deleted file mode 100755 index dc4defe..0000000 --- a/bin/bacmet.py +++ /dev/null @@ -1,32 +0,0 @@ -#!/usr/bin/python - -""" -Usage: - -join.arg.py [arg1] [num] [arg2] - -arg1: input tabular file with annotation -num: minimum a.a. length threshold -arg2: output file name - -""" -import io -import sys -import string -import pandas as pd -#import click - - -if len(sys.argv) == 0: - print (__doc__) - exit(0) - -# Change the header names to suit your input tabular file structure -arg=pd.read_csv(sys.argv[1], sep='\t', index_col=0, header=None, names=["read_ID", "db_id", "mrg", "features", "secondary_id", "phenotype", "pident", "alignmentlength", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qcovhsp"]) -argf=arg.loc[(arg['alignmentlength']>=float(sys.argv[2]))] -# -# -argf.to_csv(sys.argv[3],sep="\t", index=True) -# - - diff --git a/bin/deeparg.py b/bin/deeparg.py deleted file mode 100755 index 9225a17..0000000 --- a/bin/deeparg.py +++ /dev/null @@ -1,32 +0,0 @@ -#!/usr/bin/python - -""" -Usage: - -join.arg.py [arg1] [num] [arg2] - -arg1: input tabular file with annotation -num: minimum a.a. length threshold -arg2: output file name - -""" -import io -import sys -import string -import pandas as pd -#import click - - -if len(sys.argv) == 0: - print (__doc__) - exit(0) - -# Change the header names to suit your input tabular file structure -arg=pd.read_csv(sys.argv[1], sep='\t', index_col=0, header=None, names=["read_ID", "db_id", "features", "db", "drug", "arg","pident", "alignmentlength", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qcovhsp"]) -argf=arg.loc[(arg['alignmentlength']>=float(sys.argv[2]))] -# -# -argf.to_csv(sys.argv[3],sep="\t", index=True) -# - - diff --git a/lib/make/makefile.diamond_bacmet b/lib/make/makefile.diamond_bacmet deleted file mode 100644 index 725468c..0000000 --- a/lib/make/makefile.diamond_bacmet +++ /dev/null @@ -1,98 +0,0 @@ -# makefile.diamond_bacmet -# # -# # Library makefile to make running DIAMOND simpler against the BacMet database. -# # -# # Include the full path of this file in your Makefile ... -# # -# # Author: diego.brambilla@lnu.se - -include makefile.diamond -include makefile.misc - -# Set this to suffixes corresponding to the types of files you wish to fetch with some -# of the t$(DIAMOND_DB)ets. -ACCEPTED_SUFFIXES = - -#In case you are not using ftp, you can override this macro and paste here the https URL. -HTTP = http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta -#IMPORTANT: databases are updated recurrently, so always check the availability of the http - -# DIAMOND-BacMet-specific macros -MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) -MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/$(DIAMOND_DB).dmnd`)" >> $@.makecall -DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) -OUTFMT= 6 - -# DIAMOND_DB_PATH sets the directory where to look for generic database files -# like the protein RefSeq. The default is the current directory. -DIAMOND_DB_PATH = . - -# Specify this when you want to point to a specific, non-generic, database. -# Full path, without ".dmnd". -DIAMOND_DB = BacMet2_EXP_database - -#bacmet.py specific macro: define minimum a.a. alignemnt length threshold -ALENGTH = 37.5 - -# *** DIAMOND BacMet targets *** - -# Fetch the BacMet protein database, and unpack the tar files. -mirror_bacmet: - wget $(HTTP) - -# Creating the database -# From protein sequences in a file ending with .fasta -# (Currently, Diamond only works with protein sequence databases.) -create_db:$(subst .fasta,.dmnd,$(wildcard *.fasta)) - -# Nucleotides in fastq/fna files against the BacMet protein database -f2bacmet_db.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fna,.bacmet.m8,$(wildcard *.fna)) - -# 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) -# N.B.: depends on the python script "bacmetg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) -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 *** - -%.dmnd: %.fasta - $(MAKECALL_MAKEDB) - diamond makedb --in $< -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.bacmet.m8: %.fastq.gz - $(MAKECALL_DIAMOND_BACMET) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -#same pattern for a different target (.fna files) -%.bacmet.m8: %.fna - $(MAKECALL_DIAMOND_BACMET) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.bacmet.sort.best.m8: %.bacmet.m8 - 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 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" - python ../../bin/bacmet.py $< $(ALENGTH) $@ - -%.genelist.tsv: %.m8 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" - cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ - -#Create md5sum unique hashes for blast tabular files - -all_m8s.gz.md5s: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) - -%.m8.gz.md5: %.m8.gz - md5sum $< > $@ - diff --git a/lib/make/makefile.diamond_deeparg b/lib/make/makefile.diamond_deeparg deleted file mode 100644 index 062bf41..0000000 --- a/lib/make/makefile.diamond_deeparg +++ /dev/null @@ -1,114 +0,0 @@ -# makefile.diamond_deeparg -# # -# # Library makefile to make running DIAMOND simpler against the deeparg database. -# # -# # Include the full path of this file in your Makefile ... -# # -# # Author: diego.brambilla@lnu.se - -include makefile.diamond -include make/makefile.misc - -# Set this to suffixes corresponding to the types of files you wish to fetch with some -# of the t$(DIAMOND_DB)ets. -ACCEPTED_SUFFIXES = - -# DIAMOND-specific macros -MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) -MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/arg (`ls -lL $(DIAMOND_DB_PATH)/arg.dmnd`)" >> $@.makecall -DIAMOND_ALIGN_OPTS = --id 90 --evalue 0.00001 -f $(OUTFMT) -OUTFMT= 6 - -# DIAMOND_DB_PATH sets the directory where to look for generic database files -# like the protein RefSeq. The default is the current directory. -DIAMOND_DB_PATH = arg_db - -# Specify this when you want to point to a specific, non-generic, database. -# Full path, without ".dmnd". -DIAMOND_DB = deeparg - -#deeparg.py specific macro: define minimum a.a. alignemnt length threshold -ALENGTH = 37.5 - -# *** misc targets *** - -# compress all .fastq files in the work directory. -compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) - -# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). -# -# The target uses a python script -- interleave_fastq.py -- cloned from this -# gist: -# -# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e -# *** Warning! *** -# The aforementioned script works for bot .gzip and uncompressed files -# but it has been written in PYTHON 2!!! -# -# Dependencies will automatically be unzipped if not present in unzipped format -# thanks to the pattern above. Unzipped files will be deleted after -# interleaving. -interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) - -# The same, but not zipping -interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) - -# *** DIAMOND deeparg targets *** - -# Running the DIAMOND aligner - -# parses DIAMOND m8 matrix to get unique tab-separated hits -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) -# N.B.: depends on the python script "deeparg.py" that requires python 3.6*-3.8* (see also required dependencies at the top of the file) -deeparg_alength_filter: $(subst .argdb.m8,.argdb.filter.sort.best.m8, $(wildcard *.argdb.m8)) - -# parses DIMAOND m8 matrix to get an annotated gene list -sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb..m8)) - -# Nucleotides in fastq file against the NCBI RefSeq protein database -fastq2arg_db.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) - -# *** Make rules *** - -# Compress with gzip a .fastq file -%.fastq.gz: %.fastq - gzip -c $< > $@ - -%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz - python interleave_fastq.py $^ | gzip -c > $@ - -%.intlv.fastq: %.R1.fastq %.R2.fastq - python interleave_fastq.py $^ > $@ - -%.argdb.m8: %.fastq.gz - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -#same pattern for a different target (.fna files) -%.argdb.m8: %.fna - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.argdb.sort.best.m8: %.argdb.m8 - 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 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->filtering $><--" - python ../../bin/deeparg.py $< $(ALENGTH) $@ - -%.genelist.tsv: %.m8 - echo -E "`date +"%Y-%m-%d %H:%M:%S"`: --->sorting $><--" - cut -f 1,2 $< | sed 's/|/\t/g' | sort -u -k6,6 - > $@ - -#Create md5sum unique hashes for blast tabular files - -all_m8s.gz.md5s: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) - -%.m8.gz.md5: %.m8.gz - md5sum $< > $@ - From 1f7033b2a52c8d2a48b27bff70e1f591561e781e Mon Sep 17 00:00:00 2001 From: Daniel Lundin Date: Wed, 8 Dec 2021 16:27:47 +0100 Subject: [PATCH 08/12] Update makefile.bacmet --- lib/make/makefile.bacmet | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/make/makefile.bacmet b/lib/make/makefile.bacmet index 9c0303d..a647829 100644 --- a/lib/make/makefile.bacmet +++ b/lib/make/makefile.bacmet @@ -25,7 +25,7 @@ DIAMOND_DB_PATH = . # full path of the database file, without the ".dmnd" suffix. DIAMOND_DB = -#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +# deeparg.py specific macro: define minimum a.a. alignemnt length threshold ALENGTH = 37.5 # Set this to suffixes corresponding to the types of files you wish to fetch with some @@ -47,7 +47,7 @@ DIAMOND_ALIGN_OPTS = # Output format of the DIAMOND results. The default value refers to a tabular file. # For more format options see https://github.com/bbuchfink/diamond. -OUTFMT= 6 +OUTFMT = 6 # Download HTTP URL for the BacMet protein database # In case you are not using ftp, you can override this macro and paste here the https URL. @@ -66,12 +66,12 @@ HTTPS = https://gist.github.com/08ab264f0ea127cc242d99bc3844c743.git # file, the file that records versions, file stamps, parameters etc. # # *Don't redefine!* -MAKECALL_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall -MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall -MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall +MAKECALL_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall +MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/$(DIAMOND_DB).dmnd`)" >> $@.makecall MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) -MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) +MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) # *** misc targets *** @@ -118,7 +118,7 @@ f2bacmetdb.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fn # 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) +# 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)) @@ -138,7 +138,7 @@ sort_bacmet_genelist_tsv: $(subst .bacmet.m8,.genelist.tsv, $(wildcard *.bacmet. %.intlv.fastq: %.R1.fastq %.R2.fastq python interleave_fastq.py $^ > $@ -#Create the DIAMOND database +# Create the DIAMOND database # From protein sequences in a file ending with .faa # (Currently, Diamond only works with protein sequence databases.) %.dmnd: %.faa.gz From d809e8956c0db2465f6a41697f9a930797239b79 Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Thu, 9 Dec 2021 15:29:04 +0100 Subject: [PATCH 09/12] fixed macros and targets of makefile.blast --- lib/make/makefile.blast | 50 ++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 26 deletions(-) diff --git a/lib/make/makefile.blast b/lib/make/makefile.blast index 77f0dfb..9dca7f5 100644 --- a/lib/make/makefile.blast +++ b/lib/make/makefile.blast @@ -19,13 +19,11 @@ 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_NUCLEOTIDE_DB_PATH = . -BLAST_PROTEIN_DB_PATH = . 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 +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 +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 @@ -33,7 +31,7 @@ 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!* @@ -45,36 +43,36 @@ MAKECALL_BLAST = $(MAKECALL_BLAST_VERSION); $(MAKECALL_BLAST_NCBI_DB # 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_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_PATH)/$(BLAST_NUCLEOTIDE_DB), timestamp: $(ls -lL $(BLAST_NUCLEOTIDE_DB_PATH)/$(BLAST_NUCLEOTIDE_DB).blastdb)" >> $@.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) +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_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_PATH)/$(BLAST_PROTEIN_DB), timestamp: $(ls -lL $(BLAST_PROTEIN_DB_PATH)/$(BLAST_PROTEIN_DB).blastdb)" >> $@.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) +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_PATH)/$(BLAST_NUCLEOTIDE_DB) -query - > $@ + $(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 @@ -82,36 +80,36 @@ all_parallel_blastn6: $(subst .fna,blastn.m8,$(wildcard *.fna)) 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_PATH)/$(BLAST_PROTEIN_DB) -query - > $@ + $(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 > $@ From 892efc2484eb0ba295c077863945a28c22ff9c5b Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Thu, 9 Dec 2021 15:54:12 +0100 Subject: [PATCH 10/12] moved DIAMOND targets from makefile.deeparg and makefile.bacmet to makefile.diamond --- lib/make/makefile.bacmet | 150 +++---------------------------------- lib/make/makefile.deeparg | 151 +++----------------------------------- lib/make/makefile.diamond | 71 +++++++++++++++++- 3 files changed, 92 insertions(+), 280 deletions(-) diff --git a/lib/make/makefile.bacmet b/lib/make/makefile.bacmet index 9c0303d..7a4bfaa 100644 --- a/lib/make/makefile.bacmet +++ b/lib/make/makefile.bacmet @@ -16,39 +16,12 @@ SHELL := /bin/bash # Override in your Makefile by setting a parameter *after* the row that # includes this file, see documentation in doc/makefile.md. -# DIAMOND_DB_PATH sets the directory where to look for generic database files -# like the protein RefSeq. The default is the current directory. -DIAMOND_DB_PATH = . - -# Specify this when you want to point to a specific, non-generic, database. -# Set DIAMOND_DB in a local Makefile to the -# full path of the database file, without the ".dmnd" suffix. -DIAMOND_DB = - -#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +# 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 ALENGTH = 37.5 -# Set this to suffixes corresponding to the types of files you wish to fetch with some -# of the t$(DIAMOND_DB)ets. -ACCEPTED_SUFFIXES = - -# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment -# -# See the Diamond documentation available at GitHub: -# https://github.com/bbuchfink/diamond. -# -# Look particularly at options for sensitivity, scoring and memory and -# performance. -# -# To set options, override this by inserting a row similar to the one below but -# with whatever options you'd like to add in your Makefile *after* you include -# this file. -DIAMOND_ALIGN_OPTS = - -# Output format of the DIAMOND results. The default value refers to a tabular file. -# For more format options see https://github.com/bbuchfink/diamond. -OUTFMT= 6 - # Download HTTP URL for the BacMet protein database # In case you are not using ftp, you can override this macro and paste here the https URL. HTTP = http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta @@ -56,47 +29,20 @@ HTTP = http://bacmet.biomedicine.gu.se/download/BacMet2_EXP_database.fasta # Download HTTPS URL for the bacmet.py script, necessary for the downstream analyses HTTPS = https://gist.github.com/08ab264f0ea127cc242d99bc3844c743.git -# IMPORTANT: bacmet.py needs to be placed inside the directory that has the DIAMOND output -# so run mirror_bacmet.py into the directory that contains .bacmet.sort.best.m8 files +# IMPORTANT: after running mirror_bacmet.py, have bacmet.py inside one of the directories listed in the $PATH variable. # N.B.: the python script "bacmetg.py" requires python 3.6*-3.8* # *** Internal *** -# MAKECALL_DIAMOND_BACMET is a macro that defines what will be output to the .makecall +# 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_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall -MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_DATE = echo "`date +"%Y-%m-%d %H:%M:%S"`" > $@.makecall MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall -MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/$(DIAMOND_DB).dmnd`)" >> $@.makecall -MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) -MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) - -# *** misc targets *** - -# compress all .fastq files in the work directory. -compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) - -# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). -# -# The target uses a python script -- interleave_fastq.py -- cloned from this -# gist: -# -# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e -# *** Warning! *** -# The aforementioned script works for bot .gzip and uncompressed files -# but it has been written in PYTHON 2!!! -# -# Dependencies will automatically be unzipped if not present in unzipped format -# thanks to the pattern above. Unzipped files will be deleted after -# interleaving. -interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) - -# The same, but not zipping -interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) +MAKECALL_BACMET = $(MAKECALL_DATE); $(MAKECALL_INFILES) -# *** DIAMOND BacMet targets *** +# *** BacMet targets *** # Fetch the BacMet protein database, and unpack the tar files. mirror_bacmet: @@ -107,13 +53,6 @@ mirror_bacmet: mirror_bacmet.py: git clone $(HTTPS) -# Creating the database -# From protein sequences in a file ending with .faa -# (Currently, Diamond only works with protein sequence databases.) -create_bacmetdb:$(subst .faa,.dmnd,$(wildcard *.faa)) - -# Nucleotides or proteins in fastq/fna/faa files against the BacMet protein database -f2bacmetdb.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fna,.bacmet.m8,$(wildcard *.fna)) $(subst .faa,.bacmet.m8,$(wildcard *.faa)) # parses DIAMOND m8 matrix to get unique tab-separated hits sort_best_m8: $(subst .bacmet.m8,.bacmet.sort.best.m8, $(wildcard *.bacmet.m8)) @@ -128,82 +67,17 @@ sort_bacmet_genelist_tsv: $(subst .bacmet.m8,.genelist.tsv, $(wildcard *.bacmet. # *** Make rules *** -# Compress with gzip a .fastq file -%.fastq.gz: %.fastq - gzip -c $< > $@ - -%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz - python interleave_fastq.py $^ | gzip -c > $@ - -%.intlv.fastq: %.R1.fastq %.R2.fastq - python interleave_fastq.py $^ > $@ - -#Create the DIAMOND database -# From protein sequences in a file ending with .faa -# (Currently, Diamond only works with protein sequence databases.) -%.dmnd: %.faa.gz - $(MAKECALL_MAKEDB) - gunzip -c $< | diamond makedb -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.dmnd: %.faa - $(MAKECALL_MAKEDB) - diamond makedb --in $< -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) -%.bacmet.m8: %.fastq.gz - $(MAKECALL_DIAMOND_BACMET) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -#same pattern for a different target (.fna files) -%.bacmet.m8: %.fna - $(MAKECALL_DIAMOND_BACMET) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -#same pattern for a different target (.faa files) -%.bacmet.m8: %.faa - $(MAKECALL_DIAMOND_BACMET) - diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - %.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 - > $@ - -# Reformating the daa output -all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) - -%.blasttab.tsv.gz: %.daa - diamond view --daa $< -f 6 | gzip -c > $@ - -#Create md5sum unique hashes for DIAMOND .daa file extension -all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) - -%.daa.md5: %.daa - md5sum $< > $@ - -#Create md5sum unique hashes for blast tabular files - -all_m8s.md5sums: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) $(subst .m8,.m8.md5,$(wildcard *.m8)) - -%.m8.gz.md5: %.m8.gz - md5sum $< > $@ - -%.m8.md5: %.m8 - md5sum $< > $@ - -all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) - -%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz - md5sum $< > $@ diff --git a/lib/make/makefile.deeparg b/lib/make/makefile.deeparg index cf4b996..8fa1cdc 100644 --- a/lib/make/makefile.deeparg +++ b/lib/make/makefile.deeparg @@ -16,39 +16,12 @@ SHELL := /bin/bash # Override in your Makefile by setting a parameter *after* the row that # includes this file, see documentation in doc/makefile.md. -# DIAMOND_DB_PATH sets the directory where to look for generic database files -# like the protein RefSeq. The default is the current directory. -DIAMOND_DB_PATH = . - -# Specify this when you want to point to a specific, non-generic, database. -# Set DIAMOND_DB in a local Makefile to the -# full path of the database file, without the ".dmnd" suffix. -DIAMOND_DB = - -#deeparg.py specific macro: define minimum a.a. alignemnt length threshold +#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 ALENGTH = 37.5 -# Set this to suffixes corresponding to the types of files you wish to fetch with some -# of the t$(DIAMOND_DB)ets. -ACCEPTED_SUFFIXES = - -# DIAMOND_ALIGN_OPTS sets non-essential parameters for diamond alignment -# -# See the Diamond documentation available at GitHub: -# https://github.com/bbuchfink/diamond. -# -# Look particularly at options for sensitivity, scoring and memory and -# performance. -# -# To set options, override this by inserting a row similar to the one below but -# with whatever options you'd like to add in your Makefile *after* you include -# this file. -DIAMOND_ALIGN_OPTS = - -# Output format of the DIAMOND results. The default value refers to a tabular file. -# For more format options see https://github.com/bbuchfink/diamond. -OUTFMT= 6 - # 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) @@ -66,67 +39,32 @@ OUTFMT= 6 # Download HTTPS URL for the deeparg.py script, necessary for the downstream analyses HTTPS = https://gist.github.com/5f3b236bd36acf7a5a7fcf3cff99a56d.git -# IMPORTANT: deeparg.py needs to be placed inside the directory that has the DIAMOND output -# so run mirror_deeparg.py into the directory that contains .argdb.sort.best.m8 files +# 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* # *** Internal *** -# MAKECALL_DIAMOND_DEEPARG is a macro that defines what will be output to the .makecall +# 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_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall -MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall +MAKECALL_DATE = echo "`date +"%Y-%m-%d %H:%M:%S"`" > $@.makecall MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall -MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/$(DIAMOND_DB) (`ls -lL $(DIAMOND_DB_PATH)/*.dmnd`)" >> $@.makecall -MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) -MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) - -# *** misc targets *** - -# compress all .fastq files in the work directory. -compress_all_fastqs: $(subst .fastq,.fastq.gz, $(wildcard *.fastq)) - -# Interleave two fastq files (named %.R1.fastq and %.R2.fastq respectively). -# -# The target uses a python script -- interleave_fastq.py -- cloned from this -# gist: -# -# https://gist.github.com/rpucheq/cf50eeea794529ca6b0e -# *** Warning! *** -# The aforementioned script works for bot .gzip and uncompressed files -# but it has been written in PYTHON 2!!! -# -# Dependencies will automatically be unzipped if not present in unzipped format -# thanks to the pattern above. Unzipped files will be deleted after -# interleaving. -interleave_all_fastq.gzs: $(subst .R1.fastq.gz,.intlv.fastq.gz,$(wildcard *.R1.fastq.gz)) +MAKECALL_DEEPARG = $(MAKECALL_DATE); $(MAKECALL_INFILES) -# The same, but not zipping -interleave_all_fastqs: $(subst .R1.fastq,.intlv.fastq,$(wildcard *.R1.fastq)) - -# *** DIAMOND DeepARG targets *** +# *** DeepARG targets *** # Fetch the deeparg.py script # Uses git as a prerequisite, if not available you can manually download the file via wget mirror_deeparg.py: git clone $(HTTPS) -# Creating the database -# From protein sequences in a file ending with .faa -# (Currently, Diamond only works with protein sequence databases.) -create_deepargdb:$(subst .faa.gz,.dmnd, $(wildcard *.faa.gz)) $(subst .faa,.dmnd, $(wildcard *.faa)) - -# Nucleotides or proteins in fastq/fna/faa files against the DeepARG protein database -f2deepargdb.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) $(subst .faa,.argdb.m8,$(wildcard *.faa)) # parses DIAMOND m8 matrix to get unique tab-separated hits # N.B.: need to set OUTFMT=6 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 @@ -136,82 +74,17 @@ sort_argdb_genelist_tsv: $(subst .argdb.m8,.genelist.tsv, $(wildcard *.argdb.m8) # *** Make rules *** -# Compress with gzip a .fastq file -%.fastq.gz: %.fastq - gzip -c $< > $@ - -%.intlv.fastq.gz: %.R1.fastq.gz %.R2.fastq.gz - python interleave_fastq.py $^ | gzip -c > $@ - -%.intlv.fastq: %.R1.fastq %.R2.fastq - python interleave_fastq.py $^ > $@ - -#Create the DIAMOND database -# From protein sequences in a file ending with .faa -# (Currently, Diamond only works with protein sequence databases.) -%.dmnd: %.faa.gz - $(MAKECALL_MAKEDB) - gunzip -c $< | diamond makedb -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -%.dmnd: %.faa - $(MAKECALL_MAKEDB) - diamond makedb --in $< -d $(basename $@) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) -%.argdb.m8: %.fastq.gz - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -#same pattern for a different target (.fna files) -%.argdb.m8: %.fna - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - -# Pattern for a BLASTp against a protein target (.faa files) -%.argdb.m8: %.faa - $(MAKECALL_DIAMOND_DEEPARG) - diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) -f $(OUTFMT) - echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall - %.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 - > $@ - -# Reformating the daa output -all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) - -%.blasttab.tsv.gz: %.daa - diamond view --daa $< -f 6 | gzip -c > $@ - -#Create md5sum unique hashes for DIAMOND .daa file extension -all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) - -%.daa.md5: %.daa - md5sum $< > $@ - -#Create md5sum unique hashes for blast tabular files - -all_m8s.md5sums: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) $(subst .m8,.m8.md5,$(wildcard *.m8)) - -%.m8.gz.md5: %.m8.gz - md5sum $< > $@ - -%.m8.md5: %.m8 - md5sum $< > $@ - -all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) - -%.blasttab.tsv.gz.md5: %.blasttab.tsv.gz - md5sum $< > $@ diff --git a/lib/make/makefile.diamond b/lib/make/makefile.diamond index b485756..6c5fe53 100644 --- a/lib/make/makefile.diamond +++ b/lib/make/makefile.diamond @@ -37,6 +37,9 @@ DIAMOND_DB = # this file. # Set number of threads with --threads; supposedly works with both aligning and # makedb. +# Set the output format of the DIAMOND results with --outfmt. +# The default value is 6 and refers to a tabular file (.m8 extension). +# For more format options see https://github.com/bbuchfink/diamond. DIAMOND_ALIGN_OPTS = DIAMOND_MAKEDB_OPTS = @@ -49,8 +52,13 @@ DIAMOND_MAKEDB_OPTS = MAKECALL_VERSION = echo "`date +"%Y-%m-%d %H:%M:%S"`: $@ was made with `diamond --version`" > $@.makecall MAKECALL_PARAMS = echo " Called with parameters: $(DIAMOND_ALIGN_OPTS)" >> $@.makecall MAKECALL_DB_RFP = echo " Database: $(DIAMOND_DB_PATH)/refseq_protein (`ls -lL $(DIAMOND_DB_PATH)/refseq_protein.dmnd`)" >> $@.makecall +MAKECALL_DB_BACMET = echo " Database: $(DIAMOND_DB_PATH)/bacmet (`ls -lL $(DIAMOND_DB_PATH)/bacmet.dmnd`)" >> $@.makecall +MAKECALL_DB_DEEPARG = echo " Database: $(DIAMOND_DB_PATH)/deeparg (`ls -lL $(DIAMOND_DB_PATH)/deeparg.dmnd`)" >> $@.makecall + MAKECALL_INFILES = echo " Input files: $^ (`ls -lL $^|tr '\n' ','`)" >> $@.makecall MAKECALL_DIAMOND_RFP = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_RFP); $(MAKECALL_INFILES) +MAKECALL_DIAMOND_BACMET = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_BACMET); $(MAKECALL_INFILES) +MAKECALL_DIAMOND_DEEPARG = $(MAKECALL_VERSION); $(MAKECALL_PARAMS); $(MAKECALL_DB_DEEPARG); $(MAKECALL_INFILES) MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) # *** Targets *** @@ -63,6 +71,9 @@ MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) # RefSeq targets +# Run diamond against the NCBI RefSeq protein database for all fna.gz, fna or fastq.gz files in the directory. +fastq.gz2daas: $(subst .fastq.gz,.daa,$(wildcard *.fastq.gz)) + %.daa: %.fna.gz $(MAKECALL_DIAMOND_RFP) diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) @@ -73,9 +84,6 @@ MAKECALL_MAKEDB = $(MAKECALL_VERSION); $(MAKECALL_INFILES) diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall -# Run diamond against the database in DIAMOND_DB for all fastq.gz files in the directory. -fastq.gz2daas: $(subst .fastq.gz,.daa,$(wildcard *.fastq.gz)) - %.daa: %.fastq.gz $(MAKECALL_DIAMOND_RFP) diamond blastx -f 100 -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) @@ -115,6 +123,53 @@ fna2refseq_protein.daas: $(subst .fna.gz,.refseq_protein.daa,$(wildcard *.fna.gz diamond blastx -f 100 -d $(DIAMOND_DB_PATH)/refseq_protein --query $< -o $@ $(DIAMOND_ALIGN_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +# BacMet targets + +# Nucleotides or proteins in fastq/fna/faa files against the BacMet protein database +f2bacmetdb.m8s: $(subst .fastq.gz,.bacmet.m8,$(wildcard *.fastq.gz)) $(subst .fna,.bacmet.m8,$(wildcard *.fna)) $(subst .faa,.bacmet.m8,$(wildcard *.faa)) + +# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) +%.bacmet.m8: %.fastq.gz + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.bacmet.m8: %.fna + $(MAKECALL_DIAMOND_BACMET) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.faa files) +%.bacmet.m8: %.faa + $(MAKECALL_DIAMOND_BACMET) + diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# DeepARG targets + +# Nucleotides or proteins in fastq/fna/faa files against the DeepARG protein database +f2deepargdb.m8s: $(subst .fastq.gz,.argdb.m8,$(wildcard *.fastq.gz)) $(subst .fna,.argdb.m8,$(wildcard *.fna)) $(subst .faa,.argdb.m8,$(wildcard *.faa)) + +# Pattern for a BLASTx against a compressed nucleotide target (.fastq.gz) +%.argdb.m8: %.fastq.gz + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +#same pattern for a different target (.fna files) +%.argdb.m8: %.fna + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastx -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + +# Pattern for a BLASTp against a protein target (.faa files) +%.argdb.m8: %.faa + $(MAKECALL_DIAMOND_DEEPARG) + diamond blastp -d $(DIAMOND_DB) --query $< -o $@ $(DIAMOND_ALIGN_OPTS) + echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall + + # Reformating the daa output all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) @@ -135,11 +190,21 @@ all_blasttab.tsv.gzs: $(subst .daa,.blasttab.tsv.gz,$(wildcard *.daa)) diamond makedb --in $< -d $(basename $@) $(DIAMOND_MAKEDB_OPTS) echo "`date +"%Y-%m-%d %H:%M:%S"`: DONE" >> $@.makecall +# Creates md5sum unique hashes for DIAMOND .daa file extension all_daa.md5sums: $(subst .daa,.daa.md5,$(wildcard *.daa)) %.daa.md5: %.daa md5sum $< > $@ +# Creates md5sum unique hashes for blast tabular files +all_m8s.md5sums: $(subst .m8.gz,.m8.gz.md5,$(wildcard *.m8.gz)) $(subst .m8,.m8.md5,$(wildcard *.m8)) + +%.m8.gz.md5: %.m8.gz + md5sum $< > $@ + +%.m8.md5: %.m8 + md5sum $< > $@ + all_blasttab.md5sums: $(subst .blasttab.tsv.gz,.blasttab.tsv.gz.md5,$(wildcard *.blasttab.tsv.gz)) %.blasttab.tsv.gz.md5: %.blasttab.tsv.gz From f0a58d370fe8f937c64343ffd16fd41fa70a2ef6 Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Thu, 9 Dec 2021 16:07:28 +0100 Subject: [PATCH 11/12] removed targets with unzipped FASTQs from makefile.bowtie2 --- lib/make/makefile.bowtie2 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/lib/make/makefile.bowtie2 b/lib/make/makefile.bowtie2 index eaf222c..1b37d1d 100644 --- a/lib/make/makefile.bowtie2 +++ b/lib/make/makefile.bowtie2 @@ -48,7 +48,7 @@ MAKECALL_BOWTIE2_BUILD = $(MAKECALL_BOWTIE2_VERSION); $(MAKECALL_BOWTIE 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)) $(subst .r1.fastq,.pe.bowtie2.sorted.bam,$(wildcard *.r1.fastq)) +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)) @@ -59,7 +59,7 @@ 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)) $(subst .int.fasta,.int.bowtie2.sorted.bam,$(wildcard *.int.fasta)) +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 @@ -72,11 +72,6 @@ all_interleaved_sorted_bams: $(subst .int.fastq.gz,.int.bowtie2.sorted.bam,$(wil bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $@.err | samtools view -Sb | samtools sort > $@ @echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall -%.pe.bowtie2.sorted.bam: %.r1.fastq %.r2.fastq - @$(MAKECALL_BOWTIE2_ALIGN) - bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $@.err | samtools view -Sb | samtools sort > $@ - @echo "`date +"%Y%m%d %H:%M:%S"`: DONE" >> $@.makecall - %.pe.bowtie2.sam.gz: %.r1.fastq.gz %.r2.fastq.gz @$(MAKECALL_BOWTIE2_ALIGN) bowtie2 $(BOWTIE_ALIGN_OPTS) -x $(BOWTIE_INDEX) -1 $(word 1,$^) -2 $(word 2,$^) 2> $(basename $@).err | gzip -c > $@ @@ -102,7 +97,7 @@ 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.fasta +%.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 From 21bc41b7ac8b1249cb708bb5d686a113bacc0f4a Mon Sep 17 00:00:00 2001 From: Diego_Brambilla Date: Thu, 9 Dec 2021 16:26:42 +0100 Subject: [PATCH 12/12] instructions on how to fetch DeepARG and BacMet are in /bin --- bin/download_bacmet.txt | 13 +++++++++++++ bin/download_deeparg.txt | 29 +++++++++++++++++++++++++++++ lib/make/makefile.bacmet | 22 ++-------------------- lib/make/makefile.deeparg | 30 +++--------------------------- 4 files changed, 47 insertions(+), 47 deletions(-) create mode 100644 bin/download_bacmet.txt create mode 100644 bin/download_deeparg.txt diff --git a/bin/download_bacmet.txt b/bin/download_bacmet.txt new file mode 100644 index 0000000..65e18e3 --- /dev/null +++ b/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* diff --git a/bin/download_deeparg.txt b/bin/download_deeparg.txt new file mode 100644 index 0000000..81a4e95 --- /dev/null +++ b/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* + + + diff --git a/lib/make/makefile.bacmet b/lib/make/makefile.bacmet index 7a4bfaa..09b432b 100644 --- a/lib/make/makefile.bacmet +++ b/lib/make/makefile.bacmet @@ -19,18 +19,9 @@ SHELL := /bin/bash # 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 +# Defines the desired minimum a.a. alignemnt length threshold that bacmet.py will use, default is 37.5 ALENGTH = 37.5 -# Download HTTP URL for the BacMet protein database -# In case you are not using ftp, you can override this macro and paste here the https URL. -HTTP = 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 the downstream analyses -HTTPS = https://gist.github.com/08ab264f0ea127cc242d99bc3844c743.git -# IMPORTANT: after running mirror_bacmet.py, have bacmet.py inside one of the directories listed in the $PATH variable. -# N.B.: the python script "bacmetg.py" requires python 3.6*-3.8* # *** Internal *** @@ -42,17 +33,8 @@ 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 *** - -# Fetch the BacMet protein database, and unpack the tar files. -mirror_bacmet: - wget $(HTTP) - -# Fetch the bacmet.py script -# Uses git as a prerequisite, if not available you can manually download the file via wget -mirror_bacmet.py: - git clone $(HTTPS) +# *** 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)) diff --git a/lib/make/makefile.deeparg b/lib/make/makefile.deeparg index 8fa1cdc..79915a7 100644 --- a/lib/make/makefile.deeparg +++ b/lib/make/makefile.deeparg @@ -19,28 +19,9 @@ SHELL := /bin/bash #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 +# Defines the desired minimum a.a. alignemnt length threshold that deeparg.py will use, default is 37.5 ALENGTH = 37.5 -# 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 -# $ source 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 the downstream analyses -HTTPS = 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* # *** Internal *** @@ -52,16 +33,11 @@ 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 *** - -# Fetch the deeparg.py script -# Uses git as a prerequisite, if not available you can manually download the file via wget -mirror_deeparg.py: - git clone $(HTTPS) +# *** DeepARG targets *** # parses DIAMOND m8 matrix to get unique tab-separated hits -# N.B.: need to set OUTFMT=6 +# 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)