Low-frequency variant calling from high-quality mtDNA sequencing data

flexbar -n 16 -r /path/to/reads1.fastq.gz -t output_prefix -f i1.8 -j -z GZ -q 28 -m 50 -a /path/to/adapters.txt -ao 10 -at 1 -ae ANY > file.log.txt

# Create index for the reference genome (required only once)
bwa index -p reference reference.fa
# bwa mem alignment
bwa mem -t 15 -P -T 19 -B 3 -L 5,4 /path/to/reference/reference
/path/to/reads1_flexbar.fastq.gz > reads1.sam
# Sorting and indexing
samtools view -Sbu reads1.sam | samtools sort - -T reads1.sorted -o
reads1.sorted.bam
samtools index reads1.sorted.bam

## Create the split reference genome (required only once) ##
# Create variables
split_genome=reference_split.fa
single_line=reference_singleline.fa
split_data=reference_split.data
# Transform the reference fasta file to contain the sequence as a single line
header=$(cat /path/to/reference.fa | grep '>')
cat /path/to/reference.fa | grep -v '>' | tr -d '\n' > $single_line
# Calculate the length of the input reference genome sequence to determine the cutting position
half_split=$(cat $single_line | awk 'BEGIN {junc=0} junc=int(length($0)/2)
{print junc}')
half_split1=$(echo $half_split | awk '{print $0+1}')
full_len=$(cat $single_line | awk 'BEGIN {len=0} len=length($0) {print len}')
paste <(echo $half_split) <(echo $half_split1) <(echo $full_len) >$split_data
# Take the genome halves according to the calculated positions
gen_start=$(cat $single_line | cut -c1-$half_split)
gen_end=$(cat $single_line | cut -c${half_split1}-$full_len)
# Combine the halves in correct order and restore the fasta format
gen_split=$(paste <(echo $gen_end) <(echo $gen_start) | tr -d '\t')
paste <(echo $header) <(echo $gen_split) | tr '\t' '\n' > $split_genome

# Create index for the reference genome (required only once)
bwa index -p reference_split reference_split.fa
# bwa mem alignment
bwa mem -t 15 -P -T 19 -B 3 -L 5,4 /path/to/reference/reference_split
/path/to/reads1_flexbar.fastq.gz > reads1_junction.sam
# Sorting and indexing
samtools view -Sbu reads1_junction.sam | samtools sort - -T
reads1_junction.sorted -o reads1_junction.sorted.bam
samtools index reads1_junction.sorted.bam

# Quality filter only uniquely aligned reads for further processing (bwa mem mapping quality 0 indicates multi-mapping reads)
# First the reads aligned to normal reference genome
samtools view -bq 1 reads1.sorted.bam > reads1.accepted.bam
samtools index reads1.accepted.bam
# Similar to reads aligned to the split reference genome
samtools view -bq 1 reads1_junction.sorted.bam > reads1_junction.accepted.bam
samtools index reads1_junction.accepted.bam

# Firts calculate coverage for the normal reference alignment
bedtools genomecov -d -ibam reads1.accepted.bam -g reference.fa > coverage.txt
# Then the same for the split reference alignment
bedtools genomecov -d -ibam reads1_junction.accepted.bam -g reference_split.fa > coverage_junction.txt
## Combine the coverages in order to represent the entire mtDNA genome
# Take the middle and end points of the split junction genome for extracting correct lines
mid_point=$(cat reference_split.data | awk '{print $2}')
end_point=$(cat reference_split.data | awk '{print $3}')
# Middle part of the normal coverage file
cat coverage.txt | awk -v endpoint=

# First set the indel qualities for using --call-indels
lofreq indelqual --dindel --ref reference.fa --out reads1.indelqual.bam reads1.accepted.bam
lofreq indelqual --dindel --ref reference_split.fa --out reads1_junction.indelqual.bam reads1_junction.accepted.bam
# Generate the index for both files
samtools index reads1.indelqual.bam
samtools index reads1_junction.indelqual.bam
# Variant calling including indels
lofreq call-parallel --pp-threads 20 -f reference.fa -o reads1.nofilter.vcf -N -B -q 30 -Q 30 --call-indels --no-default-filter reads1.indelqual.bam
lofreq call-parallel --pp-threads 20 -f reference_split.fa -o reads1_junction.nofilter.vcf -N -B -q 30 -Q 30 --call-indels --no-default-filter reads1_junction.indelqual.bam
# Filtering the results if >85% of variant reads are on single strand and for minimum quality
lofreq filter --no-defaults --snvqual-thresh 70 --indelqual-thresh 70 --sb-incl-indels -B 60 -i reads1.nofilter.vcf -o reads1.sbfiltered.vcf
lofreq filter --no-defaults --snvqual-thresh 70 --indelqual-thresh 70 --sb-incl-indels -B 60 -i reads1_junction.nofilter.vcf -o reads1_junction.sbfiltered.vcf

## Combine the original and junction region results for the junction region ##
# Read in the mid and end points of the split genome junction region (generated in Step 3 for the split reference genome)
mid_point=$(cat reference_split.data | awk '{print $2}')
end_point=$(cat reference_split.data | awk '{print $3}')
# Re-coordinate the variants and take only -200 and +200 region
cat reads1_junction.sbfiltered.vcf | awk '/#CHROM/ {flag = 1; next} flag {print}' | \
awk -v midpoint=

## Make sure the snpEff config file uses mitochondrial codons for annotations (required only once)
cd ~
mkdir snpEff_data
# Modify snpEff.config file
# data.dir = ~/snpEff_data/
# Fix mitochondrial codon table usage in the .config file to the wanted genome version:
## Original rows in snpEff.config
## GRCm38.82.genome : Mus_musculus
## GRCm38.82.reference : ftp://ftp.ensembl.org/pub/release-82/gtf/
# Add MT.codonTable in between:
## GRCm38.82.genome : Mus_musculus
## GRCm38.82.MT.codonTable: Vertebrate_Mitochondrial
## GRCm38.82.reference : ftp://ftp.ensembl.org/pub/release-82/gtf/
###

# Filter variants for minimum number of supporting reads in total and on both strands
cat reads1.snvs.vcf | \
java -jar /software/snpEff/4.2/SnpSift.jar filter \