NPC_vs_ipsNeuron_expression.R

# Gabriel Hoffman
# September 16, 2020
#
# Find expression of genes in iPS-neurons and NPCs

library(edgeR)
library(limma)
library(ggplot2)
library(reshape2)
library(synapser)

synLogin()

geneExpr_rpkm = read.table(synGet('syn6613774')$path, header=TRUE, row.names=1, sep=',', check.names=FALSE)

info = read.table(synGet('syn7079773')$path, header=TRUE, row.names=1, sep=',')

geneInfo = read.table(synGet('syn7113731')$path, header=TRUE, row.names=1, sep=',')


# Calculate TPM from RPKM
tpm_from_rpkm <- function(x) {
  rpkm.sum <- colSums(x)
  return(t(t(x) / (1e-06 * rpkm.sum)))
}

geneExpr_tpm = tpm_from_rpkm( geneExpr_rpkm)

genes = c('RERE', 'RPS21', 'EEF1A2')

ensID = rownames(geneInfo)[match(genes, geneInfo$geneName)]
names(ensID) = genes


idx = which(info$cellType == 'F')
df_FB = data.frame( t(geneExpr_tpm[ensID,idx]), cellType = "FB neuron")
colnames(df_FB)[match(ensID, colnames(df_FB))] = names(ensID)
d_FB = density(log2(geneExpr_tpm[,idx]))


idx = which(info$cellType == 'N')
df_NPC = data.frame( t(geneExpr_tpm[ensID,idx]), cellType = "NPC")
colnames(df_NPC)[match(ensID, colnames(df_NPC))] = names(ensID)
d_NPC = density(log2(geneExpr_tpm[,idx]))

# plot(d_FB, col="orange")
# lines(d_NPC, col="green")


df = rbind(melt(df_FB), melt(df_NPC))

ggplot(df, aes(variable, log2(value))) + geom_point() + theme_bw() + facet_wrap(~cellType) + theme(aspect.ratio=1) + ylab(bquote(log[2]~transcripts~per~million)) + xlab("")