This software package infers population size history from a diploid sequence
using the Pairwise Sequentially Markovian Coalescent (PSMC) model. The
detailed model is described in file `psmc.tex'.

To compile the binaries, you may run

    make; (cd utils; make)

After that, you may try

    utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa
    psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa
    utils/psmc2history.pl diploid.psmc | utils/history2ms.pl > ms-cmd.sh
    utils/psmc_plot.pl diploid diploid.psmc

where `diploid.fq.gz' is typically the whole-genome diploid consensus sequence
of one human individual, which can be generated by, for example:

    samtools mpileup -C50 -uf ref.fa aln.bam | bcftools view -c - \
      | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > diploid.fq.gz

Here option -d sets and minimum read depth and -D sets the maximum. It is
recommended to set -d to a third of the average depth and -D to twice.  Program
`fq2psmcfa' transforms the consensus sequence into a fasta-like format where
the i-th character in the output sequence indicates whether there is at least
one heterozygote in the bin [100i, 100i+100).

For dipcall output, you may use the following to generate psmcfa:

    seqtk mutfa ref.fa <(gzip -dc prefix.dip.vcf.gz|utils/vcf2snp.pl -) \
      | seqtk seq -cM prefix.dip.bed -l80 | utils/fq2psmcfa - > prefix.psmcfa

Program `psmc' infers the population size history. In particular, the `-p'
option specifies that there are 64 atomic time intervals and 28 (=1+25+1+1)
free interval parameters. The first parameter spans the first 4 atomic time
intervals, each of the next 25 parameters spans 2 intervals, the 27th spans 4
intervals and the last parameter spans the last 6 time intervals. The `-p' and
`-t' options are manually chosen such that after 20 rounds of iterations, at
least ~10 recombinations are inferred to occur in the intervals each parameter
spans. Impropriate settings may lead to overfitting. The command line in the
example above has been shown to be suitable for modern humans.

The `psmc' program infers the scaled mutation rate, the recombination rate and
the free population size parameters. All these parameters are scaled to 2N0. You
may run `psmc2history.pl' combined with `history2ms.pl' to generate the ms
command line that simulates the history inferred by PSMC, or visualize the result
with `psmc_plot.pl'.

To perform bootstrapping, one has to run splitfa first to split long chromosome
sequences to shorter segments. When the `-b' option is applied, psmc will then
randomly sample with replacement from these segments. As an example, the
following command lines perform 100 rounds of bootstrapping:

    utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa
	utils/splitfa diploid.psmcfa > split.psmcfa
    psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa
	seq 100 | xargs -i echo psmc -N25 -t15 -r5 -b -p "4+25*2+4+6" \
	    -o round-{}.psmc split.fa | sh
    cat diploid.psmc round-*.psmc > combined.psmc
	utils/psmc_plot.pl -pY50000 combined combined.psmc

One probably wants to modify the "xargs" command-line to parallelize PSMC.

If you have questions about PSMC, please ask at <https://hengli.uservoice.com/>.
You do not need to register unless you also want to modify your own questions.
You may also post comments at github (if you have a github account). I want to
make the question and the answer public such that others can see them and I do
not need to answer the same question multiple times. Thank you for using PSMC.


APPENDIX I: Scaling the PSMC output
===================================

The PSMC output is scaled to the 2N_0. There are two ways of rescaling the time
and the popuation size more meaningfully.

Firstly, suppose we know the per-site per-generation mutation rate \mu, we can
compute N_0 as:

  N_0 = \theta_0 / (4\mu) / s

where \theta_0 is given at the 2nd column of "TR" lines, and s is the bin size
we use for generating the PSMC input. Knowing N_0, we can scale time to
generations and relative population size to effective size by

  T_k = 2N_0 * t_k
  N_k = N_0 * \lambda_k

where t_k and \lambda_k are given at the 3rd and 4th columns of "RS" lines,
respectively.

A problem with the above strategy is that we do not know a definite answer of
\mu and in fact it various with regions and mutation types. An alternative way
is to use per-site pairwise sequence divergence to represent time:

  d_k = 2\mu * T_k = t_k * \theta_0 / s

and use scaled mutation rate to represent population size:

  \theta_k = 4N_k * \mu = \lambda_k * \theta_0 / s

where, again, t_k and \lambda_k are given at the "RS" line, \theta_0 at the
"TR" line and s is the bin size, which defaults to 100 in fq2psmcfa.


APPENDIX II: Correcting for low coverage
========================================

For diploid genomes sequenced to low coverage, heterozygotes will be randomly
lost due to the lack of coverage of both alleles. This has the same effect as
smaller mutation rate and can be corrected. If you know the fraction of hets
missed due to low coverage, you can generate the PSMC plot with:

  psmc_plot.pl -M "sample1=0.1,sample2=0.2" prefix sample1.psmc sample2.psmc

This says that sample1 has 10% false negative rate (FNR) on hets and sample2
has 20%. The plotting script does not correct FNR for bootstrapping. If you
want to plot the result with your own scripts, you can increase \theta_0 to
\theta_0/(1-FNR).

Unfortunately, I haven't found a reliable way to estimate the background FNR
relevant to PSMC. The simple and unscientific approach is to align the PSMC
curves by eye. Probably a better solution is to downsample a high-coverage
sample to a certain coverage and measures FNR. I have not done this.

Not correcting for low coverage is the most common pitfall when using PSMC.