SNP calling exercise part 1 - Revision history

Mick at 10:00, 7 January 2026

2026-01-07T10:00:14Z

← Older revision		Revision as of 12:00, 7 January 2026
Line 131:		Line 131:

	<pre>		<pre>
	bcftools stats [INPUT_VCF]		/home/ctools/bcftools-1.23/bcftools stats [INPUT_VCF]
	</pre>		</pre>

Line 173:		Line 173:

	<pre>		<pre>
	bcftools view		/home/ctools/bcftools-1.23/bcftools view
	</pre>		</pre>

Mick at 09:11, 7 January 2026

2026-01-07T09:11:02Z

← Older revision		Revision as of 11:11, 7 January 2026
Line 90:		Line 90:

	<pre>		<pre>
	tabix -f -p vcf [INPUT_GVCF]		/home/ctools/htslib-1.20/tabix -f -p vcf [INPUT_GVCF]
	</pre>		</pre>

Mick at 10:37, 5 January 2026

2026-01-05T10:37:16Z

← Older revision		Revision as of 12:37, 5 January 2026
Line 34:		Line 34:

	<pre>		<pre>
	gatk --java-options "-Xmx10g" HaplotypeCaller \		/home/ctools/gatk-4.6.2.0/gatk --java-options "-Xmx10g" HaplotypeCaller \
	-R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa \		-R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa \
	-I [INPUT_BAM] \		-I [INPUT_BAM] \
Line 100:		Line 100:

	<pre>		<pre>
	gatk GenotypeGVCFs \		/home/ctools/gatk-4.6.2.0/gatk GenotypeGVCFs \
	-R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa \		-R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa \
	-V [INPUT_GVCF] \		-V [INPUT_GVCF] \
Line 119:		Line 119:

	<pre>		<pre>
	tabix -f -p vcf NA24694.vcf.gz		/home/ctools/htslib-1.20/tabix -f -p vcf NA24694.vcf.gz
	</pre>		</pre>

Line 147:		Line 147:

	<pre>		<pre>
	tabix [INPUT_VCF] [CHROM]:[START]-[END]		/home/ctools/htslib-1.20/tabix [INPUT_VCF] [CHROM]:[START]-[END]
	</pre>		</pre>

Line 153:		Line 153:

	<pre>		<pre>
	tabix [INPUT_VCF] [CHROM]:[POS]-[POS]		/home/ctools/htslib-1.20/tabix [INPUT_VCF] [CHROM]:[POS]-[POS]
	</pre>		</pre>

Line 201:		Line 201:

	<pre>		<pre>
	tabix NA24694.vcf.gz chr20:32011209-32011209		/home/ctools/htslib-1.20/tabix NA24694.vcf.gz chr20:32011209-32011209
	tabix NA24694.vcf.gz chr20:32044279-32044279		/home/ctools/htslib-1.20/tabix NA24694.vcf.gz chr20:32044279-32044279
	</pre>		</pre>

Mick at 13:49, 25 November 2025

2025-11-25T13:49:48Z

← Older revision		Revision as of 15:49, 25 November 2025
Line 261:		Line 261:

	<p><b>Congratulations, you finished the exercise!</b></p>		<p><b>Congratulations, you finished the exercise!</b></p>

	~~<p>You can find suggested answers here: <a href="SNP_calling_exercise_part_1_answers_">SNP_calling_exercise_part_1_answers_</a></p>~~

Mick at 13:48, 25 November 2025

2025-11-25T13:48:36Z

← Older revision		Revision as of 15:48, 25 November 2025
Line 23:		Line 23:
	<p>We will genotype chromosome 20 from a BAM file that has been pre-processed (sorted, duplicate-marked, etc.)</p>		<p>We will genotype chromosome 20 from a BAM file that has been pre-processed (sorted, duplicate-marked, etc.)</p>

	<p>The sample is a ~~<a href="~~https://en.wikipedia.org/wiki/Han_Chinese">Han Chinese~~</a>~~ male, sequenced to approximately 24.6× coverage.</p>		<p>The sample is a [https://en.wikipedia.org/wiki/Han_Chinese Han Chinese] male, sequenced to approximately 24.6x coverage.</p>

	<pre>		<pre>

Mick at 13:46, 25 November 2025

2025-11-25T13:46:43Z

Show changes

Gabre at 15:27, 15 December 2024

2024-12-15T15:27:20Z

← Older revision		Revision as of 17:27, 15 December 2024
Line 139:		Line 139:
	Contrast the number you obtained for questions 6 to the one you obtained for question 3. Do you think it's expected?		Contrast the number you obtained for questions 6 to the one you obtained for question 3. Do you think it's expected?

			'''Congratulations you finished the exercise!'''

	Please find answers [[~~SNP_calling_exercise_answers_part_1~~\|here]]		Please find answers [[SNP_calling_exercise_part_1_answers_\|here]]

Gabre at 13:36, 29 November 2024

2024-11-29T13:36:27Z

← Older revision		Revision as of 15:36, 29 November 2024
Line 138:		Line 138:

	Contrast the number you obtained for questions 6 to the one you obtained for question 3. Do you think it's expected?		Contrast the number you obtained for questions 6 to the one you obtained for question 3. Do you think it's expected?


			Please find answers [[SNP_calling_exercise_answers_part_1\|here]]

Gabre: Created page with "

Overview

First:

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

We will:

Genotype some whole-genome sequencing data.
Get acquainted with VCF files
Soft filtering
Hard filtering
Annotation of variants

----

Genotyping

We will genotype a chromosome from a BAM file that has been processed using the steps we detailed before. It i..."

2024-11-29T13:33:29Z

Created page with "<H2>Overview</H2> First: <OL> <LI>Navigate to your home directory: <LI>Create a directory called "variant_call" <LI>Navigate to the directory you just created. </OL> We will: <OL> <LI>Genotype some whole-genome sequencing data. <LI>Get acquainted with VCF files <LI>Soft filtering <LI>Hard filtering <LI> Annotation of variants </OL> ---- <H2>Genotyping</H2> We will genotype a chromosome from a BAM file that has been processed using the steps we detailed before. It i..."

New page

<H2>Overview</H2>

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

We will:
<OL>
<LI>Genotype some whole-genome sequencing data.
<LI>Get acquainted with VCF files
<LI>Soft filtering
<LI>Hard filtering
<LI> Annotation of variants
</OL>

----

<H2>Genotyping</H2>

We will genotype a chromosome from a BAM file that has been processed using the steps we detailed before. It is from a [https://en.wikipedia.org/wiki/Han_Chinese Han Chinese] Male with a depth of coverage of 24.6X.
<pre>
/home/projects/22126_NGS/exercises/snp_calling/NA24694.bam
</pre>

It has been indexed. We will first use GATK's HaplotypeCaller, the command-line look something like this:

<pre>
gatk --java-options "-Xmx10g" HaplotypeCaller -R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa -I [INPUT BAM] -L chr20 -O [OUTPUT ] --dbsnp /home/databases/databases/GRCh38/Homo_sapiens_assembly38.dbsnp138.vcf.gz -ERC GVCF
</pre>

" -R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa " specifies a reference, "-L chr20" specifies only to consider chromosome 20, "-dbsnp /home/databases/databases/GRCh38/Homo_sapiens_assembly38.dbsnp138.vcf.gz ", add annotation from SNPs which are known to exist in the human population (see [https://www.ncbi.nlm.nih.gov/snp/ dbSNP]). We should point out that the variation coming from Eurasians is more extensively represented in this database. However, the most genetically diverse populations are found in [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1740-1 Africa]. Finally, " -ERC GVCF" will generate a GVCF (Genomic Variant Call Format) which essentially contains all sites instead of just variant sites. The input is specified above, you can call the output: NA24694.gvcf.gz

The command takes between '''20 to 30 minutes to run''', so feel free to take a break. Or feel free to copy the one I generated:

<pre>
/home/projects/22126_NGS/exercises/snp_calling/NA24694.gvcf.gz
</pre>

We will have a brief look at the file:

<pre>
zcat NA24694.gvcf.gz|less -S
</pre>

You will see several lines starting with '#', this is the header portion, scroll down to see the calls. The first 5 columns represent, the name of the chromosome, the coordinate, snp ID (from dbSNP), the reference base and the alternative base. You should see a mix of sites, some with the mention <NON_REF> in the 5th column which corresponds to the alternative base, those sites are probably invariant, and others where there is a variant in the fifth column.

Prior to running GATK, we need to index it:
<pre>
tabix -f -p vcf [VCF to index]
</pre>

Next, we will use the preliminary variant to actually generate only variant sites:
<pre>
gatk GenotypeGVCFs -R /home/databases/references/human/GRCh38_full_analysis_set_plus_decoy_hla.fa -V [INPUT GVCF] -O [OUTPUT VCF] -L chr20 --dbsnp /home/databases/databases/GRCh38/Homo_sapiens_assembly38.dbsnp138.vcf.gz
</pre>

In your case, the input is the file you generated above and the output can be called NA24694.vcf.gz.
This command takes less than 5 minutes to run.

Then let us index the file with the same tabix command:
<pre>
tabix -f -p vcf [VCF to index]
</pre>

Just like "samtools index" can allow us to create an index to retrieve portions of a BAM file, tabix is another utility that allows us to retrieve portions of a GVCF file.

<H3> Get acquainted with VCF files </H3>

'''Q1'''

Using:
<pre>
bcftools stats [input vcf]
</pre>
How many SNPs do you have? (hint: find the line with: "SN 0 number of SNPs")

<pre>
tabix -f -p vcf [VCF to index]
</pre>

The index will be stored as a "tbi" file. Then, the VCF file can be queried:

<pre>
tabix [VCF to index] [CHROMOSOME]:[START]-[END]
</pre>

where [CHROMOSOME] is the chromosome name and [START] [END] are the start and end coordinates respectively.
Or for a single coordinate:
<pre>
tabix [VCF to index] [CHROMOSOME]:[COORDINATE]-[COORDINATE]
</pre>

'''Q2'''

Use the command above, determine how many total variants are in the 1 million base pair region "chr20:32000000-33000000"? (hint, remember "wc -l" to count lines).

'''Q3'''

bcftools is a nifty utility that allows us to do various operations on VCF files. Type:

<pre>
bcftools view
</pre>

It contains several options. Try to determine how many SNPs (excluding indels and multi-allelic sites) are there in the '''same region''' as above? (hint: you need to filter for a certain '''type''' of variant, hint2: be careful not to include the header in the number of lines (option: -H)).

'''Q4'''

Use either tabix or bcftools to retrieve the SNP located on chromosome 20 at coordinate 32011209. Tell me in your own words what is the genotype? What about chromosome 20, coordinate 32044279. Become familiar with the different fields in the [https://samtools.github.io/hts-specs/VCFv4.2.pdf VCF specifications]. Pay attention to "1.4 Data lines". For the genotype fields (column 9 and 10) there are more info on the GATK fields: [https://gatk.broadinstitute.org/hc/en-us/articles/360035531692-VCF-Variant-Call-Format GATK's VCF-Variant-Call-Format].

Also, answer for each site:

<OL>
<LI> what is the allele depth i.e. how many bases are there of each type?
<LI> what is the depth of coverage i.e. how many bases cover this site?
<LI> What is your genotype quality?
<LI> What are your genotype likelihoods?
</OL>

'''Q5'''

Inspect the SNPS at positions:
<pre>
chr20 32974911
chr20 64291638
</pre>

One SNP has poor quality the other has good quality. Which is which? Why do you think this is? (hint: remember the class, the more data you have, the more you are sure).

'''Q6'''

Using the same region, 32000000-33000000, how many are novel (i.e. does not exist in databases of a large number of sampled genomed) SNPs? Note that the third column is the ID of a known SNP in dbSNP. You can use the '''cut''' command to get a specific column. Also, every ID in dbSNP starts with "rs" and you could use '''grep''' to either retain such lines or filter them out (see option -v).

'''Q7'''

Contrast the number you obtained for questions 6 to the one you obtained for question 3. Do you think it's expected?