==== 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! 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