Getting nucleotide frequencies per position

Hi,

I am very new to bioinformatics and was given the task to analyze some total RNA sequencing results. So far I have converted my .fastaq files to indexed .bam using samtools index. I found the following in literature:

# Generate nucleotide counts per position
samtools index -b ./sorted_bbmerge_norm_$sample.bam

pysamstats -f /Users/kasenriemersma/Bioinformatics/Reference/$reference.fasta -t variation -D 1000000 --format csv --output $outputdir/ntcounts_$sample.csv ./sorted_bbmerge_norm_$sample.bam

I tried running this but I received an error that I was missing the BAM file operand. Is this there any easy way to get past this? Or is there another method that would be better?

Thanks for the help!