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)

                NA18504 NA18510 NA18523 NA18498 NA18499 NA18502 NAPT30
ENSG00000188976      24      50      31      69      34      58      2
ENSG00000188157     106      65     106      12     128      39      5
ENSG00000273443      40      43      54      26      48      11     21
ENSG00000217801      60      36     164      62      61      19     34
ENSG00000237330       0       1       0       0       1       1      0
ENSG00000223823       0       0       0       0       0       0      0
                NAPT91 NA3622 NA3659 NA4973 NA18358
ENSG00000188976      1      1      1      0       0
ENSG00000188157      7      7      8     27       6
ENSG00000273443      2      3     78     59      16
ENSG00000217801      8     19    139     68      31
ENSG00000237330      0      0      2      1       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   32864   11261 

# 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)

                 NA18504  NA18510  NA18523  NA18498  NA18499  NA18502
ENSG00000188976 4.584963 5.643856 4.954196 6.108524 5.087463 5.857981
ENSG00000188157 6.727920 6.022368 6.727920 3.584963 7.000000 5.285402
ENSG00000273443 5.321928 5.426265 5.754888 4.700440 5.584963 3.459432
ENSG00000217801 5.906891 5.169925 7.357552 5.954196 5.930737 4.247928
ENSG00000237330     -Inf 0.000000     -Inf     -Inf 0.000000 0.000000
ENSG00000223823     -Inf     -Inf     -Inf     -Inf     -Inf     -Inf
                  NAPT30   NAPT91   NA3622   NA3659   NA4973  NA18358
ENSG00000188976 1.000000 0.000000 0.000000 0.000000     -Inf     -Inf
ENSG00000188157 2.321928 2.807355 2.807355 3.000000 4.754888 2.584963
ENSG00000273443 4.392317 1.000000 1.584963 6.285402 5.882643 4.000000
ENSG00000217801 5.087463 3.000000 4.247928 7.118941 6.087463 4.954196
ENSG00000237330     -Inf     -Inf     -Inf 1.000000 0.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)

                   NA18504   NA18510   NA18523    NA18498   NA18499
ENSG00000188976  0.9877925  1.619259  1.208519  1.9826994  1.067390
ENSG00000188157  3.0590448  1.984237  2.926267 -0.3399951  2.915674
ENSG00000273443  1.6881262  1.411139  1.975775  0.6477686  1.539712
ENSG00000217801  2.2544180  1.167851  3.547531  1.8334849  1.872281
ENSG00000237330 -3.0080795 -2.448042 -3.008080 -3.0080795 -2.455753
ENSG00000223823 -3.0080795 -3.008080 -3.008080 -3.0080795 -3.008080
                   NA18502     NAPT30     NAPT91     NA3622     NA3659
ENSG00000188976  1.8416633 -2.0220152 -2.4030724 -2.4332254 -2.4799651
ENSG00000188157  1.2932560 -1.2206741 -0.7918098 -0.8618518 -0.8265488
ENSG00000273443 -0.3575383  0.4899323 -1.9781161 -1.7043978  2.1408291
ENSG00000217801  0.3308185  1.1354086 -0.6385042  0.3566344  2.9564118
ENSG00000237330 -2.4425685 -3.0080795 -3.0080795 -3.0080795 -2.0942143
ENSG00000223823 -3.0080795 -3.0080795 -3.0080795 -3.0080795 -3.0080795
                    NA4973    NA18358
ENSG00000188976 -3.0080795 -3.0080795
ENSG00000188157  0.7692807 -0.8079148
ENSG00000273443  1.8388114  0.3962125
ENSG00000217801  2.0369808  1.2829205
ENSG00000237330 -2.4514663 -3.0080795
ENSG00000223823 -3.0080795 -3.0080795

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)

                   NA18504   NA18510   NA18523    NA18498   NA18499
ENSG00000188976  0.9877925  1.619259  1.208519  1.9826994  1.067390
ENSG00000188157  3.0590448  1.984237  2.926267 -0.3399951  2.915674
ENSG00000273443  1.6881262  1.411139  1.975775  0.6477686  1.539712
ENSG00000217801  2.2544180  1.167851  3.547531  1.8334849  1.872281
ENSG00000237330 -3.0080795 -2.448042 -3.008080 -3.0080795 -2.455753
ENSG00000223823 -3.0080795 -3.008080 -3.008080 -3.0080795 -3.008080
                   NA18502     NAPT30     NAPT91     NA3622     NA3659
ENSG00000188976  1.8416633 -2.0220152 -2.4030724 -2.4332254 -2.4799651
ENSG00000188157  1.2932560 -1.2206741 -0.7918098 -0.8618518 -0.8265488
ENSG00000273443 -0.3575383  0.4899323 -1.9781161 -1.7043978  2.1408291
ENSG00000217801  0.3308185  1.1354086 -0.6385042  0.3566344  2.9564118
ENSG00000237330 -2.4425685 -3.0080795 -3.0080795 -3.0080795 -2.0942143
ENSG00000223823 -3.0080795 -3.0080795 -3.0080795 -3.0080795 -3.0080795
                    NA4973    NA18358
ENSG00000188976 -3.0080795 -3.0080795
ENSG00000188157  0.7692807 -0.8079148
ENSG00000273443  1.8388114  0.3962125
ENSG00000217801  2.0369808  1.2829205
ENSG00000237330 -2.4514663 -3.0080795
ENSG00000223823 -3.0080795 -3.0080795

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)

                 NA18504    NA18510  NA18523  NA18498  NA18499   NA18502
ENSG00000217801 2.248146  1.0811901 3.657635 1.792969 1.843394 0.2076129
ENSG00000186891 5.023562  5.0327842 4.883988 6.663573 5.691238 3.3594991
ENSG00000186827 3.107235  4.6526859 1.697751 1.815945 3.957501 1.8998264
ENSG00000078808 6.873762  7.0290328 7.460760 7.026971 6.609430 6.7348224
ENSG00000176022 4.767653  4.6792021 4.786920 5.258909 4.597113 4.6984962
ENSG00000184163 1.432462 -0.1665159 1.627720 1.395504 1.046753 0.3500809
                  NAPT30     NAPT91   NA3622   NA3659   NA4973  NA18358
ENSG00000217801 1.164822 -0.7903107 0.279009 3.001252 2.105148 1.303916
ENSG00000186891 4.637726  3.1601319 4.730313 6.450470 6.482165 6.485174
ENSG00000186827 0.256488  2.0075331 3.153929 0.649213 5.258424 2.178300
ENSG00000078808 6.691955  6.7661750 6.779328 6.786500 7.215767 7.064091
ENSG00000176022 4.835934  5.5508794 4.848618 5.278743 4.794461 5.564750
ENSG00000184163 2.902942  2.1962979 1.412850 1.768892 1.691914 3.225391

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

This is wierd. Normalization moves 2 samples to opposite species clusters but the samples that separate in the correlation are not those samples. 4973 and 18498 are the samples that looked funny on the original 3’ data. This may be a sample swap at the RNA stage. These samples were in the same extraction batch. It could have happened then. I will look into this more.

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:14)] %>% 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)