streptomyces/ripprep repository overview

⁠MEME searches on the output of ripper

⁠Pull the docker image

docker pull streptomyces/ripprep

Following installation, run the container using the following command, where your input accession list file is stored in /home/tom/work on Linux and MacOS systems or C:/Users/tom/work on a MS Windows system. (substitute your relevant directories in place of these):

⁠Get a container

Change to the directory where your input data is located. Then

docker run --rm -it \
-v ${PWD}:/home/mnt \
-v ${HOME}/databases/pfam:/mnt/pfam \
streptomyces/ripprep

Do not change the /home/mnt part or the /mnt/pfam part. The first refers to a directory in the container and scripts in the container expect to find this directory. The host directory you mount on /home/mnt in the container is where the output directories and files are written to. /mnt/pfam is the directory inside the container where the Pfam-A HMM database is expected to be located.

⁠The Pfam-A database

This image does not include the Pfam-A database which is needed to run the analyses. Including it would make the image very big and it will be difficult to keep the database updated to the latest release.

Your directory which appears as /mnt/pfam on the container side should contain the Pfam database inside a directory named pfam.

From inside a running container you could make the database by stepping through the commands below. This needs to be done only once everytime the Pfam-A models are updated.

# Make sure that /mnt/pfam on the container side is bound to a
# directory on the host side.
cd /mnt/pfam
wget \
'https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz'
gunzip Pfam-A.hmm.gz
hmmpress Pfam-A.hmm
cd /home/work/

The above should result in the following files in /mnt/pfam/.

Pfam-A.hmm
Pfam-A.hmm.h3f
Pfam-A.hmm.h3i
Pfam-A.hmm.h3m
Pfam-A.hmm.h3p

⁠Input files

The output of RiPPer consists of genbank files in the folder orgnamegbk. There should be one genbank file for each protein accession for which a genbank file was successfully retrieved from Genbank by RiPPer. These genbank files are the input to ripprep.

⁠Using ripprep for MEME and Pfam Searches

First we make a list of the genbank files in /home/mnt/orgnamegbk/

perl make_gbk_list.pl -outfile gbk.list -indir /home/mnt/orgnamegbk
wc -l gbk.list # Count lines in gbk.list.
head -n 6 gbk.list # Inspect the top 6 lines in gbk.list.

The variables and a function to run the searches in parallel to take advantage of multiple processors that most modern machines have.

njobs=3; # Change this to the number of processors you wish to use.
# njobs=$(nproc); # To use the number of processors on your machine.
gbklist=gbk.list
outdir=/home/mnt # Not to be changed.
errfn=${outdir}/err
ofn=${outdir}/ripprep.out

para-ripprep () {
for pf in $(seq 1 $njobs); do
 echo perl ripprep.pl -job ${pf} -jobs $njobs \
 -errfile $errfn -outfile $ofn -gbklist $gbklist
done
}

para-ripprep
# Above only prints out the commands that will be passed to
# parallel. Nothing actually runs.

Finally, run in parallel. This is done by piping the output of the function para-ripprep to parallel. Since we run more than one instance in parallel the output and error files are opened for appending. So older results in these files will not be overwritten and the newer results will be appended at the bottom. The hmmscan output goes into directories (names derived from input file names) inside the pfam_output directory. In these directories, file names are derived from protein identifiers.

rm -rf ${outdir}/pfam_output
rm  $ofn $errfn ${outdir}/bioperl.err
cp $gbklist $outdir/$gbklist
date > timing.txt # Save start time in timing.txt.
para-ripprep | parallel --jobs $njobs # The actual run.
date >> timing.txt # Append finish time to timing.txt.

Adding the header to ripprep.out ($ofn). Since multiple processes write to this output file the header can only be added after all the processes are done. The file rippmeme.header contains space separated column names.

print_header.pl -outfile $ofn -- rippmeme.header
# In this example, the command above is the same as the one below
# (commented out) because  ofn=${outdir}/ripprep.out done above.
# print_header.pl -outfile /home/mnt/ripprep.out -- rippmeme.header

Make a table containing the MEME output (in the file ripprep.out) and the top Pfam domain. In the command below the output of Pfam searches by hmmscan is assumed to be in the directory /home/mnt/pfam_output. Use to option -pfamoutdir if it is somewhere else.

perl ripprep_pfam_merge.pl -out /home/mnt/ripprep_pfam_merge.out \
-- /home/mnt/ripprep.out
print_header.pl -out /home/mnt/ripprep_pfam_merge.out -- ripppfam.header

⁠Reference for GNU Parallel

Tange, O. (2021, August 22). GNU Parallel 20210822 ('Kabul'). Zenodo. https://doi.org/10.5281/zenodo.5233953⁠

⁠Please ignore the stuff below.

On lin9, in docker/rebuild/.

source code/rebuild.sh

docker run -it -v ${PWD}:/home/mnt \
-v ${HOME}/databases/pfam:/home/mnt/pfam \
streptomyces/ripprep