Skip to content

3 Preprocessing workflow

Jump to bottom
Antton Alberdi edited this page Mar 20, 2021 · 18 revisions

Once reference genomes have been prepared with preparegenomes, you need to begin the data processing with an initial preprocessing step that will:

  1. Remove adapter remains
  2. Filter out low-quality reads
  3. Remove exact duplicate sequences
  4. Map the reads to the reference genome(s)
  5. Split reads between genomic and metagenomic fraction.

3.1- Prepare the environment

If you haven't done this previously, it is essential to declare the relevant variables before proceeding with the rest of steps.

#Declare paths
workdir=/home/projects/ku-cbd/people/antalb/holoproject
holodir=/home/projects/ku-cbd/people/antalb/holoflow

#Go to project directory
cd ${workdir}

3.2- Get the data from ERDA

Raw data should be stored in ERDA until processing it with HoloFlow to avoid storing large amounts of data in Computerome 2. Storage in C2 is VERY EXPENSIVE!!. We advice to transfer the data from ERDA to the rawdata directory.

cd ${workdir}/rawdata

Then, use sftp to connect to ERDA and the command get to transfer the files to the rawdata directory.

3.3- Prepare batches

If the data to be processed exceeds 100GB approx (e.g. one BGI lane, 1/10 of a Illumina Novaseq S4), it is desirable to run the HoloFlow PPR workflow in multiple batches to avoid memory issues. The thin nodes of C2 have a maximum memory capacity of 190GB. The qsub job submission scripts in this tutorial request 180GB of memory. Any job that requires more than that amount of memory will be automatically killed. To prepare the batches we can declare a new variable called batch.

batch=batch1

Then, for each of the batches, you can declare a different name, e.g. batch2, batch3, etc.

3.4- Prepare HoloFlow PPR files

To run HoloFlow PPR you will need to prepare two or three files.

  • PPR.txt contains the information about the sample names and data files you want to process.
  • PPR.sh is the shell script you will use to launch HoloFlow PPR.

If you need to overwrite the default configuration options, for example because you need to remove BGI adaptors rather than Illumina adaptors, you will also need to edit a config file.

  • PPR.yaml contains detailed settings of the PPR workflow.

We advice to store PRG.txt and PPR.yaml files in the input folder, and the PRG.sh file in the jobs folder.

3.4.1- Prepare PPR input file (input/PPR_${batch}.txt)

The file PPR_${batch}.txt must contain three space-separated columns with the following information:

  • Column 1: name of the sample. It must be unique.
  • Column 2: absolute path to the forward read.
  • Column 3: absolute path to the reverse read.

The document can be created and edited with the text editor nano.

nano ${workdir}/input/PPR_${batch}.txt

And it should look like this:

C05M5OP /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAH03_EKDL210001545-1a-AK6539-AK6724_HY3FCDSXY_L1_1.fq.gz /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAH03_EKDL210001545-1a-AK6539-AK6724_HY3FCDSXY_L1_2.fq.gz
C05M5HT /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAB42_EKDL210001545-1a-AK5845-AK6676_HY3FCDSXY_L1_1.fq.gz /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAB42_EKDL210001545-1a-AK5845-AK6676_HY3FCDSXY_L1_2.fq.gz
C05M5HR /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAC12_EKDL210001545-1a-AK6526-AK2651_HY3FCDSXY_L1_1.fq.gz /home/projects/ku-cbd/people/antalb/holoproject/rawdata/BHAC12_EKDL210001545-1a-AK6526-AK2651_HY3FCDSXY_L1_2.fq.gz

If using nano, type "Ctrl+X", then "Y" and finally "enter" to save the file.

Always check input file has been properly saved and review its contents.

cat ${workdir}/input/PPR_${batch}.txt

3.4.2- Prepare PPR config file (input/PPR_${batch}.yaml)

3.2.3- Prepare PPR job file (jobs/PPR_${batch}.sh)

The document can be created and edited again with the text editor nano.

nano ${workdir}/jobs/PPR_${batch}.sh

Prepare the contents in a simple text editor (Atom, Brackets, Textedit...) and paste it in the console.

python ${holodir}/preprocessing.py -f ${workdir}/input/PPR_${batch}.txt -g ${workdir}/GgalHsap.tar.gz -d ${workdir} -t 40

Check the job file has been properly saved and review its contents

cat ${workdir}/jobs/PRG.sh

3.5- Launch PPR job

Whenever you run HoloFlow it is important to unload modules to avoid conflicts with the modules that HoloFlow will automatically load. Find further explanations about it in the document Reference genome(s) preparation workflow.

Usually this job will take a few hours, so a walltime of 1 day (walltime=1:00:00:00) will suffice.

module load tools anaconda3/4.4.0
qsub -V -A ku-cbd -W group_list=ku-cbd -v "workdir=${workdir},holodir=${holodir},batch=${batch}" -d `pwd` -e ${workdir}/jobs/PPR_${batch}.err -o ${workdir}/jobs/PPR_${batch}.out -l nodes=1:ppn=40,mem=180gb,walltime=1:00:00:00 -N PPR_${batch} ${workdir}/jobs/PPR_${batch}.sh

3.6- Check results and remove intermediate files

The PPR workflow will create the following directories:

  • PPR_00-InputData: contains re-formatted input files, with simplified headers to avoid issues in downstream steps.
  • PPR_01-QualityFiltered: contains quality-filtered files.
  • PPR_02-DuplicatesRemoved: contains files with exact duplicates removed.
  • PPR_03-MappedToReference: contains ref.bam files containing the reads mapped to the reference genome(s), .fastq files containing the metagenomic fraction and .stats files containing the read filtering statistics. THESE ARE THE FILES YOU WILL NEED FOR DOWNSTREAM ANALYSES.

If the directory PPR_03-MappedToReference contains all expected files with reasonable sizes, and the HoloFlow PRG logfile does not print any error, then you can remove the intermediate files from the rest of the folders.

rm rawdata/*
rm PPR_00-InputData/${batch}*
rm PPR_01-QualityFiltered/${batch}*
rm PPR_02-DuplicatesRemoved/${batch}*

Clone this wiki locally