0
$\begingroup$

I was following the RNAseq analysis tutorial online(https://combine-australia.github.io/RNAseq-R/06-rnaseq-day1.html) but am not obtaining the same variance values for each row in the logcounts matrix. Please could you help me identify the error.

library(edgeR)
library(limma)
library(Glimma)
library(org.Mm.eg.db)
library(RColorBrewer)
library(gplots)

#Read the data into R
seqdata <- read.delim("data/GSE60450_Lactation-GenewiseCounts.txt", stringsAsFactors = FALSE)
sampleinfo <- read.delim("data/SampleInfo.txt")
head(seqdata)
dim(seqdata)
sampleinfo

##Format the data
countdata <- seqdata[,-(1:2)]
head(countdata)
rownames(countdata) <- seqdata[,1]
head(countdata)
colnames(countdata)
colnames
colnames(countdata) <- substr(colnames(countdata),start=1,stop=7)
head(countdata)
sampleinfo
table(colnames(countdata)==sampleinfo$SampleName)

##Filtering to remove lowly expressed genes
#Obtain CPMs
myCPM <- cpm(countdata)
head(myCPM)
thresh <- myCPM > 0.5
head(thresh)
table(rowSums(thresh))
keep <- rowSums(thresh) >=2
#Subset the rows of countdata to keep the more highly expressed genes 
counts.keep <- countdata[keep,]
summary(keep)
dim(counts.keep)
#Find out whether CPM threshold corresponds to a count of around 10-15
#Look at first sample
plot(myCPM[,1],countdata[,1],ylim=c(0,50),xlim=c(0,3))
abline(v=0.5)


#Convert counts to edgeList objects
y <- DGEList(counts.keep)
y
names(y)
y$samples

###Quality Control

#Check the number of reads we have for each sample in y
y$samples$lib.size
#Barplot of library size to spot any major discrepancies between the samples
barplot(y$samples$lib.size, names=colnames(y),las=2)
title("Barplot of Library Sizes")
#Get log2 CPM
logcounts <- cpm(y,log=TRUE)
#Check distribution of samples using boxplots
boxplot(logcounts, xlab="",ylab="log2 CPM", las=2, ylim=c(-10,15))
#Add a horizontal line corresponding to the median logCPM
abline(h=median(logcounts), col="blue")
title("Boxplots of logCPMs (unnormalised)")

##Heirarchical clustering with heatmaps
#Estimate the variance for each row in the logcounts matrix
var_genes <- apply(logcounts,1,var)
head(var_genes)

My incorrect variance values of the first 6 genes:497097,20671,27395 18777, 21399 and 58175 are 13.6624115, 2.7493077, 0.1581944, 0.1306781, 0.3929526 and 4.8232522 respectively.

My boxplot

My boxplot

Tutorial boxplot

Tutorial boxplot

$\endgroup$
  • $\begingroup$ Do you know if you are using the same version of the R modules as in the tutorial? $\endgroup$ – Wouter De Coster Jan 17 at 16:14
  • $\begingroup$ I can't find any information on which modules the tutorial uses. I am relatively new to R programming and haven't come across modules yet in R, but how would you find out what R modules you have? $\endgroup$ – Bio314 Jan 17 at 18:25
  • $\begingroup$ Do the rest of your outputs and graphs before that point match the tutorial? That tutorial is 3 years old, maybe a module has slightly changed how it calculates something; your values aren't that far off from the tutorial. And that log cpm business is for visualization only anyway. You use raw counts with EdgeR to find DE genes. $\endgroup$ – swbarnes2 Jan 17 at 19:26
  • $\begingroup$ The outputs and graph are identical until I check the distribution of logCPMs using a boxplot. The "boxes" of the boxplot are identical to those in the tutorial but the lower "whiskers" differ, maybe this is why my variances for the log counts don't match either. $\endgroup$ – Bio314 Jan 18 at 12:12
  • $\begingroup$ @swbarnes2 When I try to examine the hierarchical clustering using heatmap.2 function, my top 6 most variable genes ("22373" "12797" "11475" "11468" "14663" "24117") also differ from the tutorial's most variable genes ("16846" "72902" "56747" "497097" "24117" "109205"). Consequently, my heatmap is different from the tutorial's heatmap. $\endgroup$ – Bio314 Jan 21 at 15:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.