Answers.Body: <p>After seeing that you want to use DESeq2, I'd recommend you do the count filtering within R, after generating the count matrix. Note that DESeq2 includes functions for loading from HTSeq-count output files:</p>
  <p><a href="https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq-count-input" rel="nofollow noreferrer">https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq-count-input</a></p>
  <p>It is not necessary to filter out low count data when processing with DESeq2, but if desired this can be done after converting to a DESeq2 object:</p>
  <blockquote>
  <p>While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.</p>
  </blockquote>
  <pre><code>keep &lt;- apply(counts(dds), 1, function(x){all(x &gt;= 5)} # for all gene counts &gt;= 5
  keep &lt;- apply(counts(dds), 1, function(x){any(x &gt;= 5)} # for any gene counts &gt;= 5
  dds &lt;- dds[keep,]
  </code></pre>
  <p><a href="https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering" rel="nofollow noreferrer">https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering</a></p>
  <p>Note that DESeq2 cares about replicates, and will almost always pin an adjusted p-value near 1 if any conditions have 1 or 2 replicates (three is the minimum number required to identify outlier values). There's not really any way to work around that other than including more replicates for the low-replicate conditions.</p>
  <p>Another program might give different results, but that won't get around the fundamental issue: two replicates is not sufficient to determine sample variation. Either accept that p-values cannot be used as a reliable estimate of model fit this case, or add more replicates.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-10-10T16:30:57.787Z
Answers.Text: Should I use any other tool like edgeR which I have not used previously?
Answers.UserDisplayName: Aranyak Goswami
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-10-10T20:05:35.803Z
Answers.Text: @AranyakGoswami, I've added additional clarification about p-values. The replicate numbers simply aren't sufficient to get reliable distribution estimates.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-09-28T01:34:10.047Z
Answers.LastActivityDate: 2022-10-11T08:56:47.687Z
Answers.LastEditDate: 2022-10-11T08:56:47.687Z
Answers.UserDisplayName: gringer
Body: <p>I have some raw counts from HTSeq after aligning with hg38 human reference genome. I want to do filtering in a way that the filtered count files should have the same number of lines. The reason behind this is that I want to use DESeq2 for differential gene expression which requires all the count files to be of the same length to carry out the analysis.</p>
  <p>I have written a python program to select reads greater than 5 from the second column.</p>
  <p>My count files look like this</p>
  <pre><code>ENSG00000000003 20
  ENSG00000000005 0
  ENSG00000000419 123
  ENSG00000000457 35
  ENSG00000000460 56
  ENSG00000000938 0
  ENSG00000000971 0
  ENSG00000001036 164
  ENSG00000001084 83
  ENSG00000001167 68
  ENSG00000001460 8
  ENSG00000001461 13
  ENSG00000001497 80
  ENSG00000001561 0
  ENSG00000001617 39
  ENSG00000001626 0
  ENSG00000001629 146
  ENSG00000001630 4
  </code></pre>
  <p>Python program I have written is as follows</p>
  <pre><code>import sys
  
  inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]  
  
  for i in range(len(inputlist)): 
      inputfile = open(inputlist[i], 'r')         
      for line in inputfile:              
          splitline = line.strip().split('\t')            
          if int(splitline[1]) &gt;= 5:      
              if not splitline[0] in genelist:    
                  genelist.append(splitline[0])   
      inputfile.close()
      
  outputlist = open(sys.argv[9], 'w') 
  for gene in genelist:
      outputlist.write(gene + '\n')
  outputlist.close()
      
  for i in range(len(inputlist)):
      inputfile = open(inputlist[i], 'r') 
      outputfile = open(sys.argv[i+10], 'w')
          for line in inputfile:      
          splitline = line.strip().split('\t')        
          if splitline[0] in genelist:                
              outputfile.write(line)  
      inputfile.close()
      outputfile.close()
  </code></pre>
  <p><code>python filtercounts.py 'mir99-1-L1count.txt' 'mir99-1-L2count.txt' 'mir99-2-L1count.txt' 'mir99-2-L2count.txt' 'mir99-3-L1count.txt' 'mir99-3-L2count.txt' 'genelist.txt' 'mir99-1-L1filtered.txt' 'mir99-1-L2filtered.txt' 'mir99-2-L1filtered.txt'  'mir99-2-L2filtered.txt'  'mir99-3-L1filtered.txt' 'mir99-3-L2filtered.txt'</code></p>
  <p>This program has worked well for me in past , but now I see that although we see many such reads removed but there are many zero reads or reads less than 5 in the dataset.</p>
  <p>Some of my output looks like this</p>
  <pre><code>ENSG00000000003 0
  ENSG00000000419 0
  ENSG00000000457 15
  ENSG00000000460 12
  ENSG00000001036 0
  ENSG00000001084 0
  ENSG00000001167 0
  ENSG00000001460 1
  ENSG00000001461 0
  ENSG00000001497 1
  ENSG00000001617 0
  ENSG00000001629 0
  ENSG00000001631 0
  ENSG00000002016 2
  ENSG00000002330 2
  ENSG00000002549 0
  ENSG00000002822 0
  ENSG00000002834 0
  ENSG00000002919 0
  ENSG00000003056 4
  ENSG00000003147 0
  ENSG00000003249 1
  </code></pre>
  <p>Any help regarding the problematic behavior of the problem will be appreciated.</p>
  
Tags: <python><deseq2><genomics>
Title: Filtering reads greater than 5 from HT-seq count files
ThreadText: # Filtering reads greater than 5 from HT-seq count files
  
  (Asked by: username_0 on 2022-09-27)
  
  I have some raw counts from HTSeq after aligning with hg38 human reference genome. I want to do filtering in a way that the filtered count files should have the same number of lines. The reason behind this is that I want to use DESeq2 for differential gene expression which requires all the count files to be of the same length to carry out the analysis.
  
  I have written a python program to select reads greater than 5 from the second column.
  
  My count files look like this
  
  ```
  ENSG00000000003 20
  ENSG00000000005 0
  ENSG00000000419 123
  ENSG00000000457 35
  ENSG00000000460 56
  ENSG00000000938 0
  ENSG00000000971 0
  ENSG00000001036 164
  ENSG00000001084 83
  ENSG00000001167 68
  ENSG00000001460 8
  ENSG00000001461 13
  ENSG00000001497 80
  ENSG00000001561 0
  ENSG00000001617 39
  ENSG00000001626 0
  ENSG00000001629 146
  ENSG00000001630 4
  
  ```
  
  Python program I have written is as follows
  
  ```
  import sys
  
  inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]  
  
  for i in range(len(inputlist)): 
      inputfile = open(inputlist[i], 'r')         
      for line in inputfile:              
          splitline = line.strip().split('\t')            
          if int(splitline[1]) >= 5:      
              if not splitline[0] in genelist:    
                  genelist.append(splitline[0])   
      inputfile.close()
      
  outputlist = open(sys.argv[9], 'w') 
  for gene in genelist:
      outputlist.write(gene + '\n')
  outputlist.close()
      
  for i in range(len(inputlist)):
      inputfile = open(inputlist[i], 'r') 
      outputfile = open(sys.argv[i+10], 'w')
          for line in inputfile:      
          splitline = line.strip().split('\t')        
          if splitline[0] in genelist:                
              outputfile.write(line)  
      inputfile.close()
      outputfile.close()
  
  ```
  
  `python filtercounts.py 'mir99-1-L1count.txt' 'mir99-1-L2count.txt' 'mir99-2-L1count.txt' 'mir99-2-L2count.txt' 'mir99-3-L1count.txt' 'mir99-3-L2count.txt' 'genelist.txt' 'mir99-1-L1filtered.txt' 'mir99-1-L2filtered.txt' 'mir99-2-L1filtered.txt' 'mir99-2-L2filtered.txt' 'mir99-3-L1filtered.txt' 'mir99-3-L2filtered.txt'`
  
  This program has worked well for me in past , but now I see that although we see many such reads removed but there are many zero reads or reads less than 5 in the dataset.
  
  Some of my output looks like this
  
  ```
  ENSG00000000003 0
  ENSG00000000419 0
  ENSG00000000457 15
  ENSG00000000460 12
  ENSG00000001036 0
  ENSG00000001084 0
  ENSG00000001167 0
  ENSG00000001460 1
  ENSG00000001461 0
  ENSG00000001497 1
  ENSG00000001617 0
  ENSG00000001629 0
  ENSG00000001631 0
  ENSG00000002016 2
  ENSG00000002330 2
  ENSG00000002549 0
  ENSG00000002822 0
  ENSG00000002834 0
  ENSG00000002919 0
  ENSG00000003056 4
  ENSG00000003147 0
  ENSG00000003249 1
  
  ```
  
  Any help regarding the problematic behavior of the problem will be appreciated.
  
  **Tags:** python, deseq2, genomics
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2022-09-28. Score: 1
  
  After seeing that you want to use DESeq2, I'd recommend you do the count filtering within R, after generating the count matrix. Note that DESeq2 includes functions for loading from HTSeq-count output files:
  
  <https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq-count-input>
  
  It is not necessary to filter out low count data when processing with DESeq2, but if desired this can be done after converting to a DESeq2 object:
  
  > While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.
  
  ```
  keep <- apply(counts(dds), 1, function(x){all(x >= 5)} # for all gene counts >= 5
  keep <- apply(counts(dds), 1, function(x){any(x >= 5)} # for any gene counts >= 5
  dds <- dds[keep,]
  
  ```
  
  <https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering>
  
  Note that DESeq2 cares about replicates, and will almost always pin an adjusted p-value near 1 if any conditions have 1 or 2 replicates (three is the minimum number required to identify outlier values). There's not really any way to work around that other than including more replicates for the low-replicate conditions.
  
  Another program might give different results, but that won't get around the fundamental issue: two replicates is not sufficient to determine sample variation. Either accept that p-values cannot be used as a reliable estimate of model fit this case, or add more replicates.
  
  > **Comments:**
  > **username_0**: Should I use any other tool like edgeR which I have not used previously?  
  > **username_1**: @username_0, I've added additional clarification about p-values. The replicate numbers simply aren't sufficient to get reliable distribution estimates.  
  
  ---
  
  