RNAseq SNP discovery: deciding upon filters and dealing with allele expression bias

Question

I am working with non-model plant RNA samples which we have been deep sequenced and analysed using STAR aligner under default parameters.

Aim We would like to conduct SNP discovery of these samples.

Objective Our ultimate goal with this genotypic data is to search for variants (both SNPs and indels) that fall within genes from a cohort of gene families, which may be responsible for phenotypic variation that we have assessed among our sample plants.

Background The RNA alignments use a reference genome of the species I am dealing with, and following GATK best practices and produced a VCF file for each of my samples.

There are three things that give me pause,

1. Am I correct in merging the .VCF files after calling SNPs, rather than merging .BAM files and then calling SNPs?
2. What filters are required for MAF, read-depth, base call quality, missing data, snp proximity to each other and proximity to indels, etc.

• Should I produce plots to make these decisions?

3. Lastly how do I deal with specific allele expression skewing the results, my assumption so far is that other studies seem to have accounted for this by merely using a low bar for the MAF filter (such as >= 0.1), thus ensuring that as many true SNPs as possible are captured.

Any insight into any of these questions, or a point in the direction of where I might find it, would be much appreciated.