Go Back   SEQanswers > Bioinformatics > Bioinformatics

Similar Threads
Thread Thread Starter Forum Replies Last Post
differential gene expression nicole_01 RNA Sequencing 0 05-12-2014 05:04 PM
differential gene expression - Illumina nicole_01 Bioinformatics 1 10-29-2013 09:06 PM
Expression quantification/differential expression gene analysis by RNA-Seq chenjy Bioinformatics 12 08-02-2013 04:06 AM
Differential gene expression analysis colaneri Bioinformatics 15 06-14-2013 06:37 AM
Differential gene expression of gene clusters anjana.vr RNA Sequencing 1 10-28-2010 11:33 AM

Thread Tools
Old 10-20-2018, 01:14 PM   #1
Junior Member
Location: Virginia

Join Date: Aug 2018
Posts: 3
Default R Code Differential Gene Expression

Hi All,

I was wondering if you could look over my R code for differential gene expression using EdjeR. I am looking to determine differential gene expression between wild type (WT) cells and knockout cells (KO). Three biological replicates were grown for each cell line and RNA was harvested. The paired end reads were mapped using STAR. Exon counts were obtained using feature counts. The exon counts were then used for the R code. Would this be sufficient to determine differential gene expression between WT and KO?

library("edgeR") library("gdata") library("heatmaply") library("ggplot2") library("genefilter") library("methylumi")

setwd("/Users/jwlandry2/Desktop/RNA-Seq\ ESC\ Data\ Analysis/R_Studio_Data_Sets")

counts=read.delim("BPTFKD_ESC_RNA_Seq_Counts_Final2.txt", header=T, row.names="Geneid")


group <- factor(c("WT","WT","WT","KO","KO","KO"))

dge = DGEList(counts=counts,group=group)

keep <- rowSums(cpm(dge)>2) >= 3 dge <- dge[keep, , keep.lib.sizes=FALSE]

y <- calcNormFactors(dge) design <- model.matrix(~group) y <- estimateDisp(y,design) y$common.dispersion

plotMDS(y) plotBCV(y)

quasi-likelihood F-test
fit <- glmQLFit(y, design) qlf <- glmQLFTest(fit, coef=2) topTags(qlf)

output to text file
df <- qlf$table write.csv(df, 'qlf.csv')

Get normalizaed CPMs
mtx <- cpm(y, log = TRUE, normalized.lib.sizes = TRUE) mtx_to_plot <- varFilter(mtx, var.cutoff= 0.95)

Correlation matrix
IAC <- mtx_to_plot %>% cor(. , use = "pairwise.complete.obs", method = "pearson") heatmaply(IAC) %>% layout(margin=list(l=200,b=150))

plot(hclust(as.dist(1-IAC), method="ward"))
landrjos is offline   Reply With Quote

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

All times are GMT -8. The time now is 04:06 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO