Answers.Body: <p>Generally, the recommended way of doing such analysis is exploring data by a dimensionality reduction technique such as PCA, e.g. implemented in DESeq2 via <code>plotPCA</code>. I personally always also look at MA-plots (<code>plotMA</code>). This will inform how samples cluster, and if there is evidence that some unwanted technical variation is present, or that samples are simply not different, or that normalization might be an issue. Without these plots one really cannot help here.</p>
  <p>It is of course possible that there are no DEGs. No difference is not necessarily a technical error.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-03T14:15:22.567Z
Answers.Text: thank you very much for your answer! I added the PCA and MA plots to the question. Based on your answer, I assume that since there is no clear clustering the case here is most likely that there is too much variation within the replicates and no issue with the analysis itself.
Answers.UserDisplayName: Adrianna V
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-04T06:59:43.590Z
Answers.Text: It might be helpful if you project other variables (other than treatment-control, for example batch, sex, ...) in your PCA plot, especially the replicates that you define above. For example I see treatment-control pairs in your PCA plots, if these pairs belong to individual replicates/batches, you can add that to your design model as a covariate and expect to find more DE genes.
Answers.UserDisplayName: haci
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-04T10:52:38.787Z
Answers.Text: I can only second what swbarnes2 is saying. The treatment does not seem to be doing much. Has the experiment been done in batches?
Answers.UserDisplayName: ATpoint
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-04T13:34:45.503Z
Answers.Text: Yes the experiment has been done in batches, but the replicates treatment-control do not correspond to the seemingly nice pairing in the PCA plot.
Answers.UserDisplayName: Adrianna V
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-03T10:15:32.293Z
Answers.LastActivityDate: 2024-06-03T10:15:32.293Z
Answers.UserDisplayName: ATpoint
Answers.Body: <p>The reason is that, according to your PCA, your treatment doesn't do very much.  It's not what you want to hear, but that's probably the truth.  You probably didn't do anything wrong, it's just what your experiment is.</p>
  <p>I would not find it surprising that your knockdown affects only a handful of genes.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-06-03T17:57:02.533Z
Answers.LastActivityDate: 2024-06-03T17:57:02.533Z
Answers.UserDisplayName: swbarnes2
Body: <p>I am fairly new at RNA-sequencing analysis, but I have attempted to analyse the data I obtained from RNA sequencing on my own by following an online EdX class and by basing the majority of the analysis on the recommendations of Azenta in the following guide: A Quick Start Guide to RNA-Seq Data Analysis .</p>
  <p>For context, I am looking at making two different comparisons each time between human A549 cells treated with control siRNA (non-targeted siRNA) and a targeted knockdown using siRNA. I have 3 sequenced replicates for each condition.</p>
  <p>As mentioned I mainly followed the steps recommended by azenta:</p>
  <ol>
  <li><p>FASTQC to look at the quality of the samples (everything seems fine there)</p>
  </li>
  <li><p>Trimming reads with trimmomatic (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 MINLEN:36)</p>
  </li>
  <li><p>Aligning the reads to the hg19 reference genome with STAR (the relevant files for the genome assembly, fasta and gtf, were retrieved from the gencode database)</p>
  </li>
  <li><p>Calculate gene hit counts with featureCounts (in R; isGTFAnnotationFile=TRUE, countMultiMappingReads=FALSE, GTF.featureType=&quot;exon&quot;, largestOverlap=TRUE, GTF.attrType=&quot;gene_id&quot;)</p>
  </li>
  <li><p>Compare hit counts between groups using DESeq2 with LFCshrink type Apeglem.</p>
  </li>
  </ol>
  <p>Without taking fold change into account, by only looking at the p-adjusted value &lt; 0.05, we only identify 19 genes matching this parameter in one comparison and 2 in the other comparison. In each comparison, one of the genes that is identified as being differentially expressed is the one that was knocked down by our siRNA treatment.</p>
  <p>We were expecting to identify much more differentially expressed genes, but since we do identify the knocked down gene, we are unsure whether something is wrong with our pipeline or with the samples? Or whether this could be interpreted as a &quot;normal&quot; amount of differentially expressed genes?</p>
  <p>I tried different types of LFCshrink normalization which did not significantly alter the results. I also tried to run the easily accessible airway dataset through the DESeq2 pipeline, this dataset did seem to behave normally. I also tried a genome assembly based on files obtained from the UCSC database, which also did not alter the data.</p>
  <p>Any guidance or advise would be appreciated to help figure out whether something is wrong with our analysis or not! Thank you for your help.</p>
  <p><a href="https://i.sstatic.net/oJ285BFA.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/oJ285BFA.png" alt="MA plot of comparison with 19 genes with padj &lt; 0.05" /></a></p>
  <p><a href="https://i.sstatic.net/pBq1DZkf.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/pBq1DZkf.png" alt="PCA plot of comparison with 19 genes with padj &lt; 0.05" /></a></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-05-31T21:07:30.230Z
Comments.Text: Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-06-05T09:19:08.040Z
Comments.Text: What is the paired read length? Is it a stranded sequencing? Often there is no point of trimming RNA-Seq prior to mapping. STAR should be able to better figure out where to map a given pair of reads when you keep "a not so great quality but okeish sequence-wise" read ends.
Comments.UserDisplayName: darked89
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-06-05T10:01:49.070Z
Comments.Text: While I doubt this will change the results drastically, I think it makes sense to use the best available version of the genome (T2T) and the best gene annotation. At least hg38/Gencode
Comments.UserDisplayName: darked89
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-06-06T19:42:09.800Z
Comments.Text: Can you please change the Y range of that MA plot so that it includes all the data? There are currently what looks like hundreds of out-of-range genes.
Comments.UserDisplayName: gringer
Tags: <rna-seq><differential-expression><deseq2><star>
Title: Reasons for extremely low number of DESeq identified differentially expressed genes after RNAseq?
ThreadText: # Reasons for extremely low number of DESeq identified differentially expressed genes after RNAseq?
  
  (Asked by: username_0 on 2024-05-31)
  
  I am fairly new at RNA-sequencing analysis, but I have attempted to analyse the data I obtained from RNA sequencing on my own by following an online EdX class and by basing the majority of the analysis on the recommendations of Azenta in the following guide: A Quick Start Guide to RNA-Seq Data Analysis .
  
  For context, I am looking at making two different comparisons each time between human A549 cells treated with control siRNA (non-targeted siRNA) and a targeted knockdown using siRNA. I have 3 sequenced replicates for each condition.
  
  As mentioned I mainly followed the steps recommended by azenta:
  
  1. FASTQC to look at the quality of the samples (everything seems fine there)
  2. Trimming reads with trimmomatic (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 MINLEN:36)
  3. Aligning the reads to the hg19 reference genome with STAR (the relevant files for the genome assembly, fasta and gtf, were retrieved from the gencode database)
  4. Calculate gene hit counts with featureCounts (in R; isGTFAnnotationFile=TRUE, countMultiMappingReads=FALSE, GTF.featureType="exon", largestOverlap=TRUE, GTF.attrType="gene\_id")
  5. Compare hit counts between groups using DESeq2 with LFCshrink type Apeglem.
  
  Without taking fold change into account, by only looking at the p-adjusted value < 0.05, we only identify 19 genes matching this parameter in one comparison and 2 in the other comparison. In each comparison, one of the genes that is identified as being differentially expressed is the one that was knocked down by our siRNA treatment.
  
  We were expecting to identify much more differentially expressed genes, but since we do identify the knocked down gene, we are unsure whether something is wrong with our pipeline or with the samples? Or whether this could be interpreted as a "normal" amount of differentially expressed genes?
  
  I tried different types of LFCshrink normalization which did not significantly alter the results. I also tried to run the easily accessible airway dataset through the DESeq2 pipeline, this dataset did seem to behave normally. I also tried a genome assembly based on files obtained from the UCSC database, which also did not alter the data.
  
  Any guidance or advise would be appreciated to help figure out whether something is wrong with our analysis or not! Thank you for your help.
  
  [![MA plot of comparison with 19 genes with padj < 0.05](https://i.sstatic.net/oJ285BFA.png)](https://i.sstatic.net/oJ285BFA.png)
  
  [![PCA plot of comparison with 19 genes with padj < 0.05](https://i.sstatic.net/pBq1DZkf.png)](https://i.sstatic.net/pBq1DZkf.png)
  
  **Tags:** rna-seq, differential-expression, deseq2, star
  
  > **Comments:**
  > **username_1**: Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.  
  > **username_2**: What is the paired read length? Is it a stranded sequencing? Often there is no point of trimming RNA-Seq prior to mapping. STAR should be able to better figure out where to map a given pair of reads when you keep "a not so great quality but okeish sequence-wise" read ends.  
  > **username_2**: While I doubt this will change the results drastically, I think it makes sense to use the best available version of the genome (T2T) and the best gene annotation. At least hg38/Gencode  
  > **username_3**: Can you please change the Y range of that MA plot so that it includes all the data? There are currently what looks like hundreds of out-of-range genes.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_4 on 2024-06-03. Score: 2 | Accepted answer
  
  Generally, the recommended way of doing such analysis is exploring data by a dimensionality reduction technique such as PCA, e.g. implemented in DESeq2 via `plotPCA`. I personally always also look at MA-plots (`plotMA`). This will inform how samples cluster, and if there is evidence that some unwanted technical variation is present, or that samples are simply not different, or that normalization might be an issue. Without these plots one really cannot help here.
  
  It is of course possible that there are no DEGs. No difference is not necessarily a technical error.
  
  > **Comments:**
  > **username_0**: thank you very much for your answer! I added the PCA and MA plots to the question. Based on your answer, I assume that since there is no clear clustering the case here is most likely that there is too much variation within the replicates and no issue with the analysis itself.  
  > **username_5**: It might be helpful if you project other variables (other than treatment-control, for example batch, sex, ...) in your PCA plot, especially the replicates that you define above. For example I see treatment-control pairs in your PCA plots, if these pairs belong to individual replicates/batches, you can add that to your design model as a covariate and expect to find more DE genes.  
  > **username_4**: I can only second what username_6 is saying. The treatment does not seem to be doing much. Has the experiment been done in batches?  
  > **username_0**: Yes the experiment has been done in batches, but the replicates treatment-control do not correspond to the seemingly nice pairing in the PCA plot.  
  
  ---
  
  ### Answer by username_6 on 2024-06-03. Score: 2
  
  The reason is that, according to your PCA, your treatment doesn't do very much. It's not what you want to hear, but that's probably the truth. You probably didn't do anything wrong, it's just what your experiment is.
  
  I would not find it surprising that your knockdown affects only a handful of genes.
  
  
  ---
  
  