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,