In eckartbindewald/bfxapps2: R Bioinformatics Applications

knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE)

Preliminaries

if (!require(rrBLUP)) {
  install.packages("rrBLUP")
}
if (!require(VariantAnnotation)) {
  BiocManager::install("VariantAnnotation")
}

Loading Libraries

library(rrBLUP)
library(VariantAnnotation)

remote_folder="https://raw.githubusercontent.com/PacktPublishing/R-Bioinformatics-Cookbook/master/datasets/ch2" # was Chapter02
vcf_file = "small_sample.vcf"
if (!file.exists(vcf_file)) {
  download.file(file.path(remote_folder,vcf_file),vcf_file)
}
stopifnot(file.exists(vcf_file))
vcf <- readVcf(vcf_file,"hg19")
vcf

Extract Genotype, Sample and Marker Position

gts <- geno(vcf)$GT
gts

Extracting Needed Information

samples <- samples(header(vcf))
markers <- rownames(gts)
chrom <- as.character(seqnames(rowRanges(vcf)))
pos <- as.numeric(start(rowRanges(vcf)))
samples
markers
chrom
pos

Different Encodings - A Workaround

• In VCF file format, we encode mutations as "0/0", "0/1", "1/0", "1/1"
• Corresponds to AA (homozygous), Aa and aA (heterozygous) and "aa"
• GWAS function from rrBLUP package needs different encoding (numbers 1,0,0,-1 respectively)
• We need as workaround convert those two encodings

Create Custom Function

convert <- function(v){
 v <- gsub("0/0",1,v)
 v <- gsub("0/1", 0, v)
 v <- gsub("1/0", 0, v)
 v <- gsub("1/1", -1, v)
 return(v)
}

Convert Data to Matrix

gt_char <- apply(gts, convert, MARGIN=2)
genotype_matrix <- matrix(as.numeric(gt_char), nrow(gt_char))
colnames(genotype_matrix) <- samples
genotype_matrix

Build Dataframe

variant_info <- data.frame(marker = markers,
                           chrom = chrom,
                           pos = pos)
variant_info

Build Combined Variant/Genotype Dataframe

genotypes <- cbind(variant_info, as.data.frame(genotype_matrix))
genotypes

Build Phenotype Dataframe

phenotypes <- data.frame(line=samples, score=rnorm(length(samples)))
phenotypes

Run GWAS

gresult=GWAS(phenotypes, genotypes, plot=FALSE)
gresult

Visualizing Manhattan Plot

library(dplyr)
library(qqman)
gresult2 <- gresult %>% rename(SNP=marker, CHR=chrom, BP=pos) %>% mutate(CHR=as.numeric(CHR), P=10^(-score)) %>% mutate(score=NULL)
manhattan(gresult2, xlim=c(10000,20000))

Practical Exercise

if (!require(mlmm.gwas)) {
  install.packages("mlmm.gwas")
}
library(mlmm.gwas)
data("mlmm.gwas.AD")

eckartbindewald/bfxapps2 documentation built on Feb. 6, 2025, 3:22 a.m.