Differential gene expression analysis for RNA-seq after HTSeq-count from the paper "Epigenomic plasticity enables human pancreatic α to β cell reprogramming"

alpha_beta_DEG.r

  library(limma)
  library(edgeR)

  x<-read.delim('counts.csv',skip=0, sep="\t", check.names=FALSE)
  counts <- x[,c('a1','a2','a3','b1','b2','b3')]
  keep <- apply(counts, 1, max) >= 0
  x <- x[keep,]
  counts <- counts[keep,]
  design <- matrix(data=c(1,1,1,0,0,0,0,0,0,1,1,1), nrow=6, ncol=2, dimnames = list(c(), c('alpha','beta')))
 

  nf <- calcNormFactors(counts)
  y<-voom(counts, design, plot=FALSE,lib.size=colSums(counts)*nf)

  cont.matrix <- matrix(data=c(-1,1), nrow=2, ncol=1, dimnames = list(c(), c('beta')))

  fit <- lmFit(y,design)
  fit2 <- contrasts.fit(fit, cont.matrix)
  fit2 <- eBayes(fit2)

  out <- topTable(fit2, n=Inf, sort.by='none')

  out2 <- cbind(fit2$coef,
                out[, c('adj.P.Val','AveExpr')],
                x[, c('Gene','a1','a2','a3','b1','b2','b3')] )

  write.csv(out2, file="out-file.csv", row.names=FALSE,na='')