Task 1: Map Trimmed Reads Using BWA MEM

Install BWA

conda create -n bwa -c bioconda bwa
conda activate bwa

Map Reads

bwa mem GCF_018350215.1_P.leo_Ple1_pat1.1_genomic.fna \
sub_samples2/SRR836370_1_subset_val_1.fq.gz \
sub_samples2/SRR836370_2_subset_val_2.fq.gz \
> SRR836370_aligned_reads.sam

Loop for All Samples

for file1 in sub_samples2/*_val_1.fq.gz; do
  file2=${file1/_val_1.fq.gz/_val_2.fq.gz}
  sample_name=$(basename "$file1" _val_1.fq.gz)

  bwa mem GCF_018350215.1_P.leo_Ple1_pat1.1_genomic.fna \
    "$file1" \
    "$file2" \
    > "${sample_name}_aligned_reads.sam"
done

Deactivate BWA Environment

conda deactivate

Task 2: Convert SAM to BAM and Sort

Install and Activate Samtools

conda create -n samtools -c bioconda samtools
conda activate samtools

Convert and Sort

samtools sort SRR836370_aligned_reads.sam \
-o SRR836370_sorted.bam

Batch Conversion and Sorting

for file in *_aligned_reads.sam; do
  sample_name=${file%_aligned_reads.sam}
  samtools sort "$file" -o "${sample_name}_sorted.bam"
done

Deactivate Samtools Environment

conda deactivate

Task 3: Mark Duplicates Using GATK4

Install and Activate GATK4

conda create -n gatk4 -c bioconda gatk4
conda activate gatk4

Add Read Groups

for file in *_sorted.bam
do
  gatk AddOrReplaceReadGroups \
    I="$file" \
    O="${file/_sorted.bam/_sorted_RG.bam}" \
    RGID="${file/_sorted.bam/}" \
    RGPL=ILLUMINA \
    RGPU="${file/_sorted.bam/}" \
    RGSM="${file/_sorted.bam/}" \
    RGLB="${file/_sorted.bam/}"
done

Mark Duplicates

gatk MarkDuplicates \
  -I SRR836370_sorted_RG.bam \
  -O SRR836370_deduplicated.bam \
  -M SRR836370_duplication_metrics.txt \
  --REMOVE_DUPLICATES true

Loop Version

for file in *_sorted_RG.bam; do
  base=${file%_sorted_RG.bam}

  gatk MarkDuplicates \
    -I "$file" \
    -O "${base}_deduplicated.bam" \
    -M "${base}_duplication_metrics.txt" \
    --REMOVE_DUPLICATES true
done

Task 4: Index Deduplicated Files

Index a single BAM file

samtools index SRR836370_deduplicated.bam

Batch Indexing

for file in *_deduplicated.bam; do
  samtools index "$file"
done

Task 5: Estimate Coverage Using Qualimap

Install and Activate Qualimap

conda create -n qualimap -c bioconda qualimap
conda activate qualimap

Run on a Single File

qualimap bamqc \
  -bam SRR836370_deduplicated.bam \
  -outdir SRR836370_qualimap_results \
  -outformat HTML

Run on All Files

for file in *_deduplicated.bam; do
  sample_name=${file%_deduplicated.bam}

  qualimap bamqc \
    -bam "$file" \
    -outdir "${sample_name}_qualimap_results" \
    -outformat HTML
done

Check Output

cd SRR836370_qualimap_results
cat genome_results.txt

Extra Note

The reference genome, FASTQ files, metadata, and variant files are available from the Wildlife Genomics 4 Zenodo dataset.

Zenodo dataset: https://zenodo.org/records/21871083