04_ONT_alignment.md

ONT long-read alignment with Winnowmap

Once long-reads were trimmed and filtered, they were aligned to the kākāpō reference assembly (Jane's genome) for SV discovery. As with the Illumina data, reads were aligned with or without the W sex chromosome depending on the bird's sex.

trim=/kakapo-data/ONT/trimmed/
align=/kakpo-data/ONT/winnowmap/alignments/
mref=/kakapo-data/References/kakapo_no_Wchromosome.fa
fref=/kakapo-data/References/kakapo_full_ref.fa
source ~/anaconda3/etc/profile.d/conda.sh

ONT read alignment

Reads were aligned to Jane's genome using Winnowmap v2.03.

cd ${data}winnowmap

meryl count k=15 ${fref} output kakapo_full.meryl
meryl count k=15 ${mref} output kakapo_noW.meryl

meryl print greater-than distinct=0.9998 kakapo_full.meryl > kakapo_full_repetitive_k15.txt
meryl print greater-than distinct=0.9998 kakapo_noW.meryl > kakapo_noW_repetitive_k15.txt

Alignments were run separately for males and females. To save some time, loops were set to run through each trimmed length simultaneously for each individual.

for male in Bill Blades Gulliver Rangi Sass Smoko
    do
    for length in {3,4,5,10}
        do
        winnowmap -W kakapo_noW_repetitive_k15.txt -ax map-ont ${mref} ${trim}${male}_q10_${length}kbtrim.fastq.gz > ${align}sam/${male}_${length}kb.sam
        echo "Converting sam to bam for $base..."
        samtools view -@64 -T ${mref} -b ${align}sam/${male}_${length}kb.sam | samtools sort -@ 64 -o ${align}bam/${male}_${length}kb.bam
        samtools index -@32 ${align}bam/${male}_${length}kb.bam
    done
done
for female in Bella Kuia Margaret-Maree Sue
    do
    for length in {3,4,5,10}
        do
        winnowmap -W kakapo_full_repetitive_k15.txt -ax map-ont ${fref} ${trim}${female}_q10_${length}kbtrim.fastq.gz > ${align}sam/${female}_${length}kb.sam
        echo "Converting sam to bam for $base..."
        samtools view -@64 -T ${fref} -b ${align}sam/${female}_${length}kb.sam | samtools sort -@ 64 -o ${align}bam/${female}_${length}kb.bam
        samtools index -@32 ${align}bam/${female}_${length}kb.bam
    done
done

Fun way to implement an if/then loop. Wasn't actually used.

for sam in ${data}winnowmap/alignment/sam/*.sam
    do
    base=$(basename $sam .sam)
    if [[ $base == *(Bella|Kuia|Margaret-Maree|Sue)* ]]
        then
            echo "Converting sam to bam for $base..."
            samtools view -@64 -T ${fref} -b $sam | samtools sort -@ 64 -o ${data}winnowmap/alignment/bam/${base}.bam
        else
            echo "Converting sam to bam for $base..."
            samtools view -@64 -T ${mref} -b $sam | samtools sort -@ 64 -o ${data}winnowmap/alignment/bam/${base}.bam
    fi
done

Read Mapping Quality

Depths were estimated using Mosdepth v0.3.2 and Qualimap v2.2.2. Summaries were visualised using MultiQC v1.9.

For future inclusion, samples had to have a coverage of >=4x.

for bam in ${align}bam/*.bam
    do
    indiv=$(basename $bam .bam)
    echo "Running Qualimap for ${indiv}..."
    qualimap bamqc \
        -bam ${bam} \
        -nw 10000 \
        -nt 16 -c \
        -outdir ${align}stats/${indiv}.graphmap \
        --java-mem-size=8G &
    echo "Running mosdepth for $indiv"
    mosdepth --threads 24 --fast-mode --by 50 ${align}stats/${indiv} ${bam}
    wait
done
multiqc ${align}stats

After review of the plots generated by mosdepth, A and K fell below the minimum depth thresholds and were removed from future analyses.