Differential Expression Remove Bad Individuals

Raw Counts

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

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


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

head(counts_genes)

                NA18504 NA18510 NA18523 NA18499 NA18502 NAPT30 NAPT91
ENSG00000188976      24      50      31      34      58      2      1
ENSG00000188157     106      65     106     128      39      5      7
ENSG00000273443      40      43      54      48      11     21      2
ENSG00000217801      60      36     164      61      19     34      8
ENSG00000237330       0       1       0       1       1      0      0
ENSG00000223823       0       0       0       0       0      0      0
                NA3622 NA3659 NA18358
ENSG00000188976      1      1       0
ENSG00000188157      7      8       6
ENSG00000273443      3     78      16
ENSG00000217801     19    139      31
ENSG00000237330      0      2       0
ENSG00000223823      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-10-1.png

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

   Mode   FALSE    TRUE 
logical   31883   12242 

# 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-12-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-13-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-14-1.png

Convert to log2

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

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

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

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

Convert to CPM

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

                   NA18504   NA18510   NA18523   NA18499    NA18502
ENSG00000188976  0.9894703  1.620342  1.209959  1.068977  1.8425919
ENSG00000188157  3.0594442  1.985078  2.926705  2.916115  1.2946137
ENSG00000273443  1.6891590  1.412390  1.976621  1.540857 -0.3532791
ENSG00000217801  2.2551156  1.169332  3.547816  1.873190  0.3334631
ENSG00000237330 -2.9815417 -2.429988 -2.981542 -2.437603 -2.4245831
ENSG00000223823 -2.9815417 -2.981542 -2.981542 -2.981542 -2.9815417
                   NAPT30     NAPT91     NA3622     NA3659    NA18358
ENSG00000188976 -2.008556 -2.3855698 -2.4153553 -2.4615104 -2.9815417
ENSG00000188157 -1.212936 -0.7860576 -0.8558140 -0.8206567 -0.8020982
ENSG00000273443  0.492301 -1.9650589 -1.6935886  2.1415838  0.3987400
ENSG00000217801  1.136923 -0.6333308  0.3592322  2.9568407  1.2842881
ENSG00000237330 -2.981542 -2.9815417 -2.9815417 -2.0800683 -2.9815417
ENSG00000223823 -2.981542 -2.9815417 -2.9815417 -2.9815417 -2.9815417

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

Past versions of unnamed-chunk-19-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   NA18499    NA18502
ENSG00000188976  0.9894703  1.620342  1.209959  1.068977  1.8425919
ENSG00000188157  3.0594442  1.985078  2.926705  2.916115  1.2946137
ENSG00000273443  1.6891590  1.412390  1.976621  1.540857 -0.3532791
ENSG00000217801  2.2551156  1.169332  3.547816  1.873190  0.3334631
ENSG00000237330 -2.9815417 -2.429988 -2.981542 -2.437603 -2.4245831
ENSG00000223823 -2.9815417 -2.981542 -2.981542 -2.981542 -2.9815417
                   NAPT30     NAPT91     NA3622     NA3659    NA18358
ENSG00000188976 -2.008556 -2.3855698 -2.4153553 -2.4615104 -2.9815417
ENSG00000188157 -1.212936 -0.7860576 -0.8558140 -0.8206567 -0.8020982
ENSG00000273443  0.492301 -1.9650589 -1.6935886  2.1415838  0.3987400
ENSG00000217801  1.136923 -0.6333308  0.3592322  2.9568407  1.2842881
ENSG00000237330 -2.981542 -2.9815417 -2.9815417 -2.0800683 -2.9815417
ENSG00000223823 -2.981542 -2.9815417 -2.9815417 -2.9815417 -2.9815417

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-21-1.png

Past versions of unnamed-chunk-21-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-22-1.png

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

Past versions of unnamed-chunk-23-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-24-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  NA18499  NA18502   NAPT30
ENSG00000186891 5.019760 5.041382 4.900641 5.711124 3.366537 4.651440
ENSG00000078808 6.869945 7.037620 7.477403 6.629313 6.741815 6.705657
ENSG00000176022 4.763855 4.687804 4.803574 4.617006 4.705504 4.849646
ENSG00000160087 4.849052 5.381048 5.119697 4.965717 5.086619 5.207984
ENSG00000131584 6.069280 4.920816 4.539908 4.737501 5.356339 5.133882
ENSG00000169972 3.357225 3.394890 3.555318 3.773273 3.315091 3.621735
                  NAPT91   NA3622   NA3659  NA18358
ENSG00000186891 3.175319 4.739154 6.460823 6.498061
ENSG00000078808 6.781323 6.788155 6.796852 7.076976
ENSG00000176022 5.566032 4.857458 5.289102 5.577640
ENSG00000160087 5.451318 5.542696 5.337421 5.414395
ENSG00000131584 5.709962 5.380906 5.258878 5.565531
ENSG00000169972 3.552346 4.077282 3.072501 3.321316

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

Past versions of unnamed-chunk-26-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            0            1

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-28-1.png

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

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

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

Past versions of unnamed-chunk-30-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-31-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-31-2.png

Past versions of unnamed-chunk-31-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-32-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-33-1.png

Past versions of unnamed-chunk-33-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-34-1.png

Past versions of unnamed-chunk-34-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-35-1.png

Past versions of unnamed-chunk-35-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-36-1.png

Past versions of unnamed-chunk-36-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-37-1.png

Past versions of unnamed-chunk-37-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-38-1.png

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

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