Homework 6: Ancestry

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! But please run all commands to produce 
### in the end a PCA plot containing yourself compared to all 1000 Genomes 
### samples.

# 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
module load tabix_0.2.6


###############
### 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
### because we're going to conduct PCA on these SNPs and only want common ones.
### 
###------ QUESTION 3: WHY WOULD WE RETAIN ONLY 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/PCA/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