Willy Rodríguez
Menu
CV

Tutorial

Applying PSMC to Neanderthal (Altai) data

Published 2015-07-08

In 2014, Prüfer and colleagues published The complete genome sequence of a Neanderthal from the Altai Mountains in Nature. As part of that study, the authors carried out a PSMC analysis. I describe here how I was able to reproduce the analysis on 3 July 2015.

The tools I used were:

Download and compile the tools

git clone https://github.com/samtools/samtools.git
git clone https://github.com/samtools/htslib.git
git clone https://github.com/samtools/bcftools.git

cd bcftools
make
cd ../samtools
make

The read sequences of the Altai Neanderthal were aligned to the GRCh37 build of the human reference. The 1000 Genomes reference used here remains available as human_g1k_v37.fasta.gz (the 2015 tutorial used the equivalent FTP path).

Assume we have all the aligned chromosomes in a folder called bam, named chr1.bam, chr2.bam, and so on.

Uncompress and index the reference file:

gunzip ./human_g1k_v37.fasta.gz
./samtools/samtools faidx human_g1k_v37.fasta

This will produce a file human_g1k_v37.fasta.fai.

Variant calling pipeline

./samtools/samtools mpileup -C50 -uf ./human_g1k_v37.fasta ./bam/chr1.bam | \
./bcftools/bcftools call -c | ./bcftools/vcfutils.pl vcf2fq -d 10 -D 100 | \
gzip > ./chr1.fq.gz

Samtools can also open a BAM file on a remote FTP or HTTP server. In 2015 the Altai BAM files were hosted at the Max Planck Institute EVA; the file used in the original tutorial remains available at AltaiNea.hg19_1000g.1.dq.bam.

The .fq file produced by vcfutils.pl contains lowercase and uppercase letters (for example nnATgagtttAGnn), which is confusing given that a FASTA sequence is often expected to contain only uppercase characters. After some searching, lowercase appears to indicate low coverage.

Convert to PSMC input

To obtain the input sequence for PSMC, I used the fq2psmcfa tool that comes with PSMC. For each chromosome:

for i in {1..22};
  do ../psmc/utils/fq2psmcfa -q20 ./chr$i.fq.gz > chr$i.psmcfa;
done

Concatenate the chromosomes:

for i in {1..22};
  do cat ./chr$i.psmcfa >> AltaiNea.psmcfa;
done

Run PSMC

I used the default parameters. Assuming you have a compiled version of PSMC inside a folder called psmc and your data (AltaiNea.psmcfa) is in the current directory:

./psmc/psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o ./AltaiNea.psmc ./AltaiNea.psmcfa

This will produce the inferred history (it can take a while). The output is written to AltaiNea.psmc.

To plot the inferred demographic history:

./psmc/utils/psmc_plot.pl -p ./AltaiNea.pdf ./AltaiNea.psmc

This produces a PDF plot of the inferred history. Later work on the IICR discusses how such plots should be interpreted when the population is structured.