This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Biological replicates and depth seq

Hello, everyone I have a question about samples replicates of RNA sequencing depth.

When I run the heatmap in my set of experiments with DESeq2, the replicates are not always together and the position can change based on the number of genes that I consider (100 or 1000). I noticed that this effect is more evident in samples where the number of reads is very different although normalized with DESeq2. Is it possible to have variance between biological replicates also because of different sequencing depth?

Thank you

rna-seq deseq2
• 487 views • link
updated by Charles Warden Charles Warden • written by Morris_Chair Morris_Chair

DESeq2`s normalisation method accounts for sequencing depth, so, that should not be an issue. Can you please show the code that you used when processing your data? Some information on your wet-lab protocol would also help, and whether the replicates are biological or technical.. Edit: you mentioned that your replicates are biological, so, variation is expected.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin,

The experiment consists of treating the cells in 4 different conditions (knock-out) where each condition has a biological replicate, then in total I have 8 sample. I used two different cell lines treated the same way so at the end I have have 16 samples.

The script used is the following:

    samples <- read.table(file.path("/salmonSFquant/sampleinfo.txt"), header=TRUE)
files <- file.path("/salmonSFquant", samples$condition, "quant.sf")

condition <- factor(c(rep(“sample1”,4),rep(“sample2”,4), rep(“sample3”,4), rep(“sample”,4)))
model.matrix(~condition)

condition=relevel(condition, “sample4”)
all(file.exists(files))

txdb <- makeTxDbFromGFF(file="/gencode.v19.annotation.gtf")
saveDb(x=txdb, file = "gencode.v19.annotation.Txdb")
k <- keys(txdb, keytype = "TXNAME")

tx2gene <- select(txdb, k, "GENEID", "TXNAME")

k <- keys(txdb, keytype = "GENEID")
df <- select(txdb, keys = k,  columns = "TXNAME", keytype = "GENEID")
tx2gene <- df[, 2:1]

write.table(tx2gene, "tx2gene.gencode.v19_RPMI_PF382.csv")
txi.salmon <- tximport(files, type="salmon", tx2gene=tx2gene,ignoreAfterBar=TRUE  )

library(DESeq2)

coldata <- data.frame( condition=factor(c(rep("sample1",4),rep("sample2",4), rep("sample3",4), rep("sample",4))))
rownames(samples) <- colnames(txi.salmon$counts)
dds <- DESeqDataSetFromTximport(txi.salmon,coldata, design= ~condition)

rld <- rlog( dds )
head( assay(rld) )

dds$condition<- relevel(dds$condition, ref = “sample4”)

dds <-DESeq(dds)
res <- results(dds, alpha=0.5)

resultsNames(dds)
resLFC <- lfcShrink(dds, coef=4, res=res, type="apeglm")
resFile <-na.omit(resLFC)
resFile <-resFile[order(resFile$padj),]
write.csv(as.data.frame(resFile),file=“experiment_.csv", quote=F)


mytable <- counts(dds, normalized = TRUE)
write.csv(mytable, file=“experiment_norm_counts.csv")

After that I start with the plots, please let me know if you need also of that,

Thank you

• 0 views • link
written by Morris_Chair Morris_Chair

Kevin your silence scares me :D

• 0 views • link
written by Morris_Chair Morris_Chair

Nobody has any obligation at all to keep on answering.

• 0 views • link
written by WouterDeCoster WouterDeCoster

sure WonterDeCoster you're right, but Kevin is not like others .. (I'm kidding), sorry if violate the rules ..

• 0 views • link
written by Morris_Chair Morris_Chair

It's not a written rule. It's just etiquette.

• 0 views • link
written by WouterDeCoster WouterDeCoster

Can you show a PCA? 

plotPCA(rld, "condition", 1000)
• 0 views • link
written by ATpoint ATpoint

Hi ATpoint, Here you can see the PCA plot, I'm repeating one condition (the green on the right side) because there is too much variance. This are two different cell lines well separated

BioStar2

Thank you

• 0 views • link
written by Morris_Chair Morris_Chair

Any idea why these two "groups" exist, is this maybe a strong batch effect? This cannot be all based on seq. depth, I guess.

• 0 views • link
written by ATpoint ATpoint

The two groups exist because they are two different cell lines, means a different genetic background. What really matters to me is that the condition for each cell line have to be close which in the heatmap differs based on the number of genes selected

Thanks

• 0 views • link
written by Morris_Chair Morris_Chair

I was thinking that probably processing the two different cell lines in the same simultaneously is not a good idea...

• 0 views • link
written by Morris_Chair Morris_Chair

1 answer

I would usually expect more noise among genes with lower expression. However, I think it matters i) what is the total number of aligned reads per sample and ii) what is the variation between samples.

In other words, if your samples are low coverage, perhaps I understand a large difference in the ratio of total reads (such as 2 million versus 6 million reads). However, if you have more than 20-30 million reads, it is unusual to have 20 million reads for one sample and 60 million reads for another sample (in the same experiment). If the samples are from different experiments (which may have used different protocols, different strandedness, different sample degradation, etc.), then I think you could have additional difficulties beyond the normalization.

If all samples had similar reads (with say 20 million reads, or perhaps even 10 million aligned reads), perhaps focusing on the more highly expressed genes could help control for factors not completely correct (or maybe even subtly over-corrected) with normalization.

Plus, as Keven points out, if you have biological replicates, that can decrease concordance among replicates (but I think it would be best if you could control as many technical factors as possible, in order to convince yourself these are really biological differences).

• 0 views • link
written by Charles Warden Charles Warden

Hi Charles,

My samples are from the same experiment, How would you focus on high more expressed genes? still using DESeq2?

Thank you for you answer

• 0 views • link
written by Morris_Chair Morris_Chair

Hi Morris.

I would use an independently calculated FPKM expression value for quality control plots and assessing differential expression strategies (as long as you have replicates, I would preferably test out at least DESeq2 / limma-voom / and edgeR for every project).

Even though I would eventually like to be able to have a paper to show that you can't lock-down the methods for every project ahead of time (at least not if you are working with a wide variety of topics). However, I mention the topic in this discussion A: RNA-seq dispersion estimation (which I most recently cited here)

That said, I'm not sure what I should specifically say in this situation.

In general, if you wanted to filter candidates, you should try looking at log2(FPKM + 1) values instead of log2(FPKM + 0.1) values. However, again, I don't know how big a deal the total number of reads is (or how many genes even have have been detected at an expression level of FPKM of 1 or 0.1). So, I'm not sure if that is what I should recommend for your project.

• 0 views • link
written by Charles Warden Charles Warden

Log in to answer this question.

More posts like this