Differential Expression

Raw Counts

Fix header for fc files:

python fixExonFC.py /project2/gilad/briana/Comparative_APA/Human/data/RNAseq/ExonCounts/RNAseqOrthoExon.fc /project2/gilad/briana/Comparative_APA/Human/data/RNAseq/ExonCounts/RNAseqOrthoExon.fixed.fc

python fixExonFC.py /project2/gilad/briana/Comparative_APA/Chimp/data/RNAseq/ExonCounts/RNAseqOrthoExon.fc /project2/gilad/briana/Comparative_APA/Chimp/data/RNAseq/ExonCounts/RNAseqOrthoExon.fixed.fc

HumanCounts=read.table("../Human/data/RNAseq/ExonCounts/RNAseqOrthoExon.fixed.fc", header = T, stringsAsFactors = F) %>% select(-Chr,-Start,-End,-Strand, -Length)

ChimpCounts=read.table("../Chimp/data/RNAseq/ExonCounts/RNAseqOrthoExon.fixed.fc", header = T, stringsAsFactors = F) %>% select(-Chr,-Start,-End,-Strand, -Length)


counts_genes=HumanCounts %>% inner_join(ChimpCounts,by="Geneid") %>% column_to_rownames(var="Geneid")

head(counts_genes)

                NA18498 NA18504 NA18510 NA18523 NA18499 NA18502 NA4973
ENSG00000188976      48      24      50      31      34      58      2
ENSG00000188157     148     106      59     106     128      39      6
ENSG00000273443      69      40      43      54      48      11     24
ENSG00000217801      77      60      36     164      61      19     62
ENSG00000237330       1       0       1       0       1       1      0
ENSG00000223823       0       0       0       0       0       0      0
                NAPT30 NAPT91 NA3622 NA3659 NA18358
ENSG00000188976      2      1      1      1       0
ENSG00000188157      5      7      7      8       6
ENSG00000273443     21      2      3     78      16
ENSG00000217801     34      8     19    139      31
ENSG00000237330      0      0      0      2       0
ENSG00000223823      0      0      0      0       0

# Load colors 

colors <- colorRampPalette(c(brewer.pal(9, "Blues")[1],brewer.pal(9, "Blues")[9]))(100)

pal <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))
labels <- paste(metaData$Species,metaData$Line, sep=" ")

#PCA function (original code from Julien Roux)
#Load in the plot_scores function
plot_scores <- function(pca, scores, n, m, cols, points=F, pchs =20, legend=F){
  xmin <- min(scores[,n]) - (max(scores[,n]) - min(scores[,n]))*0.05
  if (legend == T){ ## let some room (35%) for a legend                                                                                                                                                 
    xmax <- max(scores[,n]) + (max(scores[,n]) - min(scores[,n]))*0.50
  }
  else {
    xmax <- max(scores[,n]) + (max(scores[,n]) - min(scores[,n]))*0.05
  }
  ymin <- min(scores[,m]) - (max(scores[,m]) - min(scores[,m]))*0.05
  ymax <- max(scores[,m]) + (max(scores[,m]) - min(scores[,m]))*0.05
  plot(scores[,n], scores[,m], xlab=paste("PC", n, ": ", round(summary(pca)$importance[2,n],3)*100, "% variance explained", sep=""), ylab=paste("PC", m, ": ", round(summary(pca)$importance[2,m],3)*100, "% variance explained", sep=""), xlim=c(xmin, xmax), ylim=c(ymin, ymax), type="n")
  if (points == F){
    text(scores[,n],scores[,m], rownames(scores), col=cols, cex=1)
  }
  else {
    points(scores[,n],scores[,m], col=cols, pch=pchs, cex=1.3)
  }
}

# Clustering (original code from Julien Roux)
cors <- cor(counts_genes, method="spearman", use="pairwise.complete.obs")


heatmap.2( cors, scale="none", col = colors, margins = c(12, 12), trace='none', denscol="white", labCol=labels, ColSideColors=pal[as.integer(as.factor(metaData$Species))], RowSideColors=pal[as.integer(as.factor(metaData$Collection))+9], cexCol = 0.2 + 1/log10(15), cexRow = 0.2 + 1/log10(15))

Past versions of unnamed-chunk-11-1.png

select <- counts_genes
summary(apply(select, 1, var) == 0) 

   Mode   FALSE    TRUE 
logical   32602   11523 

# Perform PCA

pca_genes <- prcomp(t(counts_genes), scale = F)
scores <- pca_genes$x


#Make PCA plots with the factors colored by species

### PCs 1 and 2 Raw Data
for (n in 1:1){
  col.v <- pal[as.integer(metaData$Species)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-13-1.png

### PCs 3 and 4 Raw Data

for (n in 3:3){
  col.v <- pal[as.integer(metaData$Species)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-14-1.png

Plot density for raw data:

density_plot_18504 <- ggplot(counts_genes, aes(x = NA18504)) + geom_density() + labs(title = "Density plot of raw gene counts of NA18504") + labs(x = "Raw counts for each gene")

density_plot_18504

Past versions of unnamed-chunk-15-1.png

Convert to log2

log_counts_genes <- as.data.frame(log2(counts_genes))
head(log_counts_genes)

                 NA18498  NA18504  NA18510  NA18523  NA18499  NA18502
ENSG00000188976 5.584963 4.584963 5.643856 4.954196 5.087463 5.857981
ENSG00000188157 7.209453 6.727920 5.882643 6.727920 7.000000 5.285402
ENSG00000273443 6.108524 5.321928 5.426265 5.754888 5.584963 3.459432
ENSG00000217801 6.266787 5.906891 5.169925 7.357552 5.930737 4.247928
ENSG00000237330 0.000000     -Inf 0.000000     -Inf 0.000000 0.000000
ENSG00000223823     -Inf     -Inf     -Inf     -Inf     -Inf     -Inf
                  NA4973   NAPT30   NAPT91   NA3622   NA3659  NA18358
ENSG00000188976 1.000000 1.000000 0.000000 0.000000 0.000000     -Inf
ENSG00000188157 2.584963 2.321928 2.807355 2.807355 3.000000 2.584963
ENSG00000273443 4.584963 4.392317 1.000000 1.584963 6.285402 4.000000
ENSG00000217801 5.954196 5.087463 3.000000 4.247928 7.118941 4.954196
ENSG00000237330     -Inf     -Inf     -Inf     -Inf 1.000000     -Inf
ENSG00000223823     -Inf     -Inf     -Inf     -Inf     -Inf     -Inf

density_plot_18504 <- ggplot(log_counts_genes, aes(x = 18504)) + geom_density()

density_plot_18504 + labs(title = "Density plot of log2 counts of 18504") + labs(x = "Log2 counts for each gene") + geom_vline(xintercept = 1) 

Past versions of unnamed-chunk-17-1.png

plotDensities(log_counts_genes, col=pal[as.numeric(metaData$Species)], legend="topright")

Past versions of unnamed-chunk-18-1.png

Convert to CPM

cpm <- cpm(counts_genes, log=TRUE)
head(cpm)

                  NA18498    NA18504   NA18510   NA18523   NA18499
ENSG00000188976  1.588860  0.9875607  1.619314  1.208048  1.067299
ENSG00000188157  3.172608  3.0589988  1.849193  2.925931  2.915757
ENSG00000273443  2.094212  1.6879881  1.411168  1.975384  1.539689
ENSG00000217801  2.248116  2.2543286  1.167845  3.547215  1.872294
ENSG00000237330 -2.441797 -3.0119624 -2.450603 -3.011962 -2.458352
ENSG00000223823 -3.011962 -3.0119624 -3.011962 -3.011962 -3.011962
                   NA18502     NA4973     NAPT30     NAPT91     NA3622
ENSG00000188976  1.8417566 -2.1228041 -2.0245721 -2.4061474 -2.4360717
ENSG00000188157  1.2932833 -1.1815912 -1.2226553 -0.7938422 -0.8632943
ENSG00000273443 -0.3579576  0.4766994  0.4884893 -1.9807956 -1.7064048
ENSG00000217801  0.3306479  1.7649286  1.1340509 -0.6404855  0.3555977
ENSG00000237330 -2.4451144 -3.0119624 -3.0119624 -3.0119624 -3.0119624
ENSG00000223823 -3.0119624 -3.0119624 -3.0119624 -3.0119624 -3.0119624
                    NA3659    NA18358
ENSG00000188976 -2.4830271 -3.0119624
ENSG00000188157 -0.8283452 -0.8093695
ENSG00000273443  2.1395468  0.3951397
ENSG00000217801  2.9551621  1.2819822
ENSG00000237330 -2.0968409 -3.0119624
ENSG00000223823 -3.0119624 -3.0119624

plotDensities(cpm, col=pal[as.numeric(metaData$Species)], legend="topright")

Past versions of unnamed-chunk-20-1.png

Log2 CPM

TMM/log2(CPM)

## Create edgeR object (dge) to calculate TMM normalization  
dge_original <- DGEList(counts=as.matrix(counts_genes), genes=rownames(counts_genes), group = as.character(t(labels)))
dge_original <- calcNormFactors(dge_original)

tmm_cpm <- cpm(dge_original, normalized.lib.sizes=TRUE, log=TRUE, prior.count = 0.25)
head(cpm)

                  NA18498    NA18504   NA18510   NA18523   NA18499
ENSG00000188976  1.588860  0.9875607  1.619314  1.208048  1.067299
ENSG00000188157  3.172608  3.0589988  1.849193  2.925931  2.915757
ENSG00000273443  2.094212  1.6879881  1.411168  1.975384  1.539689
ENSG00000217801  2.248116  2.2543286  1.167845  3.547215  1.872294
ENSG00000237330 -2.441797 -3.0119624 -2.450603 -3.011962 -2.458352
ENSG00000223823 -3.011962 -3.0119624 -3.011962 -3.011962 -3.011962
                   NA18502     NA4973     NAPT30     NAPT91     NA3622
ENSG00000188976  1.8417566 -2.1228041 -2.0245721 -2.4061474 -2.4360717
ENSG00000188157  1.2932833 -1.1815912 -1.2226553 -0.7938422 -0.8632943
ENSG00000273443 -0.3579576  0.4766994  0.4884893 -1.9807956 -1.7064048
ENSG00000217801  0.3306479  1.7649286  1.1340509 -0.6404855  0.3555977
ENSG00000237330 -2.4451144 -3.0119624 -3.0119624 -3.0119624 -3.0119624
ENSG00000223823 -3.0119624 -3.0119624 -3.0119624 -3.0119624 -3.0119624
                    NA3659    NA18358
ENSG00000188976 -2.4830271 -3.0119624
ENSG00000188157 -0.8283452 -0.8093695
ENSG00000273443  2.1395468  0.3951397
ENSG00000217801  2.9551621  1.2819822
ENSG00000237330 -2.0968409 -3.0119624
ENSG00000223823 -3.0119624 -3.0119624

pca_genes <- prcomp(t(tmm_cpm), scale = F)
scores <- pca_genes$x

for (n in 1:2){
  col.v <- pal[as.integer(metaData$Species)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-22-1.png

Past versions of unnamed-chunk-22-2.png

# Plot library size

boxplot_library_size <- ggplot(dge_original$samples, aes(x=metaData$Species, y = dge_original$samples$lib.size, fill = metaData$Species)) + geom_boxplot()
 
boxplot_library_size + labs(title = "Library size by Species") + labs(y = "Library size") + labs(x = "Species") + guides(fill=guide_legend(title="Species"))

Past versions of unnamed-chunk-23-1.png

plotDensities(tmm_cpm, col=pal[as.numeric(metaData$Species)], legend="topright")

Past versions of unnamed-chunk-24-1.png

Filter low expressed gene

Filter based on log2 cpm

filter log2(cpm >1) in at least 10 of the samples (2/3)

#filter counts
keep.exprs=rowSums(tmm_cpm>1) >8

counts_filtered= counts_genes[keep.exprs,]




plotDensities(counts_filtered, col=pal[as.numeric(metaData$Species)], legend="topright")

Past versions of unnamed-chunk-25-1.png

labels <- paste(metaData$Species, metaData$Line, sep=" ")
dge_in_cutoff <- DGEList(counts=as.matrix(counts_filtered), genes=rownames(counts_filtered), group = as.character(t(labels)))
dge_in_cutoff <- calcNormFactors(dge_in_cutoff)

cpm_in_cutoff <- cpm(dge_in_cutoff, normalized.lib.sizes=TRUE, log=TRUE, prior.count = 0.25)
head(cpm_in_cutoff)

                 NA18498  NA18504    NA18510  NA18523  NA18499   NA18502
ENSG00000217801 2.256718 2.234174  1.0782676 3.660677 1.851409 0.2045780
ENSG00000186891 6.449880 5.009564  5.0298795 4.887036 5.699309 3.3564937
ENSG00000186827 5.245424 3.093249  4.6497808 1.700761 3.965562 1.8968145
ENSG00000078808 7.201366 6.859761  7.0261290 7.463811 6.617502 6.7318204
ENSG00000176022 4.796475 4.753656  4.6762970 4.789968 4.605179 4.6954931
ENSG00000184163 1.692704 1.418512 -0.1694645 1.630729 1.054723 0.3470491
                  NA4973    NAPT30     NAPT91   NA3622    NA3659  NA18358
ENSG00000217801 1.778595 1.1604206 -0.7907848 0.275830 2.9986515 1.300800
ENSG00000186891 6.668817 4.6333309  3.1597896 4.727162 6.4478748 6.482073
ENSG00000186827 1.778595 0.2520802  2.0071797 3.150776 0.6465898 2.175191
ENSG00000078808 7.024203 6.6875607  6.7658411 6.776178 6.7839057 7.060989
ENSG00000176022 5.263240 4.8315395  5.5505445 4.845467 5.2761471 5.561648
ENSG00000184163 1.381145 2.8985450  2.1959470 1.409687 1.7662843 3.222286

hist(cpm_in_cutoff, xlab = "Log2(CPM)", main = "Log2(CPM) values for genes meeting the filtering criteria", breaks = 100 )

Past versions of unnamed-chunk-27-1.png

Voom transformation:

Species <- factor(metaData$Species)
design <- model.matrix(~ 0 + Species)
head(design)

  SpeciesChimp SpeciesHuman
1            1            0
2            1            0
3            1            0
4            1            0
5            1            0
6            1            0

colnames(design) <- gsub("Species", "", dput(colnames(design)))

c("SpeciesChimp", "SpeciesHuman")

Voom creates a random effect.

# Voom with individual as a random variable

cpm.voom<- voom(counts_filtered, design, normalize.method="quantile", plot=T)

Past versions of unnamed-chunk-29-1.png

boxplot(cpm.voom$E, col = pal[as.numeric(metaData$Species)],las=2)

Past versions of unnamed-chunk-30-1.png

plotDensities(cpm.voom, col =  pal[as.numeric(metaData$Species)], legend = "topleft") 

Past versions of unnamed-chunk-31-1.png

Looks like i still have a skew on the lower side of the distribution.

# PCA 

pca_genes <- prcomp(t(cpm.voom$E), scale = T)
scores <- pca_genes$x


eigsGene <- pca_genes$sdev^2
proportionG = eigsGene/sum(eigsGene)

plot(proportionG)

Past versions of unnamed-chunk-32-1.png

for (n in 1:2){
  col.v <- pal[as.integer(metaData$Species)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-32-2.png

Past versions of unnamed-chunk-32-3.png

#Clustering (original code from Julien Roux)
cors <- cor(cpm.voom$E, method="spearman", use="pairwise.complete.obs")


heatmap.2( cors, scale="none", col = colors, margins = c(12, 12), trace='none', denscol="white", labCol=labels, ColSideColors=pal[as.integer(as.factor(metaData$Species))], RowSideColors=pal[as.integer(as.factor(metaData$Species))+9], cexCol = 0.2 + 1/log10(15), cexRow = 0.2 + 1/log10(15))

Past versions of unnamed-chunk-33-1.png

One thing I can do is look at the correlation between the PCs and other factors in the data.

# PCA 

pca_genes <- prcomp(t(cpm.voom$E), scale = F)
scores <- pca_genes$x

for (n in 1:2){
  col.v <- pal[as.integer(metaData$Collection)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-34-1.png

Past versions of unnamed-chunk-34-2.png

metaData$Extraction=as.factor(metaData$Extraction)

for (n in 1:2){
  col.v <- pal[as.integer(metaData$Extraction)]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-35-1.png

Past versions of unnamed-chunk-35-2.png

It does not look like batch (who collected or extraction date batch)

cols = brewer.pal(9, "Blues")
palC = colorRampPalette(cols)


metaData$UndilutedAverageorder = findInterval(metaData$UndilutedAverage, sort(metaData$UndilutedAverage))
for (n in 1:2){
  col.v <- palC(nrow(metaData))[metaData$UndilutedAverageorder]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-36-1.png

Past versions of unnamed-chunk-36-2.png

metaData$BioAConcorder = findInterval(metaData$BioAConc, sort(metaData$BioAConc))
for (n in 1:2){
  col.v <- palC(nrow(metaData))[metaData$BioAConcorder]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-37-1.png

Past versions of unnamed-chunk-37-2.png

metaData$RinConcorder = findInterval(metaData$Rin, sort(metaData$Rin))
for (n in 1:2){
  col.v <- palC(nrow(metaData))[metaData$RinConcorder]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-38-1.png

Past versions of unnamed-chunk-38-2.png

The samples do not cluster by collection concentration, RNA rin score or RNA concentration.

metaData$AssignedOrthoorder = findInterval(metaData$AssignedOrtho, sort(metaData$AssignedOrtho))
for (n in 1:2){
  col.v <- palC(nrow(metaData))[metaData$AssignedOrthoorder]
  plot_scores(pca_genes, scores, n, n+1, col.v)
}

Past versions of unnamed-chunk-39-1.png

Past versions of unnamed-chunk-39-2.png

They also do not cluster by number of reads mapping to ortho exons.

PCA heatmap: Code from Michelle Ward:

x.pca <- pca_genes

tech_factors <- metaData
tech_factors_sum <- tech_factors[,c(2:15)] %>% select(-CollectionDate)

p_comps <- 1:6
pc_cov_cor <- matrix(nrow = ncol(tech_factors_sum), ncol = length(p_comps),
                     dimnames = list(colnames(tech_factors_sum), colnames(x.pca$x)[p_comps]))
for (pc in p_comps) {
  for (covariate in 1:ncol(tech_factors_sum)) {
    lm_result <- lm(x.pca$x[, pc] ~ tech_factors_sum[, covariate])
    r2 <- summary(lm_result)$r.squared
    pc_cov_cor[covariate, pc] <- r2
  }
}

pc_cov_pval <- matrix(nrow = ncol(tech_factors_sum), ncol = length(p_comps),
                      dimnames = list(colnames(tech_factors_sum), colnames(x.pca$x)[p_comps]))

for (pc in p_comps) {
  for (covariate_2 in 1:ncol(tech_factors_sum)) {
    lm_result_2 <- lm(x.pca$x[, pc] ~ tech_factors_sum[, covariate_2])
    pval <- anova(lm_result_2)$'Pr(>F)'[1]
    pc_cov_pval[covariate_2, pc] <- pval
  }
}

PCs <- c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6")
Tech_fac <- colnames(tech_factors_sum)
#Tech_fac <- c("Species",   "Individual", "O2.",  "Condition" , "Sex", "RIN" , "CO2", "Purity_high", "Purity_med" ,
              #"Expt_Batch", "RNA_Batch", "Library_Batch", "Seq_pool", "Episomal_integration" )

heatmap.2(as.matrix(pc_cov_cor[Tech_fac,PCs]),col=brewer.pal(4, "Greens"), trace="none",
          Rowv=FALSE, Colv=FALSE, key=T, main="Cor. PCs & tech factors", dendrogram="none",
          key.title=NA, cexRow=0.9, cexCol=0.9)

Past versions of unnamed-chunk-40-1.png

log10_pc_cov_pval <- -log(pc_cov_pval)
heatmap.2(as.matrix(log10_pc_cov_pval[Tech_fac,PCs]), col=brewer.pal(9, "Greens"), trace="none",
          Rowv=FALSE, Colv=FALSE, key=T, main="-log10 pval of cor. PCs & tech factors", dendrogram="none",
          key.title=NA, cexRow=0.9, cexCol=0.9)

Past versions of unnamed-chunk-40-2.png