Nextflow workflow report

erne-create --methyl-hash --output-prefix genome --fasta genome.fa --fasta lambda.fa

mkdir fastq fastq/logs
        cutadapt -j 2 -a AGATCGGAAGAGC -A AGATCGGAAGAGC \
        -q 20 -m 36 -O 3 \
        -o fastq/merge.sampleA_1.fastq.gz \
        -p fastq/merge.sampleA_2.fastq.gz sampleA_1.fastq.gz sampleA_2.fastq.gz \
        > fastq/logs/cutadapt.sampleA.merge.log 2>&1

mkdir fastq fastq/logs
        cutadapt -j 2 -a AGATCGGAAGAGC -A AGATCGGAAGAGC \
        -q 20 -m 36 -O 3 \
        -o fastq/input.sampleB_1.fastq.gz \
        -p fastq/input.sampleB_2.fastq.gz sampleB_1.fastq.gz sampleB_2.fastq.gz \
        > fastq/logs/cutadapt.sampleB.input.log 2>&1

mkdir fastq
        echo -e "1\n2" | xargs -I{} sh -c 'cat *$1 > fastq/sampleB_$1.fastq.gz' -- {}

mkdir fastq fastq/fastqc fastq/logs
    fastqc sampleB_1.fastq.gz sampleB_2.fastq.gz -threads 2 -outdir=fastq/fastqc > fastq/logs/fastqc.sampleB.log

mkdir sampleB sampleB/bam sampleB/bam/logs
        segemehl.x -i genome.ctidx -j genome.gaidx \
        -d genome.fa lambda.fa \
        -q merge.sampleB_1.fastq.gz -p merge.sampleB_2.fastq.gz -o raw.segemehl.sam -I 500 -A 90 -s -t 8 -F 1 -H 1 -D 1 \
        > sampleB/bam/logs/raw.segemehl.log 2>&1 || exit $?

        grep -v "^@RG" raw.segemehl.sam | samtools addreplacerg -r 'ID:sampleB' -r 'SM:sampleB' - |
        samtools view -Sb - > sampleB/bam/raw.segemehl.bam

mkdir sampleA sampleA/bam sampleA/bam/logs

        erne-bs5 --reference genome.ebm --query1 input.sampleA_1.fastq.gz --query2 input.sampleA_2.fastq.gz --fragment-size-min 0 --fragment-size-max 500 \
        --errors -1 --threads 6 --output unsorted.erne-bs5.bam --print-all \
        > sampleA/bam/logs/raw.erne-bs5.log 2>&1 || exit $?

        samtools view -h unsorted.erne-bs5.bam | grep -v "^@RG" | samtools addreplacerg -r 'ID:sampleA' -r 'SM:sampleA' - |
        samtools sort -T deleteme -m 241591910 -@ 8 \
        -o sampleA/bam/raw.erne-bs5.bam -

mkdir sampleA sampleA/bam
        correct_sam_tlens -i raw.segemehl.bam > tlens.segemehl.bam || exit $?
        samtools sort -T deleteme -m 966367642 -@ 2 \
        -o sorted.segemehl.bam tlens.segemehl.bam
        samtools index sorted.segemehl.bam
        
        filter_sam_xf_tag.py -c 0.03 sorted.segemehl.bam sampleA/bam/proc.segemehl.bam || exit $?
        ln -s sampleA/bam/proc.segemehl.bam sampleA.bam

mkdir sampleA sampleA/bam
        ln -s sampleA/bam/proc.erne-bs5.bam sampleA.bam

        if [ -z "$(tail -c 1 genome.fa)" ]
        then
        cat genome.fa lambda.fa > fasta.tmp  
        else
        cat genome.fa > fasta.tmp; echo >> fasta.tmp; cat lambda.fa >> fasta.tmp      
        fi

        samtools index raw.erne-bs5.bam
        correct_sam_format.py -i 500 -t 4 -T . raw.erne-bs5.bam corrected.erne-bs5.bam || exit $?

        samtools sort -T deleteme -m 3382286746 -@ 4 \
        -no unsorted.erne-bs5.bam corrected.erne-bs5.bam
        correct_sam_tlens -i unsorted.erne-bs5.bam > tlens.erne-bs5.bam || exit $?

        samtools sort -T deleteme -m 3382286746 -@ 4 \
        -o sorted.erne-bs5.bam tlens.erne-bs5.bam
        samtools index sorted.erne-bs5.bam

        filter_sam_erne.py -c 0.03 -t 4 -T . fasta.tmp sorted.erne-bs5.bam sampleA/bam/proc.erne-bs5.bam

mkdir sampleA sampleA/bam
    ln -s sampleA/bam/merged.bam sampleA.bam
    samtools merge -c sampleA/bam/merged.bam proc.erne-bs5.bam proc.segemehl.bam

mkdir sampleA sampleA/bam
    ln -s sampleA/bam/subset.bam sampleA.bam

    # reset alignments whose mates have erroneously aligned to lambda
    samtools view -h merged.bam |
    awk 'BEGIN {OFS="\t"} {if(($3=="J02459.1" && $7=="=") || $7=="J02459.1")     {if($2==97){$2=73} else { if($2==81){$2=89}     else { if($2==161){$2=137} else{ if($2==145){$2=153} }}}; $7="*"; $8=0}; {print $0}}' |
    
    samtools sort -T deleteme -m 966367642 -@ 2 \
    -o sort.bam - || exit $?
    samtools index sort.bam

    # remove lambda from the header and split bams
    samtools view -H sort.bam | awk '$1~/^@/ && $2!="SN:J02459.1"' > header.sam || exit $?
    awk '{printf("%s\t0\t%s\n",$1,$2);}' genome.fa.fai > fai.tmp
    samtools view -bL fai.tmp sort.bam | samtools reheader header.sam - > sampleA/bam/subset.bam || exit $?
    samtools view sort.bam J02459.1 | samtools view -bt lambda.fa.fai - > sampleA/bam/J02459.1.bam

mkdir sampleA sampleA/bam

    filter_sam_uniqs.py J02459.1.bam sampleA/bam/lambda.unique.bam sampleA/bam/lambda.multimapped.bam
    ln -s sampleA/bam/unique.bam sampleA.bam

mkdir sampleA sampleA/bam

    filter_sam_uniqs.py subset.bam sampleA/bam/unique.bam sampleA/bam/multimapped.bam
    ln -s sampleA/bam/unique.bam sampleA.bam