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, sys.argv, sys.argv, sys.argv, sys.argv, sys.argv, sys.argv, sys.argv] for i in range(len(inputlist)): inputfile = open(inputlist[i], 'r') for line in inputfile: splitline = line.strip().split('\t') if int(splitline) >= 5: if not splitline in genelist: genelist.append(splitline) inputfile.close() outputlist = open(sys.argv, '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 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.