Repository navigation
Unix VI: Disease causing mutations
Course: HT17 Unix with applications to next generation sequencing data (SC00025)
In this exercise, we will prioritize possible disease causing mutations
Copy Samples.vcf.tar.gz from /home/marcela/unix/files and untar it. You will get five VCF files. Here is a quick reference to the VCF format.
Inspect the S1.vcf file and answer:
Q1. How many SNPs do you have per chromosome?
Q2. Are there any insertions or deletions? If so, how many?
The 7th column (FILTER) in a VCF file indicates which SNPs have passed certain thresholds
(if you grep FILTER
on your VCF file, you can see which filters have been used)
Q3. How many SNPs have passed all these filters?
The 10th column in a VCF file, has information about the predicted genotype of the mutation:
- 0/1 indicates Heterozygous, while
- 1/1 indicates Homozygous
Save the homozygous mutations that have passed all thresholds into a new file (S1.pass.homozygous.vcf),
including only the following columns: CHROM, POS, REF and ALT, so it looks like:
chr1 131263 CCAGT C
chr1 162463 T C
chr1 714427 G A
chr1 715348 T G
chr1 715996 A T
chr1 723891 G C
Sometimes we have a lot of samples to process in the same way. In this cases, we can use a for-loop to speed up things:
for i in *vcf
do statement
done
Create filtered files for the other samples (S2, S3, S4 and S5) using a for-loop.
Q4. What is your command line?
It is common to share this lists among your peers. So let's practice to pack and compress the files you have created:
tar -cvzf SNP_pass_homozygous.tar.gz *pass.homozygous.vcf
Q5. How many confident homozygous SNPs in chr14 are found in all samples?
Keep working with S1.vcf. A common next step is to annotate our candidate SNPs to give priority to those
SNPs in coding regions, especially those which nucleotide changes make an impact the aminoacid
composition of the protein. As usual, we need to modify our file to get the correct format to use
for the annotation. We will use the ANNOVAR package to annotate the SNPs.
First use /home/marcela/bin/convert2annovar.pl and convert your vcf file to annovar format.
Use the following arguments:
* format: vcf4
* filter: pass
* coverage: 50
Now, annotate the output file with /home/marcela/bin/table_annovar.pl using the following settings:
* database-location: /home/marcela/bin/humandb/
* buildver: hg19
* protocol: refGene
* operation: g
Don't forget to understand what these options mean!
Several files will be created, for this exercise focus on the one with
extension annovar.hg19_multianno.txt. Have a look at the file, if you need more
info to understand it, go to the ANNOVAR webpage.
Q6. What are the top most mutated genes in S1 that harbor aminoacid changes?
Created by Marcela Dávila