Questions regarding the detection of variants (SNPs and Indels) in mammal Greetings,

I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:

1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3). 

2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).

3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15). 

4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).

5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).

6- Calling variants again using BCFtools.

7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools. 

8- Filtering variants using BCFtools.

9- Select variants SNPs and Indels in two separated files using GATK.

10- Predicting effects using SnpEff.

11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.

Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence? 

In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not? 

Your response is highly appreciated,

Kind regards,
This is close to what many people do. May I add:


between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this. 
filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold.
For that purpose, some time ago I made this and this.