We already have sequence reads aligned to a reference genome. How can we extract the final reads that map to a specific gene? To make it clear I want to make a consensus sequence from the reads aligned to each gene. SAMTOOLS suggests the command:

samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq

But it gives many wearied characters after + sign of sequence.

I will thankful if anyone kindly guides me.