Fastq-to-bam

public 1yr ago Version: Version 1 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

Fastq-to-BAM @ NCI-Gadi is a genome alignment workflow that takes raw FASTQ files, aligns them to a reference genome and outputs analysis ready BAM files. This workflow is designed for the National Computational Infrastructure's (NCI) Gadi supercompter, leveraging multiple nodes on NCI Gadi to run all stages of the workflow in parallel, either massively parallel using the scatter-gather approach or parallel by sample. It consists of a number of stages and follows the BROAD Institute's best practice recommendations.

Infrastructure_deployment_metadata: Gadi (NCI)

Code Snippets

set -e

# Wont work if in run_parallel script
postalt=./bwa-postalt.js

ref=../Reference/hs38DH.fasta
outdir=../Align_split
hladir=../HLA_fastq
logdir=Logs/Align
errdir=Logs/Align_error_capture

fqpair=`echo $1 | cut -d ',' -f 1`
fq1=$(ls ${fqpair}*R1.f*q.gz) #Must check regex for each batch
fq2=$(ls ${fqpair}*R2.f*q.gz)
sampleID=`echo $1 | cut -d ',' -f 2`
centre=`echo $1 | cut -d ',' -f 3`
lib=`echo $1 | cut -d ',' -f 4`
platform=`echo $1 | cut -d ',' -f 5`
flowcell=`echo $1 | cut -d ',' -f 6`
lane=`echo $1 | cut -d ',' -f 7`

outPrefix=$(basename $fqpair)
err=${errdir}/${outPrefix}.err

#bwakit emits default sort order (queryname, but no SO tag in headers) or option to sort
#by coordinate with samtools but does not allow the -n flag to specify sort order by name.
#Sambamba requires queryname sorted with SO tag. Below code does it by force.

echo fqpair:$fqpair fq1:$fq1 fq2:$fq2 sampleID:$sampleID centre:$centre lib:$lib platform:$platform flowcell:$flowcell lane:$lane outPrefix:$outPrefix err:$err ref:$ref NCPUS:$NCPUS

seqtk mergepe $fq1 $fq2 | bwa mem -p -t $NCPUS \
        -R "@RG\tID:${flowcell}.${lane}_${sampleID}_${lib}\tPL:${platform}\tPU:${flowcell}.${lane}\tSM:${sampleID}\tLB:${sampleID}_${lib}\tCN:${centre}" \
        -M ${ref} - 2> ${logdir}/${outPrefix}.log.bwamem | k8 ${postalt} -p ${hladir}/${outPrefix} ${ref}.alt \
        | samtools sort -n -@ $NCPUS -o ${outdir}/${outPrefix}.aln.bam

if ! samtools quickcheck ${outdir}/${outPrefix}.aln.bam
then
        printf "Corrupted or missing BAM\n" > $err
fi

Shell SAMtools BWA seqtk From line 24 of fastq-to-bam-v2/align.sh

fastq=`echo $1 | cut -d ',' -f 1`
out=`echo $1 | cut -d ',' -f 2`
logfile=`echo $1 | cut -d ',' -f 3`

echo "$(date): Running fastQC. Fastq:${fastq}, Output:${out}, Log file: ${logfile}, NCPUS:${NCPUS}" >> ${logfile} 2>&1

fastqc -t ${NCPUS} --extract -o ${out} ${fastq} >> ${logfile} 2>&1

Shell FastQC From line 24 of fastq-to-bam-v2/fastqc.sh

fqpair=`echo $1 | cut -d ',' -f 1`
file=$(basename $fqpair)
fq1=$(ls ${fqpair}*R1*f*q.gz) #Must check regex for each batch.
fq2=$(ls ${fqpair}*R2*f*q.gz)
log=./Logs/Split_fastq/${file}.log


fastp -i ${fq1} \
        -I ${fq2} \
        -AGQL \
        -w $NCPUS \
        -S 2000000 \
        -d 0 \
        --out1 ../Split_fastq/${file}_R1.fastq.gz \
        --out2 ../Split_fastq/${file}_R2.fastq.gz 2>${log}