I am running DESeq2 comparing two cell types and exploring differential expression in a feature set with a bit more than 3 million features.

I am getting padj=1 or NA as the adjusted p-value.

Yet if I remove rows with the total sum less than 30 for example, the answer changes a lot, and I get many more adjusted p-values.

Does this have to do with the sheer number of features? Any suggestions? I thought that the independent filtering would allow me to bypass removing rows manually and that it would find the threshold with most power, but maybe so many features are problematic for this.

Any ideas?

The independent filtering is designed to automatically adjust but you should probably do some manual filtering up-front. There may be some aspects of the filtering steps that don't have high enough resolution for 3 million features (I think the default filtering lays down a grid to assess power -- code lives in genefilter package), that is they aren't finding the optimal filter here.