Ancient DNA exercise - Revision history

Mick at 10:18, 9 January 2026

2026-01-09T10:18:44Z

← Older revision		Revision as of 12:18, 9 January 2026
Line 37:		Line 37:

	<pre>		<pre>
	/home/ctools/apps/mapdamage2/2.2.3/venv/bin~~/mapDamage~~/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/modern/NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.bam -r /home/databases/references/human/hs37d5.fa --no-stats		/home/ctools/apps/mapdamage2/2.2.3/venv/bin/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/modern/NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.bam -r /home/databases/references/human/hs37d5.fa --no-stats
	</pre>		</pre>

Line 56:		Line 56:
	Run mapDamage on one of the ancient individuals		Run mapDamage on one of the ancient individuals
	<pre>		<pre>
	/home/ctools/apps/mapdamage2/2.2.3/venv/bin~~/mapDamage~~/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ancient/allentoft_2015/RISE559.sort.rmdup.realign.md.100k.bam -r /home/databases/references/human/hs37d5.fa --no-stats		/home/ctools/apps/mapdamage2/2.2.3/venv/bin/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ancient/allentoft_2015/RISE559.sort.rmdup.realign.md.100k.bam -r /home/databases/references/human/hs37d5.fa --no-stats
	</pre>		</pre>

Mick at 10:07, 9 January 2026

2026-01-09T10:07:16Z

← Older revision		Revision as of 12:07, 9 January 2026
Line 37:		Line 37:

	<pre>		<pre>
	mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/modern/NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.bam -r /home/databases/references/human/hs37d5.fa --no-stats		/home/ctools/apps/mapdamage2/2.2.3/venv/bin/mapDamage/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/modern/NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.bam -r /home/databases/references/human/hs37d5.fa --no-stats
	</pre>		</pre>

Line 56:		Line 56:
	Run mapDamage on one of the ancient individuals		Run mapDamage on one of the ancient individuals
	<pre>		<pre>
	mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ancient/allentoft_2015/RISE559.sort.rmdup.realign.md.100k.bam -r /home/databases/references/human/hs37d5.fa --no-stats		/home/ctools/apps/mapdamage2/2.2.3/venv/bin/mapDamage/mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ancient/allentoft_2015/RISE559.sort.rmdup.realign.md.100k.bam -r /home/databases/references/human/hs37d5.fa --no-stats
	</pre>		</pre>

Line 62:		Line 62:
	<pre>		<pre>
	cd RISE559.sort.rmdup.realign.md.100k.mapDamage/		cd RISE559.sort.rmdup.realign.md.100k.mapDamage/
	Length_plot.pdf &		Length_plot.pdf
	Fragmisincorporation_plot.pdf &		Fragmisincorporation_plot.pdf
	</pre>		</pre>

Mick at 13:38, 6 January 2026

2026-01-06T13:38:57Z

← Older revision		Revision as of 15:38, 6 January 2026
Line 256:		Line 256:
	We can plot the first two principal components using the custom R script plotPca.R		We can plot the first two principal components using the custom R script plotPca.R

	The three positional arguments are the eigenvector file, sample info file and prefix for the output		The three positional arguments are the eigenvector file, sample info file and prefix for the output <b>(view the pdf either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R modern.eigenvec world.sampleInfo.txt modern		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R modern.eigenvec world.sampleInfo.txt modern
	~~evince~~ modern.pca.plot.pdf &		modern.pca.plot.pdf
	</pre>		</pre>

Line 266:		Line 266:
	'''Q14:''' which populations are most differentiated along PC2?		'''Q14:''' which populations are most differentiated along PC2?

	We repeat the exercise on a subset of European populations:		We repeat the exercise on a subset of European populations <b>(view the pdf either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
	/home/ctools/plink --bfile RISE507.world --keep-clusters eur.poplist --within RISE507.world.cluster --pca header tabs --out eur		/home/ctools/plink --bfile RISE507.world --keep-clusters eur.poplist --within RISE507.world.cluster --pca header tabs --out eur
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R eur.eigenvec world.sampleInfo.txt eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R eur.eigenvec world.sampleInfo.txt eur
	~~evince~~ eur.pca.plot.pdf &		eur.pca.plot.pdf
	</pre>		</pre>

Line 278:		Line 278:


	Now, let us examine how the cluster of ancient individuals compared to the modern ones:		Now, let us examine how the cluster of ancient individuals compared to the modern ones <b>(view the pdf either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
	/home/ctools/plink --bfile RISE507.world --pca header tabs --out ancient.world		/home/ctools/plink --bfile RISE507.world --pca header tabs --out ancient.world
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.world.eigenvec world.sampleInfo.txt ancient.world		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.world.eigenvec world.sampleInfo.txt ancient.world
	~~evince~~ ancient.world.pca.plot.pdf &		ancient.world.pca.plot.pdf
	</pre>		</pre>

Line 326:		Line 326:
	'''Q17:''' which ancient individuals don't cluster close to any modern individuals? what could be a plausible reason?		'''Q17:''' which ancient individuals don't cluster close to any modern individuals? what could be a plausible reason?

	Repeat the exercise but remove the non-European modern individuals:		Repeat the exercise but remove the non-European modern individuals <b>(view the pdf either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
	/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --remove-clusters noneur.poplist --pca header tabs --out ancient.eur		/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --remove-clusters noneur.poplist --pca header tabs --out ancient.eur
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.eur.eigenvec world.sampleInfo.txt ancient.eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.eur.eigenvec world.sampleInfo.txt ancient.eur
	~~evince~~ ancient.eur.pca.plot.pdf &		ancient.eur.pca.plot.pdf
	</pre>		</pre>

	'''Q18:''' which populations are most differentiated along PC1? what could be a plausible reason?		'''Q18:''' which populations are most differentiated along PC1? what could be a plausible reason?

	As a final exercise, we now project the ancient individual on PCs inferred from modern Europeans:		As a final exercise, we now project the ancient individual on PCs inferred from modern Europeans<b>(view the pdf either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
	/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --pca-clusters eur.poplist --remove-clusters noneur.poplist --pca header tabs --out ancient_proj.eur --maf 0.01		/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --pca-clusters eur.poplist --remove-clusters noneur.poplist --pca header tabs --out ancient_proj.eur --maf 0.01
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient_proj.eur.eigenvec world.sampleInfo.txt ancient_proj.eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient_proj.eur.eigenvec world.sampleInfo.txt ancient_proj.eur
	~~evince~~ ancient_proj.eur.pca.plot.pdf &		ancient_proj.eur.pca.plot.pdf
	</pre>		</pre>

Mick at 13:34, 6 January 2026

2026-01-06T13:34:33Z

← Older revision		Revision as of 15:34, 6 January 2026
Line 211:		Line 211:
	Now we will randomly sample a DNA fragment at each position and output the results in VCF format (custom python script):		Now we will randomly sample a DNA fragment at each position and output the results in VCF format (custom python script):
	<pre>		<pre>
	zcat RISE507.pileup.gz \| python2 /home/projects/22126_NGS/exercises/adna/02_popgen/get_haploid_vcf_from_pileup.py -r -s RISE507 \|bgzip -c > RISE507.vcf.gz		zcat RISE507.pileup.gz \| /home/ctools/Python-2.7.18/bin/python2.7 /home/projects/22126_NGS/exercises/adna/02_popgen/get_haploid_vcf_from_pileup.py -r -s RISE507 \|/home/ctools/htslib-1.20/bgzip -c > RISE507.vcf.gz
	</pre>		</pre>
	This is done because the coverage is insufficient to ensure proper genotyping.		This is done because the coverage is insufficient to ensure proper genotyping.

Mick at 13:13, 6 January 2026

2026-01-06T13:13:30Z

← Older revision		Revision as of 15:13, 6 January 2026
Line 40:		Line 40:
	</pre>		</pre>

	Examine the output (either via mobaxterm or by downloading it locally):		Examine the output <b>(either via mobaxterm or by downloading it locally):</b>

	<pre>		<pre>
Line 59:		Line 59:
	</pre>		</pre>

	Examine the output (either via mobaxterm or by downloading it locally)		Examine the output <b>(either via mobaxterm or by downloading it locally)</b>
	<pre>		<pre>
	cd RISE559.sort.rmdup.realign.md.100k.mapDamage/		cd RISE559.sort.rmdup.realign.md.100k.mapDamage/

Mick at 13:12, 6 January 2026

2026-01-06T13:12:10Z

← Older revision		Revision as of 15:12, 6 January 2026
Line 40:		Line 40:
	</pre>		</pre>

	Examine the output:		Examine the output (either via mobaxterm or by downloading it locally):

	<pre>		<pre>
	cd NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.mapDamage/		cd NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.mapDamage/
	~~okular~~ Length_plot.pdf &		Length_plot.pdf
	~~okular~~ Fragmisincorporation_plot.pdf &		Fragmisincorporation_plot.pdf
	cd ..		cd ..
	</pre>		</pre>
Line 59:		Line 59:
	</pre>		</pre>

	Examine the output		Examine the output (either via mobaxterm or by downloading it locally)
	<pre>		<pre>
	cd RISE559.sort.rmdup.realign.md.100k.mapDamage/		cd RISE559.sort.rmdup.realign.md.100k.mapDamage/
	~~okular~~ Length_plot.pdf &		Length_plot.pdf &
	~~okular~~ Fragmisincorporation_plot.pdf &		Fragmisincorporation_plot.pdf &
	</pre>		</pre>

Mick at 13:08, 6 January 2026

2026-01-06T13:08:41Z

← Older revision		Revision as of 15:08, 6 January 2026
Line 158:		Line 158:

	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile world --missing --out world		/home/ctools/plink --bfile world --missing --out world
	</pre>		</pre>

Line 222:		Line 222:
	We convert to plink binary format:		We convert to plink binary format:
	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --vcf RISE507.vcf.gz --make-bed --double-id --out RISE507		/home/ctools/plink --vcf RISE507.vcf.gz --make-bed --double-id --out RISE507
	</pre>		</pre>

	Try to merge the sample with the reference panel		Try to merge the sample with the reference panel
	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile world --bmerge RISE507 --out RISE507.merge		/home/ctools/plink --bfile world --bmerge RISE507 --out RISE507.merge
	</pre>		</pre>

Line 236:		Line 236:
	We will remove the failing SNPs and try again		We will remove the failing SNPs and try again
	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507 --exclude RISE507.merge.missnp --make-bed --out RISE507.merge2		/home/ctools/plink --bfile RISE507 --exclude RISE507.merge.missnp --make-bed --out RISE507.merge2
	/home/ctools~~/plink-1.9~~/plink --bfile world --bmerge RISE507.merge2 --out RISE507.world		/home/ctools/plink --bfile world --bmerge RISE507.merge2 --out RISE507.world
	</pre>		</pre>

Line 251:		Line 251:
	First, we will have a look at the modern populations in the reference panel:		First, we will have a look at the modern populations in the reference panel:
	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507.world --keep-clusters modern.poplist --within RISE507.world.cluster --pca header tabs --out modern		/home/ctools/plink --bfile RISE507.world --keep-clusters modern.poplist --within RISE507.world.cluster --pca header tabs --out modern
	</pre>		</pre>

Line 269:		Line 269:

	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507.world --keep-clusters eur.poplist --within RISE507.world.cluster --pca header tabs --out eur		/home/ctools/plink --bfile RISE507.world --keep-clusters eur.poplist --within RISE507.world.cluster --pca header tabs --out eur
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R eur.eigenvec world.sampleInfo.txt eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R eur.eigenvec world.sampleInfo.txt eur
	evince eur.pca.plot.pdf &		evince eur.pca.plot.pdf &
Line 281:		Line 281:

	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507.world --pca header tabs --out ancient.world		/home/ctools/plink --bfile RISE507.world --pca header tabs --out ancient.world
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.world.eigenvec world.sampleInfo.txt ancient.world		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.world.eigenvec world.sampleInfo.txt ancient.world
	evince ancient.world.pca.plot.pdf &		evince ancient.world.pca.plot.pdf &
Line 329:		Line 329:

	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507.world --within RISE507.world.cluster --remove-clusters noneur.poplist --pca header tabs --out ancient.eur		/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --remove-clusters noneur.poplist --pca header tabs --out ancient.eur
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.eur.eigenvec world.sampleInfo.txt ancient.eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.eur.eigenvec world.sampleInfo.txt ancient.eur
	evince ancient.eur.pca.plot.pdf &		evince ancient.eur.pca.plot.pdf &
Line 339:		Line 339:

	<pre>		<pre>
	/home/ctools~~/plink-1.9~~/plink --bfile RISE507.world --within RISE507.world.cluster --pca-clusters eur.poplist --remove-clusters noneur.poplist --pca header tabs --out ancient_proj.eur --maf 0.01		/home/ctools/plink --bfile RISE507.world --within RISE507.world.cluster --pca-clusters eur.poplist --remove-clusters noneur.poplist --pca header tabs --out ancient_proj.eur --maf 0.01
	Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient_proj.eur.eigenvec world.sampleInfo.txt ancient_proj.eur		Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient_proj.eur.eigenvec world.sampleInfo.txt ancient_proj.eur
	evince ancient_proj.eur.pca.plot.pdf &		evince ancient_proj.eur.pca.plot.pdf &

WikiSysop: Created page with "

Overview

Adapted from Martin Sikora. First:

Navigate to your home directory:
Create a directory called "adna"
Navigate to the directory you just created.

We will try to # Authenticate ancient DNA # do some basic population genetics

Data authentication

Authentication involves making sure that the DNA that you have extracted from my fossil and sequenced is indeed from the fossil and not some modern contaminant. A big differe..."

2024-03-19T15:35:43Z

Created page with "<H2>Overview</H2> Adapted from Martin Sikora. First: <OL> <LI>Navigate to your home directory: <LI>Create a directory called "adna" <LI>Navigate to the directory you just created. </OL> We will try to # Authenticate ancient DNA # do some basic population genetics <h2> Data authentication</h2> Authentication involves making sure that the DNA that you have extracted from my fossil and sequenced is indeed from the fossil and not some modern contaminant. A big differe..."

New page

<H2>Overview</H2>

Adapted from Martin Sikora.

First:
<OL>
<LI>Navigate to your home directory:
<LI>Create a directory called "adna"
<LI>Navigate to the directory you just created.
</OL>

We will try to
# Authenticate ancient DNA
# do some basic population genetics

<h2> Data authentication</h2>

Authentication involves making sure that the DNA that you have extracted from my fossil and sequenced is indeed from the fossil and not some modern contaminant. A big difference between modern DNA and ancient DNA is the presence of chemical damage due to the passage of time.

<h3> Direct measurements of the rate of chemical damage</h3>

First, create a directory:
<pre>
mkdir 01_authentication
cd 01_authentication
</pre>

We will characterize DNA damage patterns using mapDamage, a software to estimate the rate of nucleotide substitution. In this section, we will examine some example BAM files for the presence of DNA damage patterns typical of ancient DNA.

We have a set of 10 modern and 26 ancient individuals (subsampled to 100k reads)
<pre>
find /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ -name "*bam"
</pre>

First, run mapDamage on one of the modern individuals:

<pre>
mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/modern/NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.bam -r /home/databases/references/human/hs37d5.fa --no-stats
</pre>

Examine the output:

<pre>
cd NA20786.mapped.ILLUMINA.bwa.TSI.low_coverage.20130415.100k_ss.mapDamage/
okular Length_plot.pdf &
okular Fragmisincorporation_plot.pdf &
cd ..
</pre>

'''Q1:''' which fragment length occurs most frequently?

'''Q2:''' what is the frequency of 5' C>T and 3' G>A substitutions ()

Run mapDamage on one of the ancient individuals
<pre>
mapDamage -i /home/projects/22126_NGS/exercises/adna/01_authentication/bam/ancient/allentoft_2015/RISE559.sort.rmdup.realign.md.100k.bam -r /home/databases/references/human/hs37d5.fa --no-stats
</pre>

Examine the output
<pre>
cd RISE559.sort.rmdup.realign.md.100k.mapDamage/
okular Length_plot.pdf &
okular Fragmisincorporation_plot.pdf &
</pre>

'''Q3:''' At what fragment length does the distribution show its peak?

'''Q4:''' what are the frequencies of 5' C>T (red line) and 3' G>A substitutions (blue line)?

'''Q5:''' which bases are enriched at 5' flanking position?

'''Q6:''' does your sample look ancient? if not, what might be the reason?

<H2> Population genetics </H2>

Create a new subdirectory and navigate to it:
<pre>
cd ..
mkdir 02_popgen
cd 02_popgen
</pre>

<H3>Explore the reference panel dataset</H3>

Pur reference panel dataset is in binary PLINK format, a widely used format in genetic studies (see documentation [https://www.cog-genomics.org/plink/1.9/ here]). We need to access the following files:

<pre>
ls /home/projects/22126_NGS/exercises/adna/02_popgen/plink/
</pre>

However, instead of copying them, we will create symbolic links using the ln command, these acts as placeholders and tell the operating system to pretend that there is an actual file there. This saves considerable disk space compared to copying over the files.

<pre>
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/plink/world.bed .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/plink/world.bim .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/plink/world.cluster .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/plink/world.fam .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/plink/world.sampleInfo.txt .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/eur.poplist .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/modern.poplist .
ln -s /home/projects/22126_NGS/exercises/adna/02_popgen/noneur.poplist .
</pre>

The PLINK binary format consists of 3 files:

{| class="wikitable"
| '''file'''
| '''description'''
|-
| world.bed
| | genotype data in binary format ('''not to be confused with genomic intervals bed file but it is confusing''')
|-
| world.bim
| metadata for the variants, 1 line per variant
|-
| world.fam
| metadata for the samples, 1 line per sample
|-

We also have the following files than contain extra information:

{| class="wikitable"
| '''file'''
| '''description'''
|-
|world.cluster
| pre-defined population groupings for samples (for plink)
|-
| world.sampleInfo.txt
| additional sample metadata (for plotting etc)
|}

Let us explore the metadata files:

<pre>
head world.fam
head world.bim
head world.cluster
head world.sampleInfo.txt
</pre>

'''Q7:''' How many samples / SNPs are in our dataset?

'''Q8:''' what populations are in our reference panel and what sample size do they have (trick: forgo the header using "tail -n+2", you need "sort" and uniq (prints 1 instance per repeated line), to tell "uniq" to count and print how many lines were repeated "-c"?

Calculate basic summary statistics (a simple description of the data) for the dataset:

<pre>
/home/ctools/plink-1.9/plink --bfile world --missing --out world
</pre>

'''Q9:''' are you getting the same number of variants and individuals as you did via UNIX command lines?

The world.imiss file lists the number and fraction of missing genotypes for each sample

'''Q10:''' what fraction of SNPs have a missing genotype for the Tyrolean Iceman?

<H3>Genotype and merge an ancient individual</H3>

In this section, we will merge our ancient data with the reference panel to prepare our dataset for downstream analysis genotypes for our ancient data will be obtained by randomly sampling a read from the alignments (BAM files) at the reference dataset SNP positions.

We are going to use a low-coverage individual from [https://pubmed.ncbi.nlm.nih.gov/26062507/ Allentoft et al (RISE507)], this data was obtained from an ~5100-year-old individual from the Early Bronze Age [https://en.wikipedia.org/wiki/Afanasievo_culture Afanasievo culture] in the Altai Mountains region

<pre>
ls /home/projects/22126_NGS/exercises/adna/02_popgen/bam/
</pre>

First, we need to extract a genomic interval bed file for the SNP positions of the reference panel:
<pre>
awk '{print $1"\t"($4-1)"\t"$4}' world.bim | gzip > world.snps.bed.gz
</pre>

awk is a command to create small programs. In this example, we tell it, print the first columns, the fourth column minus 1 and the fourth column again.

Inspect the results:

<pre>
zcat world.snps.bed.gz | head
</pre>

Create a read pileup file for the reference panel SNP positions (might take a few minutes)

<pre>
samtools mpileup -f /home/databases/references/human/hs37d5.fa -B -l world.snps.bed.gz /home/projects/22126_NGS/exercises/adna/02_popgen/bam/RISE507.sort.rmdup.realign.md.bam |gzip > RISE507.pileup.gz
</pre>

Examine the output:

<pre>
zcat RISE507.pileup.gz |head
</pre>

'''Q11''': how many SNPs of the reference panel are covered in RISE507?

Now we will randomly sample a DNA fragment at each position and output the results in VCF format (custom python script):
<pre>
zcat RISE507.pileup.gz | python2 /home/projects/22126_NGS/exercises/adna/02_popgen/get_haploid_vcf_from_pileup.py -r -s RISE507 |bgzip -c > RISE507.vcf.gz
</pre>
This is done because the coverage is insufficient to ensure proper genotyping.

Let us inspect the result:
<pre>
zcat RISE507.vcf.gz |grep -v "^#" |head
</pre>

We convert to plink binary format:
<pre>
/home/ctools/plink-1.9/plink --vcf RISE507.vcf.gz --make-bed --double-id --out RISE507
</pre>

Try to merge the sample with the reference panel
<pre>
/home/ctools/plink-1.9/plink --bfile world --bmerge RISE507 --out RISE507.merge
</pre>

You should get an error.

'''Q12''': how many SNPs failed the merge? What is the likely reason?

We will remove the failing SNPs and try again
<pre>
/home/ctools/plink-1.9/plink --bfile RISE507 --exclude RISE507.merge.missnp --make-bed --out RISE507.merge2
/home/ctools/plink-1.9/plink --bfile world --bmerge RISE507.merge2 --out RISE507.world
</pre>

Make a cluster file for subsetting
<pre>
awk '{print $1,$2,$1}' RISE507.world.fam > RISE507.world.cluster
</pre>

<H3>Investigate the genetic affinities of the ancient sample using PCA</H3>

In this section, we will try to place our sample within a PCA of a set of modern and ancient individuals.

First, we will have a look at the modern populations in the reference panel:
<pre>
/home/ctools/plink-1.9/plink --bfile RISE507.world --keep-clusters modern.poplist --within RISE507.world.cluster --pca header tabs --out modern
</pre>

We can plot the first two principal components using the custom R script plotPca.R

The three positional arguments are the eigenvector file, sample info file and prefix for the output

<pre>
Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R modern.eigenvec world.sampleInfo.txt modern
evince modern.pca.plot.pdf &
</pre>

'''Q13:''' which populations are most differentiated along PC1?
'''Q14:''' which populations are most differentiated along PC2?

We repeat the exercise on a subset of European populations:

<pre>
/home/ctools/plink-1.9/plink --bfile RISE507.world --keep-clusters eur.poplist --within RISE507.world.cluster --pca header tabs --out eur
Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R eur.eigenvec world.sampleInfo.txt eur
evince eur.pca.plot.pdf &
</pre>

'''Q15:''' which populations are most differentiated along PC1?
'''Q16:''' which populations are most differentiated along PC2?

Now, let us examine how the cluster of ancient individuals compared to the modern ones:

<pre>
/home/ctools/plink-1.9/plink --bfile RISE507.world --pca header tabs --out ancient.world
Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.world.eigenvec world.sampleInfo.txt ancient.world
evince ancient.world.pca.plot.pdf &
</pre>

Here are some references if you want to read more about the different ancient samples:

{| class="wikitable"
| '''sample'''
| '''link'''
|-
| UstIshim
| [https://en.wikipedia.org/wiki/Ust%27-Ishim_man]
|-
| Loschbour
| [https://en.wikipedia.org/wiki/Loschbour_man] [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4170574/]
|-
| Brana
| [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4269527/]
|-
| NE1
| [https://www.pnas.org/content/113/2/368]
|-
|Stuttgart
| [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4170574/]
|-
| Iceman
| [https://www.iceman.it/en/the-iceman/]
|-
|Karelia
| [https://en.wikipedia.org/wiki/Karelians]
|-
| Samara
| [https://en.wikipedia.org/wiki/Samara_culture]
|-
| MA1
| [https://en.wikipedia.org/wiki/Mal%27ta%E2%80%93Buret%27_culture]
|-
| RISE507
|[https://pubmed.ncbi.nlm.nih.gov/26062507/]
|}

'''Q17:''' which ancient individuals don't cluster close to any modern individuals? what could be a plausible reason?

Repeat the exercise but remove the non-European modern individuals:

<pre>
/home/ctools/plink-1.9/plink --bfile RISE507.world --within RISE507.world.cluster --remove-clusters noneur.poplist --pca header tabs --out ancient.eur
Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient.eur.eigenvec world.sampleInfo.txt ancient.eur
evince ancient.eur.pca.plot.pdf &
</pre>

'''Q18:''' which populations are most differentiated along PC1? what could be a plausible reason?

As a final exercise, we now project the ancient individual on PCs inferred from modern Europeans:

<pre>
/home/ctools/plink-1.9/plink --bfile RISE507.world --within RISE507.world.cluster --pca-clusters eur.poplist --remove-clusters noneur.poplist --pca header tabs --out ancient_proj.eur --maf 0.01
Rscript /home/projects/22126_NGS/exercises/adna/02_popgen/plotPca.R ancient_proj.eur.eigenvec world.sampleInfo.txt ancient_proj.eur
evince ancient_proj.eur.pca.plot.pdf &
</pre>

'''Q19:''' where does our study individual cluster now?

'''Q20:''' How do you explain that an individual that is found closer to the modern-day Chinese border is closer to modern Europeans than he is to the Han Chinese?

Please find answers [[Ancient_DNA_exercise_answers|here]]