Filtering reads greater than 5 from HT-seq count files

Question

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.

Answer 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.