Sign inSign up

streptomyces/riprep

By streptomyces

Updated 5 months ago

MEME search for repeats within proteins.

Image
Data science
0

839

streptomyces/riprep repository overview

MEME searches on the output of norod.sh

Pull the docker image

docker pull streptomyces/riprep

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

# Example usage on Linux
cd /home/tom/work/
docker run -it -v ${PWD}:/home/mnt streptomyces/riprep

# Example usage on MS Windows.
docker run -it -v C:/Users/tom/work:/home/mnt streptomyces/riprep

If you already have the Pfam database somewhere then you can bind mount that directory to /home/mnt/pfam on the container side when you start a container. For example, below is how I do it during testing.

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

Do not change the /home/mnt part. This 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. You can place your input list in the mounted host directory on the host side and access it in /home/mnt/ on the container side. See the example Run on your own list below.

Running norod.sh inside the container

If you have a set of already processed genbank files then you can go direction to MEME Searches below. i.e. it is not necessary to run norod.sh.

Following the docker run command above, to ensure that norod.sh is working correctly, you can run a small test analysis on the accessions that are included in a test file named minitest.txt. Use the following command:

./norod.sh minitest.txt
Run on your own list

Use the following command to analyse your own list, substituting in a relevant filename for te_accessions.txt:

./norod.sh /home/mnt/te_accessions.txt
Output files

The output consists 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.

MEME and Pfam Searches

MEME and Pfam searches on the files written by norod.sh in /home/mnt/orgnamegbk/. These files do not have to come from the output of norod.sh. Genbank files produced by any variant of RiPPER will do.

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 Pfam database is expected to be in /home/mnt/pfam and it is expected to be named Pfam-A.hmm. You can place it in this directory either from the host side or, in a running container, do the following.

cd /home/mnt/
mkdir pfam
cd 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/
ls /home/mnt/pfam # To confirm the presence of the HMM database files.

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

njobs=4; # 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}/riprep.out

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

para-riprep
# 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-riprep 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-riprep | parallel --jobs $njobs # The actual run.
date >> timing.txt # Append finish time to timing.txt.

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

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

Make a table containing the MEME output (in the file riprep.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 riprep_pfam_merge.pl -out /home/mnt/riprep_pfam_merge.out \
-- /home/mnt/riprep.out
print_header.pl -out /home/mnt/riprep_pfam_merge.out -- rip_pfam.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/riprep/.

source code/builds.sh

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

cd; cd dockertest/riprep
docker run -it -v ${PWD}:/home/mnt streptomyces/riprep

Tag summary

Content type

Image

Digest

sha256:1ecf0dd62

Size

1.2 GB

Last updated

5 months ago

docker pull streptomyces/riprep