bwa as aligner

RQCP.sh

@@ -18,47 +18,78 @@ echo "________________________________________________________________________"
 echo ""
 echo "##### HARDWARE #####"
 echo "--------------------"
+
 physicalcpu=$(sysctl -n hw.physicalcpu)   # Get physical cpu
 echo "Physical CPU: ${physicalcpu}"       # Print physical cpu
+
 logicalcpu=$(sysctl -n hw.logicalcpu)     # Get logical cpu
 echo "Logical CPU: ${logicalcpu}"         # Print logical cpu
+
 hwmemsize=$(sysctl -n hw.memsize)         # Get memory size
 ramsize=$(expr $hwmemsize / $((1024**3))) # 1024**3 = GB
 echo "System Memory: ${ramsize} GB"       # Print RAM size
+
 timestamp=$(date +'%Y-%m-%d_%H-%M')
 echo "Date: ${timestamp}"
+
 echo "________________________________________________________________________"
 
 ##### Working directory #####
 echo ""
 echo "##### WORKING DIRECTORY #####"
 echo "-----------------------------"
+
 workdir=${0%/*} # for MacOS users, running script with double-click on "RQCP.sh"
 #workdir="." # for Linux users (todo), running script with terminal command "bash ./RQCP.sh"
 echo "CWD: ${workdir}"
-echo "________________________________________________________________________"
 
-##### Get project name ##### 
-echo ""
-echo "##### GET PROJECT NAME #####"
-echo "----------------------------"
-project="project" 
-read -p "Please enter a project name: " project
-echo "Project name: ${project}"
 echo "________________________________________________________________________"
 
+###### Get threads number #####
+#echo ""
+#echo "##### GET THREADS NUMBER ####"
+#echo "--------------------------"
+#
+#threads=$logicalcpu 
+#read -p "Please enter threads number (default, all: ${logicalcpu}): " input_threads
+#if [[ $input_threads -gt 0 ]] && [[ $input_threads -le $logicalcpu ]] 2> /dev/null
+#then
+#    threads=$input_threads
+#fi
+#echo "Pipeline running with ${threads} threads"
+#
+#echo "________________________________________________________________________"
+#
+###### Get project name ##### 
+#echo ""
+#echo "##### GET PROJECT NAME #####"
+#echo "----------------------------"
+#
+#project="project"
+#read -p "Please enter a project name: " input_project
+#if [[ $input_project != $project ]] 2> /dev/null
+#then
+#    project=$input_project
+#fi
+#echo "Project name: ${project}"
+#
+#echo "________________________________________________________________________"
+
 ###### Create links for raw samples #####
 echo ""
 echo "##### CREATE RAW READ LINKS #####"
 echo "---------------------------------"
+
 mkdir ${workdir}/results/ 2> /dev/null
 mkdir ${workdir}/results/reads/ 2> /dev/null
 mkdir ${workdir}/results/reads/raw_links/ 2> /dev/null
+
 # Create Symbolic links for raw reads
+
 for fastq in ${workdir}/resources/reads/*.fastq.gz ; do
-    ln -s $fastq ${workdir}/results/reads/raw_links/$(basename $fastq) 2> /dev/null ;
+    ln -s $fastq ${workdir}/results/reads/raw_links/$(basename $fastq) ;
 done
-echo "Silence mode"
+
 echo "________________________________________________________________________"
 
 ###### Rename links #####
@@ -66,54 +97,114 @@ echo "________________________________________________________________________"
 echo ""
 echo "##### RENAME FASTQ FILES #####"
 echo "------------------------------"
-rename -v 's/_S\d+_/_/' ${workdir}/results/reads/raw_links/*.fastq.gz 2> /dev/null                # Remove barcode-ID like {_S001_}
-rename -v 's/_L\d+_/_/' ${workdir}/results/reads/raw_links/*.fastq.gz 2> /dev/null                # Remove line-ID ID like {_L001_}
-rename -v 's/_\d+.fastq.gz/.fastq.gz/' ${workdir}/results/reads/raw_links/*.fastq.gz 2> /dev/null # Remove end-name ID like {_001}.fastq.gz
-echo "Silence mode"
+
+#rename -v 's/_S\d+_/_/' ${workdir}/results/reads/raw_links/*.fastq.gz                 # Remove barcode-ID like {_S10_}
+#rename -v 's/_L00\d_/_/' ${workdir}/results/reads/raw_links/*.fastq.gz                # Remove line-ID ID like {_L001_}
+rename -v 's/_S\d+_L00\d_/_/' ${workdir}/results/reads/raw_links/*.fastq.gz            # Remove barcode-ID and line-ID like {_S10_L001}
+rename -v 's/_001.fastq.gz/.fastq.gz/' ${workdir}/results/reads/raw_links/*.fastq.gz  # Remove end-name ID like {_001}.fastq.gz
+
+# Create dir for QC-tools and Conf-backup #####
+mkdir ${workdir}/results/reads/cutadapt_removed/
+mkdir ${workdir}/results/reads/sickle-trim_paired-end_filtered/
+mkdir ${workdir}/results/reads/fastq-join_single-end/
+mkdir ${workdir}/results/reads/all_merged_single-end/
+cp -R ${workdir}/config/ ${workdir}/results/config/
+
 echo "________________________________________________________________________"
 
-###### Create dir for QC tools #####
-mkdir ${workdir}/results/reads/cutadapt_removed/ 2> /dev/null
-mkdir ${workdir}/results/reads/sickle-trim_paired-end_filtered/ 2> /dev/null
-mkdir ${workdir}/results/reads/fastq-join_single-end/ 2> /dev/null
-mkdir ${workdir}/results/reads/all_merged_single-end/ 2> /dev/null
 
 ###### Extract bwa indexes for small genomes #####
 echo ""
 echo "##### BOWTIE2 INDEXES EXTRACTION #####"
 echo "---------------------------------"
-#tar --keep-old-files -xzvf ${workdir}/resources/databases/bowtie2_small_genomes.tar.gz -C ${workdir}/resources/databases/ # 2> /dev/null
-tar --keep-old-files -xzvf ${workdir}/resources/databases/bowtie2_large_genomes.tar.gz -C ${workdir}/resources/databases/ # 2> /dev/null
+
+#tar --keep-old-files -xzvf ${workdir}/resources/databases/bowtie2_small_genomes.tar.gz -C ${workdir}/resources/databases/
+#tar --keep-old-files -xzvf ${workdir}/resources/databases/bowtie2_large_genomes.tar.gz -C ${workdir}/resources/databases/
+tar --keep-old-files -xzvf ${workdir}/resources/databases/bwa_small_genomes.tar.gz -C ${workdir}/resources/databases/
+tar --keep-old-files -xzvf ${workdir}/resources/databases/bwa_large_genomes.tar.gz -C ${workdir}/resources/databases/
+
 echo "________________________________________________________________________"
 
 ###### Call snakemake pipeline #####
 echo ""
 echo "##### SNAKEMAKE PIPELIE #####"
 echo "-----------------------------"
+
 echo ""
 echo "Unlocking working directory:"
-snakemake --directory ${workdir} --use-conda --keep-going --rerun-incomplete -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --cores --unlock # unlock first, if previous error
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --printshellcmds \
+    --keep-going \
+    --rerun-incomplete \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --cores \
+    --unlock # unlock first, if previous error
+
 echo ""
 echo "Dry run:"
-snakemake --directory ${workdir} --use-conda --keep-going --rerun-incomplete -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --cores --dryrun # then dry run, check error like missing file
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --printshellcmds \
+    --keep-going \
+    --rerun-incomplete \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --cores \
+    --dryrun # then dry run, check error like missing file
+
 echo ""
 echo "Let's go!"
-snakemake --directory ${workdir} --use-conda --keep-going --rerun-incomplete -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --cores          # at last, real run 
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --printshellcmds \
+    --keep-going \
+    --rerun-incomplete \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --cores 1 # at last, real run 
+
 echo "________________________________________________________________________"
 
 ###### Create usefull graphs and summary ######
 echo ""
 echo "##### SNAKEMAKE PIPELINE GRAPHS ######"
 echo "--------------------------------------"
+
 mkdir ${workdir}/results/graphs/ 2> /dev/null
-snakemake --directory ${workdir} --use-conda -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --dag       | dot -Tpdf > ${workdir}/results/graphs/DAGraph.pdf
-snakemake --directory ${workdir} --use-conda -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --rulegraph | dot -Tpdf > ${workdir}/results/graphs/RuleGraph.pdf
-snakemake --directory ${workdir} --use-conda -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --filegraph | dot -Tpdf > ${workdir}/results/graphs/FileGraph.pdf
-snakemake --directory ${workdir} --use-conda -s ${workdir}/workflow/rules/reads_quality_control_pipeline.smk --summary   >             ${workdir}/results/graphs/Summary.txt
+
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --dag | dot -Tpdf > ${workdir}/results/graphs/DAGraph.pdf
+
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --rulegraph | dot -Tpdf > ${workdir}/results/graphs/RuleGraph.pdf
+
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --filegraph | dot -Tpdf > ${workdir}/results/graphs/FileGraph.pdf
+
+snakemake \
+    --directory ${workdir} \
+    --use-conda \
+    --snakefile ${workdir}/workflow/rules/reads_quality_control_pipeline.smk \
+    --summary > ${workdir}/results/graphs/Summary.txt
+
 echo "________________________________________________________________________"
 
 ###### Rename results directory with timestamp and project name ######
 echo ""
 echo "##### SCRIPT END ######"
 echo "-----------------------"
-mv ${workdir}/results/ ${workdir}/results_${timestamp}_${project}/ 2> /dev/null
+
+#mv ${workdir}/results/ ${workdir}/results_${timestamp}_${project}/ 2> /dev/null
+
+echo "________________________________________________________________________"

config/config.yaml

----
 ###############################################################################
 # Author: Nicolas Fernandez
 # Affiliation: IRD_TransVIHMI
@@ -36,9 +35,6 @@ datasets:
     #- 'fastq-join_single-end'
     #- 'fastq-join_paired-end'
     #- 'all_merged'
-  quality-tool:
-    - 'fastqc'
-    - 'fastq-screen'
 
 ## CUTADAPT ------------------------------------------------------------------------------------------
 cutadapt:
@@ -73,10 +69,10 @@ fastq-join:
 ## FASTQSCREEN ---------------------------------------------------------------------------------------
 fastq-screen:
   config: 'config/fastq-screen.conf' # path to the fastq-screen configuration file
-  subset: '100000' # Don't use the whole sequence file, but create a temporary dataset of this specified number of read (default: '100000', set '0' for all dataset)
+  subset: '1000' # Don't use the whole sequence file, but create a temporary dataset of this specified number of read (default: '100000', set '0' for all dataset)
   aligner:
-    #- 'bwa'     # Burrows-Wheeler Aligner, for mapping low-divergent sequences against a large reference genome
+    - 'bwa'     # Burrows-Wheeler Aligner, for mapping low-divergent sequences against a large reference genome
     #- 'bowtie'  # Bowtie, an ultrafast, memory-efficient short read aligner
-    - 'bowtie2' # Bowtie2, an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences
+    #- 'bowtie2' # Bowtie2, an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences
 
 ## ---------------------------------------------------------------------------------------------------

config/fastq-screen.conf

@@ -7,11 +7,12 @@
 ## Uncomment the relevant line below and set the appropriate location.
 ## Please note, this path should INCLUDE the executable filename.
 
-BOWTIE	  /usr/local/bowtie
-#BOWTIE2	  /usr/local/bowtie2-2.3.4.1/bowtie2
-BOWTIE2	  .snakemake/conda/db59fc2b/bin/bowtie2
-#BWA	  /usr/local/bwa
-BWA	  .snakemake/conda/db59fc2b/bin/bwa
+##----- BOWTIE -----
+#BOWTIE	  ${CONDA_PREFIX}/bin/bowtie
+##----- BOWTIE 2 -----
+#BOWTIE2	       ${CONDA_PREFIX}/bin/bowtie2
+##----- BWA -----
+#BWA	  ${CONDA_PREFIX}/bin/bwa
 
 ############################################
 ## Bismark (for bisulfite sequencing only) #
@@ -47,43 +48,55 @@ THREADS 1
 ##
 ##----------
 ## Human h38
+DATABASE	Human	 resources/databases/bwa/Human/Homo_sapiens_h38
 #DATABASE	Human	 resources/databases/bowtie2/Human/Homo_sapiens_h38
 ##
 ## Mouse m39
+DATABASE	Mouse	 resources/databases/bwa/Mouse/Mus_musculus_m39
 #DATABASE	Mouse	 resources/databases/bowtie2/Mouse/Mus_musculus_m39
 ##
 ## Arabidopsis thaliana - sequence from NCBI: ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.3_TAIR10/GCF_000001735.3_TAIR10_genomic.fna.gz
+DATABASE	Arabidopsis	 resources/databases/bwa/Arabidopsis/Arabidopsis_thaliana_t10
 #DATABASE	Arabidopsis	 resources/databases/bowtie2/Arabidopsis/Arabidopsis_thaliana_t10
 ##
 ## Ecoli - sequence available from EMBL accession U00096.2
-DATABASE	Ecoli	 resources/databases/bowtie2/Ecoli/Echerichia_coli_U00096
+DATABASE	Ecoli	 resources/databases/bwa/Ecoli/Echerichia_coli_U00096
+#DATABASE	Ecoli	 resources/databases/bowtie2/Ecoli/Echerichia_coli_U00096
 ##
 ##-----
 ## PhiX - sequence available from Refseq accession NC_001422.1
-DATABASE	PhiX	 resources/databases/bowtie2/Phix/PhiX
+DATABASE	PhiX	 resources/databases/bwa/Phix/PhiX
+#DATABASE	PhiX	 resources/databases/bowtie2/Phix/PhiX
 ##
 ## Adapters - sequence derived from the FastQC contaminats file found at: www.bioinformatics.babraham.ac.uk/projects/fastqc
-DATABASE	Adapters        resources/databases/bowtie2/Adapters/Adapters
+DATABASE	Adapters        resources/databases/bwa/Adapters/Adapters
+#DATABASE	Adapters        resources/databases/bowtie2/Adapters/Adapters
 ##
 ## Vector - Sequence taken from the UniVec database: http://www.ncbi.nlm.nih.gov/VecScreen/UniVec.html, BUT, without phi-X174
-DATABASE	Vectors		 resources/databases/bowtie2/Vectors/UniVec_wo_phi-X174
+DATABASE	Vectors		 resources/databases/bwa/Vectors/UniVec_wo_phi-X174
+#DATABASE	Vectors		 resources/databases/bowtie2/Vectors/UniVec_wo_phi-X174
 ##
 ##-----------
 ## Gorilla g4	
+DATABASE	Gorilla	 resources/databases/bwa/Gorilla/Gorilla_gorilla_g4
 #DATABASE	Gorilla	 resources/databases/bowtie2/Gorilla/Gorilla_gorilla_g4
 ##
 ## Chimpanzee
+DATABASE	Chimpanzee	 resources/databases/bwa/Chimpanzee/Pan_troglodytes_t3
 #DATABASE	Chimpanzee	 resources/databases/bowtie2/Chimpanzee/Pan_troglodytes_t3
 ##
 ## Bat 10
+DATABASE	Bat	 resources/databases/bwa/Bat/Pteropus_vampyrus_v1
 #DATABASE	Bat	 resources/databases/bowtie2/Bat/Pteropus_vampyrus_v1
 ##
 ## HIV - HXB2
-DATABASE	HIV	 resources/databases/bowtie2/HIV/HXB2
+DATABASE	HIV	 resources/databases/bwa/HIV/HXB2
+#DATABASE	HIV	 resources/databases/bowtie2/HIV/HXB2
 ##
 ## Ebola - sequence available from NCBI accession NC_002549
-DATABASE	Ebola	 resources/databases/bowtie2/Ebola/ZEBOV
+DATABASE	Ebola	 resources/databases/bwa/Ebola/ZEBOV
+#DATABASE	Ebola	 resources/databases/bowtie2/Ebola/ZEBOV
 ##
 ## SARS-CoV-2 - sequence from Whuhan available from NCBI accession NC_045512.2
-DATABASE	SARS-CoV-2	 resources/databases/bowtie2/SARS-CoV-2/Wuhan-WIV04-2019
-
+DATABASE	SARS-CoV-2	 resources/databases/bwa/SARS-CoV-2/Wuhan-WIV04-2019
+#DATABASE	SARS-CoV-2	 resources/databases/bowtie2/SARS-CoV-2/Wuhan-WIV04-2019

workflow/rules/reads_quality_control_pipeline.smk

@@ -5,7 +5,7 @@
 # Aim: Reads Quality Control
 # Date: 2021.04.30
 # Run: snakemake --use-conda -s reads_quality_control_pipeline.smk --cores 
-# Latest modification: 2021.10.14
+# Latest modification: 2021.10.19
 # Todo: ND
 
 ###############################################################################
@@ -17,7 +17,7 @@ configfile: "config/config.yaml"
 # WILDCARDS #
 SAMPLE, = glob_wildcards("results/reads/raw_links/{sample}_R1.fastq.gz")
 QCDIR = config["datasets"]["quality-directory"]
-
+ 
 # ENVIRONMENT #
 CUTADAPT = config["conda"]["cutadapt"]      # Cutadapt
 SICKLETRIM = config["conda"]["sickle-trim"] # Sickle-trim
@@ -69,7 +69,7 @@ rule multiqc:
        "results/reports/multiqc/{qcdir}.log"
     input:
         fastqc = "results/quality/{qcdir}/fastqc/",
-        fastqscreen = "results/quality/{qcdir}/fastqscreen/"
+        fastqscreen = "results/quality/{qcdir}/fastq-screen/"
     output:
         multiqc = directory("results/quality/{qcdir}/multiqc/")
     log:
@@ -79,7 +79,7 @@ rule multiqc:
         "--quiet "                   # -q: Only show log warning
         "--outdir {output.multiqc} " # -o: Create report in the specified output directory
         "{input.fastqc} "            # Input FastQC files
-        "{input.fastqscreen} "       # Input Fastq-Screen files
+        "{input.fastqscreen} "       # Input Fastq-Screen
         #"--no-ansi "                 # Disable coloured log
         "&> {log}"                   # Add redirection for log
 
@@ -100,7 +100,7 @@ rule fastqscreen:
     input:
         fastq = "results/reads/{qcdir}/"
     output:
-        fastqscreen = directory("results/quality/{qcdir}/fastqscreen/")
+        fastqscreen = directory("results/quality/{qcdir}/fastq-screen/")
     log:
         "results/reports/fastq-screen/{qcdir}.log"
     shell:
@@ -268,6 +268,6 @@ rule cutadapt:
         "--paired-output {output.reverse} " # -p: Write second read in a pair to FILE
         "{input.forward} "                  # Input forward reads R1.fastq
         "{input.reverse} "                  # Input reverse reads R2.fastq
-        "&> {log} "                         # Add redirection for log
+        "&> {log}"                          # Add redirection for log
 
 ###############################################################################