Difference between revisions of "Homework 6: Ancestry"
Jump to navigation
Jump to search
(→Notes) |
|||
| Line 76: | Line 76: | ||
### Add in rsIDs from dbSNP. PLINK needs these to reconcile | ### Add in rsIDs from dbSNP. PLINK needs these to reconcile | ||
### positions/alleles | ### positions/alleles | ||
| − | vcf-add-rsid -vcf chrALL.filtered.PASS.beagled.HG00096.vcf.gz -db /Users/scvr9332/reference_data/dbsnp.144.b37.vcf.gz | bgzip -c > chrALL.filtered.PASS.beagled.HG00096.rsID.vcf.gz | + | vcf-add-rsid -vcf chrALL.filtered.PASS.beagled.HG00096.vcf.gz -db /Users/scvr9332/reference_data/dbsnp_144/dbsnp.144.b37.vcf.gz | bgzip -c > chrALL.filtered.PASS.beagled.HG00096.rsID.vcf.gz |
### The previous command keeps only variants with rsIDs, otherwise | ### The previous command keeps only variants with rsIDs, otherwise | ||
### plink throws an error that there are >1 variants with ID = "." | ### plink throws an error that there are >1 variants with ID = "." | ||
Revision as of 18:31, 15 November 2015
Notes
Code for Homework
# PSYC 7102 -- Statistical Genetics
# Homework #6: Ancestry
# Due: November 25th
# PLEASE CONTACT DAVID TO GET YOUR NEW VCF.
# Overview: Using your variant calls contained in your VCF files,
# combined with the existing 1000 Genomes variant calls from the 1000
# Genomes project, you will conduct a principal components analysis
# of the 1000 Genomes samples using PLINK and estimate your position in
# the 1000G PCA space
# This homework is a little bit different because the commands are
# slightly complicated, so I've written them all for you. In addition
# to specific conceptual questions, I ask you to run the commands
# below (changing directories appropriately as needed) and then
# explain to me what the command did and why we did it. I don't need
# super technical details but simply a solid conceptual overview of
# the purpose of the command.
### ONLY EXPLAIN COMMANDS WHERE I SPECIFICALLY REQUEST IT! YOU DO NOT
### HAVE TO EXPLAIN EVERY COMMAND!
# For many questions you'll want to run analyses by chromosome. To do
# this, please create an interactive PBS session like this:
qsub -I -l walltime=23:00:00 -q short -l mem=10gb -l nodes=1:ppn=22
### Load apigenome, plink
module load apigenome_0.0.2
module load plink_latest
###############
### STEP 1 ###
###############
# Here I take a minute to describe the various files you'll use to get started
# The 1000 Genomes vcfs from the 1000 Genomes Consortium
/Users/scvr9332/1000Gp3/variant_calls/ALL.chr[1-22].phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
# A vcf file with, depending on your level of comfort, your own
# genotypes or 1kg sample HG000096. All examples that follow use
# HG00096 for convenience, which can be obtained as follows. To use
# your own sample just replace the $10 with $NF, assuming your
# sample is in the last column of the vcf.
(zcat chr1.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr2.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr3.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr4.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr5.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr6.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr7.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr8.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr9.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr10.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr11.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr12.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr13.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr14.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr15.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr16.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr17.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr18.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr19.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr20.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr21.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}'; \
zgrep -v '#' chr22.filtered.PASS.beagled.vcf.gz | gawk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}') | bgzip -c > chrALL.filtered.PASS.beagled.HG00096.vcf.gz
### Add in rsIDs from dbSNP. PLINK needs these to reconcile
### positions/alleles
vcf-add-rsid -vcf chrALL.filtered.PASS.beagled.HG00096.vcf.gz -db /Users/scvr9332/reference_data/dbsnp_144/dbsnp.144.b37.vcf.gz | bgzip -c > chrALL.filtered.PASS.beagled.HG00096.rsID.vcf.gz
### The previous command keeps only variants with rsIDs, otherwise
### plink throws an error that there are >1 variants with ID = "."
###
###------ QUESTION 1: WHAT DOES THIS COMMAND DO? (1 point)
zgrep -E '#|\srs[0-9]' chrALL.filtered.PASS.beagled.HG00096.rsID.vcf.gz | bgzip -c > chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz
##############
### STEP 2 ###
##############
### Whole bunch of data cleaning. Retain in the 1000 Genomes
### Consortium VCFs only biallelic SNPs that also exist in your vcf
### files. PLINK 1.9 can't handle all variant types gracefully.
# Get list of your SNPs
zcat chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz | cut -f3 | grep '^rs' > MyrsIDs.txt
#
# Reduce the 1000 Genomes VCFs to bi-allelic SNPs (that's all PLINK
# can handle anyway). BE SURE TO LOG INTO A COMPUTE NODE BEFORE
# RUNNING ANY OF THESE LOOPS!!!!!
###
###------ QUESTION 2: WHAT DOES THIS COMMAND (THE ENTIRE FOR LOOP) DO? (2 points)
for i in {1..22}; do
zgrep -E '^#|\s[ACTG]\s[ACTG]\s' /Users/scvr9332/1000Gp3/variant_calls/ALL.chr$i.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | bgzip -c > chr$i.1000G.biallelic.vcf.gz &
done
### Retain in the 1000 Genomes VCF only your SNPs that are also fairly common
###
###------ QUESTION 3: WHY WOULD WE REMOVE COMMON SNPS, OTHER THAN IT
###------ MAKES EVERY COMMAND LATER FASTER? (2 points)
for i in {1..22}; do
plink --vcf chr$i.1000G.biallelic.vcf.gz \
--extract MyrsIDs.txt \
--maf .05 \
--make-bed \
--out chr$i.1000G.biallelic.common.plink &
done
### Remove any remaining duplicates from the 1000G reference files
### A) Identify duplicate sites
###
###------ QUESTION 4: WALK ME THROUGH WHAT THIS COMMAND (THE ENTIRE
###------ FOR LOOP) IS DOING. (3 points)
for i in {1..22}; do
zcat chr$i.1000G.biallelic.vcf.gz | cut -f3 | sort | uniq -c | grep -E '\s2\s' | grep '^rs' | gawk '{print $2}' > dups$i.txt &
done
### B) Remove the duplicate sites
for i in {1..22}; do
plink --bfile chr$i.1000G.biallelic.common.plink \
--exclude dups$i.txt \
--make-bed \
--out chr$i.1000G.biallelic.common.plink.nodups &
done
rm chr*.1000G.biallelic.vcf.gz
rm chr*.1000G.biallelic.common.plink.bed
rm chr*.1000G.biallelic.common.plink.bim
rm chr*.1000G.biallelic.common.plink.fam
## Identify variants in LD
for i in {1..22}; do
plink --bfile chr$i.1000G.biallelic.common.plink.nodups \
--indep-pairwise 100k 5 .01 \ #------ QUESTION 5: EXPLAIN WHAT THIS ARGUMENT IS DOING (1 point)
--maf .05 \
--out chr$i.snpsToExtract &
done
## Extract LD-independent variants
for i in {1..22}; do
plink --bfile chr$i.1000G.biallelic.common.plink.nodups \
--extract chr$i.snpsToExtract.prune.in \
--make-bed \
--out chr$i.1000G.biallelic.common.plink.nodups.LDpruned &
done
rm chr*.1000G.biallelic.common.plink.nodups.bed
rm chr*.1000G.biallelic.common.plink.nodups.bim
rm chr*.1000G.biallelic.common.plink.nodups.fam
rm dups*
##############
### STEP 3 ###
##############
### Merge the new and improved 1000 Genomes plink files into a single
### file with all chromosomes. This is easier to do when they're vcf
### files
### Example Command:
for i in {1..22}; do
plink --bfile chr$i.1000G.biallelic.common.plink.nodups.LDpruned --recode-vcf --out tmp$i &
done
(cat tmp1.vcf; grep -v '#' tmp2.vcf; grep -v '#' tmp3.vcf; grep -v '#' tmp4.vcf; grep -v '#' tmp5.vcf; \
grep -v '#' tmp6.vcf; grep -v '#' tmp7.vcf; grep -v '#' tmp8.vcf; grep -v '#' tmp9.vcf; grep -v '#' tmp10.vcf; \
grep -v '#' tmp11.vcf; grep -v '#' tmp12.vcf; grep -v '#' tmp13.vcf; grep -v '#' tmp14.vcf; grep -v '#' tmp15.vcf; \
grep -v '#' tmp16.vcf; grep -v '#' tmp17.vcf; grep -v '#' tmp18.vcf; grep -v '#' tmp19.vcf; grep -v '#' tmp20.vcf; \
grep -v '#' tmp21.vcf; grep -v '#' tmp22.vcf) | bgzip -c > chrALL.1000G.biallelic.common.plink.nodups.LDpruned.vcf.gz
###----- QUESTION 6: HOW MANY VARIANTS ARE IN THE OUTPUT FILE FROM THE
###----- COMMAND ABOVE? (1 point)
plink --vcf chrALL.1000G.biallelic.common.plink.nodups.LDpruned.vcf.gz \
--make-bed \
--out chrALL.1000G.biallelic.common.plink.nodups.LDpruned
rm tmp*
rm chr[0-9]*.bed
rm chr[0-9]*.bim
rm chr[0-9]*.fam
rm chr[0-9]*.nosex
rm chr[0-9]*.log
rm chr[0-9]*.prune.in
rm chr[0-9]*.prune.out
##############
### STEP 4 ###
##############
# Conduct a PCA on these new and improved 1000 Genomes plink files.
plink --bfile chrALL.1000G.biallelic.common.plink.nodups.LDpruned --pca var-wts
# Use R to organize the results from the PCA into a file that can then be
# imported into PLINK using the --score command
echo "### Extract per - SNP weights for each of the first 20 PCs and get the
### correct effect allele for later scoring
###
PCA <- read.table('plink.eigenvec.var', header = F, col.names = c('CHROM','RS',paste('PC',1:20,sep='')))
bim <- read.table('chrALL.1000G.biallelic.common.plink.nodups.LDpruned.bim', header = F, col.names = c('CHROM','RS','unknown','POS','A1','A2'))
x <- merge (PCA, bim)
write.table(subset(x, select = c('RS','A1', paste('PC', 1:20, sep = '' ))) , file = 'PCA-score-file.txt', quote=F)
### End: extract per - SNP weights for each of the first 20
" | R --vanilla
##############
### STEP 5 ###
##############
### Score yourself (or 1000 Genomes sample HG00096) on the PCs
### Calculate the score for each person by taking the PC-weighted
### average for each site in the score file.
###
###------ QUESTION 7: EXPLAIN WHAT THE FOLLOWING COMMANDS ARE DOING (2 points)
plink --vcf chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz --score PCA-score-file.txt 2 3 4 header --out PC1
plink --vcf chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz --score PCA-score-file.txt 2 3 5 header --out PC2
plink --vcf chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz --score PCA-score-file.txt 2 3 6 header --out PC3
plink --vcf chrALL.filtered.PASS.beagled.HG00096.rsIDsOnly.vcf.gz --score PCA-score-file.txt 2 3 7 header --out PC4
plink --bfile chrALL.1000G.biallelic.common.plink.nodups.LDpruned --score PCA-score-file.txt 2 3 4 header --out PC1.1kg
plink --bfile chrALL.1000G.biallelic.common.plink.nodups.LDpruned --score PCA-score-file.txt 2 3 5 header --out PC2.1kg
plink --bfile chrALL.1000G.biallelic.common.plink.nodups.LDpruned --score PCA-score-file.txt 2 3 6 header --out PC3.1kg
plink --bfile chrALL.1000G.biallelic.common.plink.nodups.LDpruned --score PCA-score-file.txt 2 3 7 header --out PC4.1kg
##############
### STEP 6 ###
##############
### Plot the results!
###
### To run the following command you'll have to install the R package
### 'car'. To do this on vieques please run the following. ONly need
### to run it once
mkdir /Users/YourUniqueName/R
export R_LIBS=/Users/scvr9332/R ## replace my unique name with yours
module load R_3.2.2
R
### Then in this R session run: 'install.packages('car', repos='http://cran.r-project.org')'
### If you log out and back in, you'll have to again run this command again to get the car library loaded:
export R_LIBS=/Users/scvr9332/R
###------ QUESTION 8: PLEASE INTERPRET THE RESULTING GRAPH. WHAT IS
###------ BEING PLOTTED? WHAT DOES IT SAY ABOUT THE
###------ ANCESTRY OF YOUR SAMPLE (OR WHATEVER SAMPLE YOU
###------ USED) (2 points)
echo "
### Read in scores for PC1 and PC2 in 1000g
kg_PC1 <- read.table('PC1.1kg.profile', header=T)
kg_PC2 <- read.table('PC2.1kg.profile', header=T)
kg_PC3 <- read.table('PC3.1kg.profile', header=T)
kg_PC4 <- read.table('PC4.1kg.profile', header=T)
### 'Self-reported' ancestry of 1000g participants
kg_sf <- read.table('/Users/scvr9332/20130502.sequence.index', header=T, sep='\t', fill=T, stringsAsFactors=F)
sample_ids <- unique(data.frame(IID=kg_sf$SAMPLE_NAME, POPULATION=kg_sf$POPULATION, stringsAsFactors=F))
sample_ids$CONTINENT <- ifelse(sample_ids$POPULATION=='CHB' |
sample_ids$POPULATION=='JPT' |
sample_ids$POPULATION=='CHS' |
sample_ids$POPULATION=='CDX' |
sample_ids$POPULATION=='KHV' |
sample_ids$POPULATION=='CHD', 'EAS', sample_ids$POPULATION)
sample_ids$CONTINENT <- ifelse(sample_ids$CONTINENT=='CEU' |
sample_ids$CONTINENT=='TSI' |
sample_ids$CONTINENT=='GBR' |
sample_ids$CONTINENT=='FIN' |
sample_ids$CONTINENT=='IBS', 'EUR', sample_ids$CONTINENT)
sample_ids$CONTINENT <- ifelse(sample_ids$CONTINENT=='YRI' |
sample_ids$CONTINENT=='LWK' |
sample_ids$CONTINENT=='GWD' |
sample_ids$CONTINENT=='MSL' |
sample_ids$CONTINENT=='ESN', 'AFR', sample_ids$CONTINENT)
sample_ids$CONTINENT <- ifelse(sample_ids$CONTINENT=='ASW' |
sample_ids$CONTINENT=='ACB' |
sample_ids$CONTINENT=='MXL' |
sample_ids$CONTINENT=='PUR' |
sample_ids$CONTINENT=='CLM' |
sample_ids$CONTINENT=='PEL', 'ADM', sample_ids$CONTINENT)
sample_ids$CONTINENT <- ifelse(sample_ids$CONTINENT=='GIH' |
sample_ids$CONTINENT=='PJL' |
sample_ids$CONTINENT=='BEB' |
sample_ids$CONTINENT=='STU' |
sample_ids$CONTINENT=='ITU', 'SAS', sample_ids$CONTINENT)
sample_ids$CONTINENT_numeric <- ifelse(sample_ids$CONTINENT == 'EAS', 1, sample_ids$CONTINENT)
sample_ids$CONTINENT_numeric <- ifelse(sample_ids$CONTINENT == 'EUR', 1, sample_ids$CONTINENT_numeric)
sample_ids$CONTINENT_numeric <- ifelse(sample_ids$CONTINENT == 'AFR', 1, sample_ids$CONTINENT_numeric)
sample_ids$CONTINENT_numeric <- ifelse(sample_ids$CONTINENT == 'ADM', 1, sample_ids$CONTINENT_numeric)
sample_ids$CONTINENT_numeric <- ifelse(sample_ids$CONTINENT == 'SAS', 1, sample_ids$CONTINENT_numeric)
kg_PC1 <- merge(kg_PC1, sample_ids)
My_PC1 <- read.table('PC1.profile', header=T)
My_PC2 <- read.table('PC2.profile', header=T)
My_PC3 <- read.table('PC3.profile', header=T)
My_PC4 <- read.table('PC4.profile', header=T)
kg <- data.frame(PC1 = -kg_PC1$SCORE,
PC2 = -kg_PC2$SCORE,
PC3 = -kg_PC3$SCORE,
PC4 = -kg_PC4$SCORE,
ancestry = kg_PC1$CONTINENT)
all <- rbind(kg, data.frame(PC1=-My_PC1$SCORE,
PC2=-My_PC2$SCORE,
PC3=-My_PC3$SCORE,
PC4=-My_PC4$SCORE,
ancestry=rep('Me', nrow(My_PC1))))
### Scatterplot matrix
library(car)
pdf('ancestry_Scatterplot.pdf')
scatterplotMatrix(~PC1+PC2+PC3+PC4 | ancestry, data=all, smoother=FALSE, reg.line=FALSE, cex=.7)
dev.off()
" | R --vanilla