BIOM262 - eQTL mapping

Genotype data and sample information

We will use genotyping data from the 1000 genomes project. Specifically we will use genotypes from Yoruba individuals from Nigeria

Read Samples

Read in the identifiers of the samples. Each sample is an immortalized lymphoblastoid cell line derived from a different Yoruba individual.

# use 'as.character' to make sure we have the sample names as a vector of character strings (rather than
# a factor)
sample.tab <- read.table("samples.txt", header=F)
sample.names <- as.character(sample.tab$V1)
head(sample.names)

[1] "NA18520" "NA19225" "NA19143" "NA19201" "NA19144" "NA19147"

# just cell line sample names, NA = nucleic acid

Read Genotypes

Read in SNP genotypes for each of the samples. For simplicity, we will only use SNPs on chromosome 22. The first 4 columns give information about the SNP, the remaining columns give the genotypes for each sample. The genotypes are coded as the number of copies of allele2 that each individual carries and can therefore be 0, 1, or 2. For convenience, also create a ‘genotype matrix’ that only contains columns with genotype counts.

genotypes.tab <- read.table("chr22_genotypes.txt", header=T)
head(genotypes.tab)

# all the individuals have geneotypes
# 0 = homozygous for reference allele (most common)
# 1 = heterozygous for reference allele (aka non-refernece allele or not the same in rerence genome)
# 2 = homozygous for non-refernce allele
# if there are 2 bases in reference genome (allele 1) but 1 in sample (allele 2), its an indel

genotype.matrix <- genotypes.tab[,sample.names]
head(genotype.matrix)

NA

Analyze genotype data

Here are a few questions to get people familiar with the genotype data

Q: How many samples are there?

n.samples <- length(sample.names)
n.samples

[1] 53

Q: what is the chromosomal position and alleles of SNP rs55926024?

snp.row <- genotypes.tab$SNP.ID == "rs55926024"
snp.row

   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [38] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
  [75] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [112] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [149] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [186] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [223] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [260] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [297] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [334] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [371] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [408] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [445] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [482] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [519] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [556] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [593] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [630] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [667] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [704] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [741] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [778] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [815] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [852] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [889] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [926] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [963] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[1000] FALSE
 [ reached getOption("max.print") -- omitted 263161 entries ]

Q: How many individuals are heterozygous for SNP rs55926024?

sum(genotype.matrix[snp.row,] == 1)

[1] 22

# just count number of columns where someone has a "1"

Q: What is the frequency of allele2 for rs55926024 in our samples?

total.alleles <- n.samples * 2
# multiply times 2 since looking at a diploid (2 allele copies)
sum(genotype.matrix[snp.row,]) / total.alleles

[1] 0.3584906

# since less frequent allele, sometimes called minor or (minor allele frequency = MAF (?))

Gene expression data

We will use RNA-seq gene expression data from the GEUVADIS project. The data are from the same Yoruba cell lines that we have genotypes for. The reads have already been mapped to the genome and gene expression has been estimated using the transcripts-per-million (TPM) measure. (For an explanation of TPM see http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/).

Read in mapped RNA-seq read counts for each gene. The first columns 5 give information about the gene. The subsequent columns that start with “NA” (e.g. NA18520) contain TPM estimates for each sample in the study.

# VCF = file that genetoying info is provided in. Has a lot of information in it


tpm.expr.tab <- read.table("gene_expr_tpm.txt", header=T)
head(tpm.expr.tab)

NA

Analyze gene expression values

Q: How many genes are in the gene expression table?

nrow(tpm.expr.tab)

[1] 63677

This is a lot of ‘genes’! The Ensembl ‘ENSG’ gene list contains many pseudo-genes, and low-expression lncRNAs etc.

Calculate mean expression across samples, then make a histogram of log mean TPM values for genes that have mean TPM > 0.0.

# so let's just get the ones that have the biggest expression
# take a mean acorss each row (gene) from the different cell lines
# extract just the expression data:
all.tpm.matrix <- tpm.expr.tab[,sample.names]

# compute mean expression and plog log mean TPM
# gene expression has a huge range, so if you want to look at that distrubution, better to use Log. So first filter out gene with 0 exprssion (we don't care about those)
mean.expr <- apply(all.tpm.matrix, 1, mean)
f <- mean.expr > 0.0
hist(log(mean.expr[f]), breaks=100)

# x axis = log TPM
# so things below 0 are probalby not expressed

The mean expression distribution is bimodal, with a large number of low-expression genes.

Q: How many genes have mean TPM values > 1.0?

# Take genes with expression value of TPM >1
sum(mean.expr > 1)

[1] 16279

This is a reasonable number of genes to consider expressed in this cell type (lymphoblastoid cell lines). So lets discard all genes with TPM < 1.0 from further analyses.

mean.expr <- apply(all.tpm.matrix, 1, mean)
f <- mean.expr > 1
hist(log(mean.expr[f]), breaks=100)

Other data pre-processing

There are additional steps we could take when processing the gene expression and genotype data.

• Correct for population structure by doing principal components analysis of genotype matrix (PCs could be used as co-variates in linear model)
• Correct for expression batch effects by regressing out PCs (principle components), performing quantile normalization

identifying eQTLs

Get a set of genes that are on chr22 and that have TPM > 1.0

f <- (tpm.expr.tab$CHROM == "chr22") & (mean.expr > 1.0)
chr22.tpm.expr.tab <- tpm.expr.tab[f,]
chr22.tpm.matrix <- tpm.expr.tab[f,sample.names] 

head(chr22.tpm.matrix)

# now every cell line has genes ewith expression values

Plot the expression of gene ENSG00000218537 (MIF-AS1) as a function of the genotype of the SNP rs5760108.

# gene expression for one gene
# gene ENSG, SNP rs###
expr <- unlist(chr22.tpm.matrix[chr22.tpm.expr.tab$GENE.NAME == "ENSG00000218537",])
genos <- unlist(genotype.matrix[genotypes.tab$SNP.ID == "rs5760108",])

plot(jitter(genos), expr, xlab="genotype", ylab="gene expression (TPM)")

# jitter adds noise so you can see more details of your plot "avoid overplotting"

Fit a linear model, using expression as the response and genotype as the predictor

mdl <- lm(expr ~ genos)
summary(mdl)

Call:
lm(formula = expr ~ genos)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.44888 -0.73025 -0.04898  0.66507  2.46593 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.1976     0.2372   0.833    0.409    
genos         1.2513     0.1841   6.797 1.14e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9117 on 51 degrees of freedom
Multiple R-squared:  0.4753,    Adjusted R-squared:  0.465 
F-statistic:  46.2 on 1 and 51 DF,  p-value: 1.138e-08

Get the slope (often called effect size or Beta) and p-value from the regression, and add the fit line to the genotype vs. expression plot.

intercept <- mdl$coefficients[1]
beta <- mdl$coefficients[2]
s <- summary(mdl)
p.val <- s$coefficients["genos","Pr(>|t|)"]

plot(jitter(genos), expr, xlab="genotype", ylab="gene expression (TPM)")
abline(a=intercept, b=beta, col="red")

Now lets test all SNPs that are close (within 50kb) of the gene ENSG00000218537 for association.

First get the genotypes of all of the nearby SNPs.

larger windows better for larger sample sizes

Why do we test SNPs only around the gene, and not throughout the whole genome?

We want to enrich for near the gene because that’s where the regulatory elements are

Whole genome would be nice, but that is very challenging. some people do it– and they find trans-acting eQTLs

• cis-acting eQTLs: only act on the same allele/chomosotme

• Trans-acting eQTLs: will act on same chromo and another chromosome

• Most variation and expression of genes is probably trans

# cis-alleles within 50kb of the gene to identify regulatory elements to test all the SNPs that are -50 and +50kb of the gene
gene.row <- which(chr22.tpm.expr.tab$GENE.NAME == "ENSG00000218537")
start <- chr22.tpm.expr.tab$GENE.START[gene.row] - 50000
end <- chr22.tpm.expr.tab$GENE.END[gene.row] + 50000

snp.rows <- which(genotypes.tab$CHROM.POS >= start & genotypes.tab$CHROM.POS <= end)

genotypes <- genotype.matrix[snp.rows,]
n.snps <- nrow(genotypes)

n.snps

[1] 1081

Next perform linear regression of the genoypes of each SNP against the expression of the gene and make a histogram of the p-values.

p.vals <- rep(NA, n.snps)

for(i in seq(1, n.snps)) {
  genos <- unlist(genotypes[i,])
  mdl <- lm(expr ~ genos)
  s <- summary(mdl)
  p.vals[i] <- s$coefficients["genos","Pr(>|t|)"]
}

hist(p.vals)

Plot the genomic position of the SNPs (x-axis) against their -log10 p-values (y-axis)

pos <- genotypes.tab$CHROM.POS[snp.rows]
plot(pos, -log10(p.vals), ylab="-log10(p-value)", xlab="chr22 position")
# add the position of the gene
lines(x=c(chr22.tpm.expr.tab$GENE.START[gene.row], chr22.tpm.expr.tab$GENE.END[gene.row]), y=c(0,0), 
         col="blue", lwd=4)