Merge pull request #151 from RoanKanninga/master

multiple fixes

README.md

@@ -1,9 +1,7 @@
 <h1>NGS_DNA pipeline</h1>
-
 <h2>Manual</h2>
 Find manual on installation and use at https://molgenis.gitbooks.io/molgenis-pipelines/
 
-
 <h2>Summary</h2>
 The sequencer is producing reads (in FastQ format) and are aligned to the hg19 reference genome with BWA (Li & Durbin <sup>1</sup>).
 Sambamba (Tarasov et al.<sup>2</sup>)  is processing the aligned reads and then we applied GATK (McKenna et al. <sup>3</sup>) duplicate removal,

generate_template.sh

@@ -2,21 +2,26 @@
 
 module load NGS_DNA/3.4.1
 module list 
-HOST=$(hostname)
+HOST=$(hostname -s)
 thisDir=$(pwd)
 
-ENVIRONMENT_PARAMETERS="parameters_${HOST%%.*}.csv"
+ENVIRONMENT_PARAMETERS="parameters_${HOST}.csv"
 TMPDIRECTORY=$(basename $(cd ../../ && pwd ))
 GROUP=$(basename $(cd ../../../ && pwd ))
 
 PROJECT=projectXX
 WORKDIR="/groups/${GROUP}/${TMPDIRECTORY}"
 RUNID=runXX
+SPECIES="homo_sapiens"
+BUILD="b37"
 
 ## Normal user, please leave BATCH at _chr
 ## For expert modus: small batchsize (6) fill in '_small'  or per chromosome fill in _chr
 BATCH="_chr"
-samplesheet=${WORKDIR}/generatedscripts/${PROJECT}/${PROJECT}.csv
+
+GENSCRIPTS="${WORKDIR}/generatedscripts/${PROJECT}"
+
+samplesheet=${GENSCRIPTS}/${PROJECT}.csv
 mac2unix $samplesheet
 
 python ${EBROOTNGS_DNA}/scripts/samplesize.py ${samplesheet} $thisDir
@@ -25,6 +30,7 @@ SAMPLESIZE=$(cat externalSampleIDs.txt | uniq | wc -l)
 python ${EBROOTNGS_DNA}/scripts/gender.py $samplesheet
 var=$(cat ${samplesheet}.tmp | wc -l)
 
+
 if [ $var != 0 ]
 then
     	mv ${samplesheet}.tmp ${samplesheet}
@@ -44,44 +50,45 @@ then
      rm .compute.properties
 fi
 
-if [ -f ${WORKDIR}/generatedscripts/${PROJECT}/out.csv  ];
+if [ -f ${GENSCRIPTS}/out.csv  ];
 then
-    	rm -rf ${WORKDIR}/generatedscripts/${PROJECT}/out.csv
+    	rm -rf ${GENSCRIPTS}/out.csv
 fi
 
-echo "tmpName,${TMPDIRECTORY}" > ${WORKDIR}/generatedscripts/${PROJECT}/tmpdir_parameters.csv 
+echo "tmpName,${TMPDIRECTORY}" > ${GENSCRIPTS}/tmpdir_parameters.csv 
 
-perl ${EBROOTNGS_DNA}/scripts/convertParametersGitToMolgenis.pl ${WORKDIR}/generatedscripts/${PROJECT}/tmpdir_parameters.csv > \
-${WORKDIR}/generatedscripts/${PROJECT}/tmpdir_parameters_converted.csv
+perl ${EBROOTNGS_DNA}/scripts/convertParametersGitToMolgenis.pl ${GENSCRIPTS}/tmpdir_parameters.csv > \
+${GENSCRIPTS}/tmpdir_parameters_converted.csv
 
 perl ${EBROOTNGS_DNA}/scripts/convertParametersGitToMolgenis.pl ${EBROOTNGS_DNA}/parameters.csv > \
-${WORKDIR}/generatedscripts/${PROJECT}/out.csv
+${GENSCRIPTS}/out.csv
 
 perl ${EBROOTNGS_DNA}/scripts/convertParametersGitToMolgenis.pl ${EBROOTNGS_DNA}/parameters_${GROUP}.csv > \
-${WORKDIR}/generatedscripts/${PROJECT}/group_parameters.csv
+${GENSCRIPTS}/group_parameters.csv
 
 perl ${EBROOTNGS_DNA}/scripts/convertParametersGitToMolgenis.pl ${EBROOTNGS_DNA}/${ENVIRONMENT_PARAMETERS} > \
-${WORKDIR}/generatedscripts/${PROJECT}/environment_parameters.csv
+${GENSCRIPTS}/environment_parameters.csv
+
 
 sh $EBROOTMOLGENISMINCOMPUTE/molgenis_compute.sh \
--p ${WORKDIR}/generatedscripts/${PROJECT}/out.csv \
--p ${WORKDIR}/generatedscripts/${PROJECT}/group_parameters.csv \
--p ${WORKDIR}/generatedscripts/${PROJECT}/environment_parameters.csv \
--p ${WORKDIR}/generatedscripts/${PROJECT}/tmpdir_parameters_converted.csv \
+-p ${GENSCRIPTS}/out.csv \
+-p ${GENSCRIPTS}/group_parameters.csv \
+-p ${GENSCRIPTS}/environment_parameters.csv \
+-p ${GENSCRIPTS}/tmpdir_parameters_converted.csv \
 -p ${EBROOTNGS_DNA}/batchIDList${BATCH}.csv \
--p ${WORKDIR}/generatedscripts/${PROJECT}/${PROJECT}.csv \
+-p ${GENSCRIPTS}/${PROJECT}.csv \
 -w ${EBROOTNGS_DNA}/create_in-house_ngs_projects_workflow.csv \
--rundir ${WORKDIR}/generatedscripts/${PROJECT}/scripts \
+-rundir ${GENSCRIPTS}/scripts \
 --runid ${RUNID} \
 -o "workflowpath=${WORKFLOW};\
-outputdir=scripts/jobs;mainParameters=${WORKDIR}/generatedscripts/${PROJECT}/out.csv;\
-group_parameters=${WORKDIR}/generatedscripts/${PROJECT}/group_parameters.csv;\
+outputdir=scripts/jobs;mainParameters=${GENSCRIPTS}/out.csv;\
+group_parameters=${GENSCRIPTS}/group_parameters.csv;\
 groupname=${GROUP};\
 ngsversion=$(module list | grep -o -P 'NGS_DNA(.+)');\
-environment_parameters=${WORKDIR}/generatedscripts/${PROJECT}/environment_parameters.csv;\
-tmpdir_parameters=${WORKDIR}/generatedscripts/${PROJECT}/tmpdir_parameters_converted.csv;\
+environment_parameters=${GENSCRIPTS}/environment_parameters.csv;\
+tmpdir_parameters=${GENSCRIPTS}/tmpdir_parameters_converted.csv;\
 batchIDList=${EBROOTNGS_DNA}/batchIDList${BATCH}.csv;\
-worksheet=${WORKDIR}/generatedscripts/${PROJECT}/${PROJECT}.csv" \
+worksheet=${GENSCRIPTS}/${PROJECT}.csv" \
 -weave \
 --generate
 

parameters.csv

@@ -199,6 +199,7 @@ projectBatchGenotypedVariantCalls,${projectPrefix}.batch-${batchID}.variant.call
 
 ### Protocols 19, 20, 21, 21, 22, 23 and 24 (SnpEff,CmdLineAnnotator, VariantAnnotator, MergeBatches,Gavin VariantFiltration, IndelFiltration, SplitIndelsSNPs) ###
 projectVariantCallsSnpEff_Annotated,${projectPrefix}.batch-${batchID}.variant.calls.snpeff.vcf
+projectVariantCallsSnpEff_SummaryHtml,${projectPrefix}.batch-${batchID}.variant.calls.snpeff_Summary.html
 projectVariantCallsSnpEff_ExAC_Annotated,${projectPrefix}.batch-${batchID}.variant.calls.snpeff.exac.vcf
 projectVariantCallsSnpEff_ExAC_GoNL_Annotated,${projectPrefix}.batch-${batchID}.variant.calls.snpeff.exac.gonl.vcf
 projectVariantCallsSnpEff_ExAC_GoNL_CADD_Annotated,${projectPrefix}.batch-${batchID}.variant.calls.snpeff.exac.gonl.cadd.vcf

parameters_umcg-atd.csv

@@ -0,0 +1 @@
+workDir,${root}/groups/umcg-atd/

protocols/Convading_XHMM_MakeControlGroup/Convading_MakeControlGroup.sh

@@ -8,7 +8,7 @@ sh Convading_MakeControlGroup.sh -i /path/to/bams -w /path/to/workdir/ -p ONCO_v
 
 Arguments
         Required:
-	-i|--bamsfolder		path to where all bams are
+	-i|--bamsfolder		path to where all bams are (WARNING: do not use inputBams as the directory name)
 	-w|--workdir		path to working directory
         -p|--panel	    	name of panel (e.g. CARDIO_v2, ONCO_v3)
 	Optional:

protocols/CopyPrmTmpData.sh

@@ -1,4 +1,8 @@
 #MOLGENIS walltime=02:00:00 mem=4gb
+
+set -e 
+set -u 
+
 #string tmpName
 #string allRawNgsTmpDataDir
 #string allRawNgsPrmDataDir
@@ -21,6 +25,7 @@
 #list barcode
 #list lane
 
+
 n_elements=${internalSampleID[@]}
 max_index=${#internalSampleID[@]}-1
 

protocols/CreateExternSamplesProjects.sh

@@ -34,8 +34,11 @@
 #list lane
 #string ngsUtilsVersion
 
+set -e 
+set -u
+
 umask 0007
-module load ${ngsUtilsVersion}
+module load $ngsUtilsVersion
 module load $ngsversion
 
 module list
@@ -67,40 +70,36 @@ for ((samplenumber = 0; samplenumber <= max_index; samplenumber++))
 do
 	if [[ ${seqType[samplenumber]} == "SR" ]]
 	then
-  		if [[ ${barcode[samplenumber]} == "None" ]]
+		if [[ ${barcode[samplenumber]} == "None" ]]
 		then
 			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz
 			ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz.md5
-  		else
-      			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz
-      			ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz.md5
+		else
+			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz
+			ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}.fq.gz.md5
 		fi
 	elif [[ ${seqType[samplenumber]} == "PE" ]]
 	then
 		if [[ ${barcode[samplenumber]} == "None" ]]
-    		then
-    			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz
+		then
+			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz
 			ln -sf ${externalFastQ_2[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz
 			ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz.md5
-        		ln -sf ${externalFastQ_2[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz.md5
-		else          
-        		ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz
-        		ln -sf ${externalFastQ_2[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz
-        		ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz.md5
-        		ln -sf ${externalFastQ_2[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz.md5
-    		fi
- 	fi
+			ln -sf ${externalFastQ_2[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz.md5
+		else
+			ln -sf ${externalFastQ_1[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz
+			ln -sf ${externalFastQ_2[samplenumber]} ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz
+			ln -sf ${externalFastQ_1[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_1.fq.gz.md5
+			ln -sf ${externalFastQ_2[samplenumber]}.md5 ${projectRawTmpDataDir}/${sequencingStartDate[samplenumber]}_${sequencer[samplenumber]}_${run[samplenumber]}_${flowcell[samplenumber]}_L${lane[samplenumber]}_${barcode[samplenumber]}_2.fq.gz.md5
+		fi
+	fi
 done
 
 cd $ROCKETPOINT
 
 echo "before splitting"
 echo $(pwd)
 
-#
-# TODO: array for each sample:
-#
-
 #
 # Create subset of samples for this project.
 #

protocols/CreateInhouseProjects.sh

@@ -1,4 +1,6 @@
 #MOLGENIS walltime=02:00:00 mem=4gb
+set -e
+set -u
 #string tmpName
 #list seqType
 #string project

protocols/GenderCheck.sh

@@ -98,15 +98,15 @@ if ($0 ~ /^#/){
 				print "Unknown"
 			}
 			else if ($2/$1 < 0.65 ){
-				print $2," divided by ",$1," is "$2""/"$1", this is less than 0.65 and this means that the coverage on chromosome X is 0.65 times less than the average coverage of the entire genome, this will most likely be a male";
+				print $2," divided by ",$1," is "$2"/"$1", this is less than 0.65 and this means that the coverage on chromosome X is 0.65 times less than the average coverage of the entire genome, this will most likely be a male";
 				print "Male"
 
 			}else if ($2/$1 > 0.85 ){
-				print $2," divided by ",$1," is "$2""/"$1", this is more than 0.85 and this means that the coverage on chromosome X is almost the same as the average coverage of the entire genome, this will most likely be a female";
+				print $2," divided by ",$1," is "$2"/"$1", this is more than 0.85 and this means that the coverage on chromosome X is almost the same as the average coverage of the entire genome, this will most likely be a female";
 				print "Female"
 			}
 			else{
-				print $2," divided by ",$1," is "$2""/"$1", this is in between the 0.65 and 0.85, we are not sure what the sex is based on the coverage on chromosome X"
+				print $2," divided by ",$1," is "$2"/"$1", this is in between the 0.65 and 0.85, we are not sure what the sex is based on the coverage on chromosome X"
 				print "Unknown" 
 			}
 		}' ${checkSexMeanCoverage} >> ${whichSex}

protocols/QCReport.sh

@@ -108,6 +108,29 @@ cat > ${intermediateDir}/${project}_QCReport.rhtml <<'_EOF'
 <div class="page" style="text-align:center;">
 <font STYLE="font-size: 45pt;">
 <br>
+Preprocessing:
+<br>During the first preprocessing steps of the pipeline, PhiX reads are inserted in each sample to create control SNPs in the dataset.
+<br>Subsequently, Illumina encoding is checked and QC metrics are calculated using a FastQC tool Andrews S. (2010)<sup>1</sup>
+<br>
+<br>Alignment to a reference genome:
+<br>The bwa-mem command from Burrows-Wheeler Aligner(BWA) (Li & Durbin <sup>2</sup>)  is used to align the sequence data to a reference genome resulting in a SAM (Sequence Alignment Map) file.
+<br>The reads in the SAM file are sorted with Sambamba(Tarasov et al.<sup>3</sup>). resulting in a sorted BAM file. When multiple lanes were used during sequencing,
+<br>all lane BAMs were merged into a sample BAM using Sambamba. The (merged) BAM file is marked for duplicates of the same read pair using Sambamba.
+
+<br>Variant discovery:
+<br>The GATK (McKenna et al. <sup>4</sup>) HaplotypeCaller estimates the most likely genotypes and allele frequencies in an alignment using a Bayesian likelihood model
+<br>for every position of the genome regardless of whether a variant was detected at that site or not. This information can later be used in the project based genotyping step.
+<br>A joint analysis has been performed of all the samples in the project. This leads to a posterior probability of a variant allele at a site. 
+<br>SNPs and small Indels are written to a VCF file, along with information such as genotype quality, allele frequency, strand bias and read depth for that SNP/Indel.
+<br>Based on quality thresholds from the GATK "best practices" (Van der Auwera et al.<sup>5</sup>) the SNPs and indels are filtered and marked as Lowqual or Pass resulting in a final VCF file.
+
+<br>Reference
+<br>1. Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data. Available online at:http://www.bioinformatics.babraham.ac.uk/projects/fastqc
+<br>2. Li Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform.
+<br>3. Sambamba: Fast processing of NGS alignment formats 
+<br>4. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data
+<br>5. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline
+
 <br>
 <br>
 <br>

protocols/SnpEff.sh

@@ -1,4 +1,4 @@
-#MOLGENIS walltime=23:59:00 mem=6gb ppn=8
+#MOLGENIS walltime=23:59:00 mem=5gb ppn=2
 
 #Parameter mapping
 #string tmpName
@@ -7,6 +7,7 @@
 #string tempDir
 #string intermediateDir
 #string projectVariantCallsSnpEff_Annotated
+#string projectVariantCallsSnpEff_SummaryHtml
 #string projectBatchGenotypedVariantCalls
 #string project
 #string logsDir 
@@ -29,13 +30,13 @@ then
         ##
 	#
 	#Run snpEff
-	java -XX:ParallelGCThreads=4 -Djava.io.tmpdir=${tempDir} -Xmx4g -jar \
+	java -XX:ParallelGCThreads=2 -Djava.io.tmpdir=${tempDir} -Xmx4g -jar \
 	$EBROOTSNPEFF/snpEff.jar \
 	-v hg19 \
 	-csvStats ${tmpProjectVariantCallsSnpEff_Annotated}.csvStats.csv \
 	-noLog \
 	-lof \
-	-stats ${intermediateDir}/ \
+	-stats ${projectVariantCallsSnpEff_SummaryHtml} \
 	-canon \
 	-ud 0 \
 	-c $EBROOTSNPEFF/snpEff.config \

protocols/VariantCalling.sh

@@ -1,4 +1,4 @@
-#MOLGENIS walltime=23:59:00 mem=14gb ppn=1
+#MOLGENIS walltime=23:59:00 mem=13gb ppn=1
 
 #Parameter mapping
 #string tmpName

test/test_pipeline.sh

@@ -56,7 +56,9 @@ cp generate_template.sh ${workfolder}/generatedscripts/PlatinumSubset/generate_t
 ## Grep used version of molgenis compute out of the parameters file
 fgrep "computeVersion," parameters.csv > ${workfolder}/generatedscripts/PlatinumSubset/mcVersion.txt
 
-perl -pi -e 's|module load NGS_DNA/3.4.1|EBROOTNGS_DNA=/groups/umcg-gaf/tmp04/tmp/NGS_DNA/|' ${workfolder}/generatedscripts/PlatinumSubset/generate_template.sh
+NGS_DNA_VERSION=NGS_DNA/3.4.1
+module load $NGS_DNA_VERSION
+perl -pi -e "s|module load $NGS_DNA_VERSION|EBROOTNGS_DNA=/groups/umcg-gaf/tmp04/tmp/NGS_DNA/|" ${workfolder}/generatedscripts/PlatinumSubset/generate_template.sh
 perl -pi -e 's|PROJECT=projectXX|PROJECT=PlatinumSubset|' ${workfolder}/generatedscripts/PlatinumSubset/generate_template.sh
 perl -pi -e 's|RUNID=runXX|RUNID=run01|' ${workfolder}/generatedscripts/PlatinumSubset/generate_template.sh
 perl -pi -e 's|ngsversion=.*|ngsversion="test";\\|' ${workfolder}/generatedscripts/PlatinumSubset/generate_template.sh
@@ -71,7 +73,7 @@ cd ${workfolder}/generatedscripts/PlatinumSubset/
 sh generate_template.sh 
 
 cd scripts
-perl -pi -e 's|module load \$ngsversion|EBROOTNGS_DNA=/groups/umcg-gaf/tmp04/tmp/NGS_DNA/\nmodule load Molgenis-Compute/\${computeVersion}|' *.sh
+perl -pi -e 's|module load \$ngsversion|EBROOTNGS_DNA=/groups/umcg-gaf/tmp04/tmp/NGS_DNA/\n|' *.sh
 sh submit.sh
 
 cd ${workfolder}/projects/PlatinumSubset/run01/jobs/