During my analysis I figured out that indicator of sex is also coverage of sex chromosomes. I found out that I can guess sex from two tests. First test is weather the average coverage of chrX is greater than minimum of average coverages of autosomal chromosomes then the sample is probably female. Second test is weather the average coverage of chrY is 10 times smaller than average coverage of autosomal chromosomes then the sample is probably female. Third test is to use +vcf2sex or similar plugin. I am also posting my code to guess male/female samples (19 autosomal chromosomes) written in bash I expect that chromosomes are named (1, X, Y) not (chr1, chrX, chrY).

coveragestats - Script to print filename and coverage statistics from samtools indexed files

echo $*
samtools idxstats $* 2>/dev/null | awk '{ tmp=($3)/($2) ; printf"%s %0.5f\n", $1, tmp }' 2>/dev/null | head -21

autosomalmin - Script to print minimum average coverage of autosomal chromosomes

coveragestats $* | cut -d' ' -f2 | head -20 | tail -19 | awk 'NR == 1 || $1 < min {min = $1}END{printf "%f\n", min}'

xcoverage - Script to print average coverage of chrX

coveragestats $* | grep ^X | cut -d" " -f2

ycoverage - Script to print average coverage of chrY

coveragestats $* | grep ^Y | cut -d" " -f2

maleorfemale - Script to perform first two described sex tests

min=$(autosomalmin $*)
xcov=$(xcoverage $*)
ycov=$(ycoverage $*)
test=$(echo "$xcov"'>'"$min" | bc -l)
if [ $test -eq 1 ]
then
    echo "$*" is female
else
    echo "$*" is male
fi
test=$(echo "$ycov"'<('"$min"/10\) | bc -l)
if [ $test -eq 1 ]
then
        echo "$*" is female
else
        echo "$*" is male
fi

This worked for me and I have pretty low number of false detections.

You want to take the pseudo autosomal regions into account on the X and Y chromosome (https://en.wikipedia.org/wiki/Pseudoautosomal_region). Best would be to check coverage outside of those regions.

You also want to discard reads with low mapping quality.

I once wrote a simple pysam funtion for this to compare the Y chromosome with chr22. You will have to adapt, I can help you if required. It becomes very clear which samples are male and which female. This example is for human data, the coordinates are for hg19. It just works on the (sorted and indexed) bam file.

    yreadcount = 0
    for read in workfile.fetch(region='chrY:2781479-56887902'):
            if not read.mapping_quality < 5 and not read.reference_end == None and not read.reference_start == None:
                    yreadcount += 1
    refreadcount = 0
    for read in workfile.fetch('chr22'):
            if not read.mapping_quality < 5 and not read.reference_end == None and not read.reference_start == None:
                    refreadcount += 1
    print("{} {} ({} vs. {})".format((yreadcount/float(refreadcount)), bam, yreadcount, refreadcount))

1. People generally look at coverage on sex chromosomes or sex-linked contigs. In eutherians, females won't have any coverage on the Y, and males will have approximately half of the coverage on the X relative to the autosomes (assuming WGS).
2. It appears that vcf2sex requires you to specify the ploidy of each region and then uses heterozygous calls in these regions to determine sex.
3. Ploidy is organism-specific. See #2. I would assume that you would specify the regions' ploidy based on the genetic structure (i.e., 2 for autosomes, etc.). I am not aware that this is released formally for genome builds, but it ought to be intuitive.

I wouldn't go through using a program to do this, though. It's somewhat unclear what vcf2sex really does. If you already have the BAM/VCF, I'd try to infer sex using differential coverage of sex chromosomes vs. autosomes.

If the .bam and .bai files are available, one can use e.g., https://github.com/edawson/check-sex

This is similar to what @kulvait has already mentioned!