# What are the ways to process a list of differentially expressed genes?

We are studying six different human macrophage/dendritic cell types isolated from healthy skin. They all differ from each other in a few cell surface markers.

We are interested in the characteristics of each cell type ("marker" genes or biological processes), the differences between them (especially cell surface proteins), and their relations (e.g. "which are the closest relatives?") as can be inferred from the transcriptome, from an immunological point of view. The wider context is HIV-infection, thus HIV infection-related differences (and similarities) are of particular interest.

One naive approach is to contrast two of the most similar cell types (as inferred from previous knowledge, the sorting strategy and the PCA plot), and produce a list of differentially expressed (DE) genes. I have now a list of ~500 DE genes between two cell types, together with fold changes, expression means etc., produced with DESeq2 on bulk RNA-seq data.

What are the best practices to process this list to answer some of the above?

I've found PPI analysis on NetworkAnalyst, and GO and pathway analysis on InnateDB useful.

What other options are available?

I'm particularly interested in approaches using R.

• The question you should be asking isn't "what can we do with this?", but rather "what biological question do we want to answer?" You sound like you're aimlessly fishing for anything, which is likely not worth your time. Jun 8, 2017 at 12:11
• What do you mean by "process"? Process the list in order to do what? What question are you asking? Please edit your post and clarify what you are trying to do. As it stands, a valid answer would be that you can use wc to count the number of lines in the file. That is a type of "processing" but, obviously, not very helpful for you. The more specific you make your question, the higher the chances of your getting a useful answer. Jun 8, 2017 at 12:11
• Hi Peter, thanks for your question and welcome to Bioinformatics Stack Exchange. Bioinformatics is a very large subject area, frequently with many different solutions to the same problem. It would be helpful if you could give more story / context surrounding your problem to guide answers. For example, what organism are you working with? Are these conditions treatment vs control, or separate treatments (without a control)? What kind of results are you interested in getting out of this analysis?
– gringer
Jun 8, 2017 at 12:13
• Thank you all for the feedback. I understand that each dataset has to be studied according to the questions of interest, just wanted to get a few ideas. Also tried to make the question general so that it will be useful for more people, but I now added the specifics. Jun 8, 2017 at 12:53
• I think the question is still too broad. What knowledge do you expect to learn? Could you narrow more to one of your interests (Process, pathways, or relations between cells lines)? Which is your "immunological point of view" ?
– llrs
Jun 8, 2017 at 14:33

You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were unsurprisingly also very different.

In the first study, I did a de-novo transcriptome assembly of Schmidtea mediterranea (Smed) from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a β-catenin knockout mutant population compared to wildtype. The other researchers had some very impressive tools available at their disposal, allowing them to subsequently visualise the positional expression of differentially-expressed genes in the knockout and wildtype populations after amputation. They found a few candidate genes that had graded expression along the anterior-posterior axis in the wildtype population, and this expression disappeared in the knockout mutants. One of these genes was a teashirt (tsh) family member, and a good candidate for further discovery based on the researchers' existing knowledge of the gene pathways associated with cellular regeneration. To cut a very long story short, they experimentally validated the teashirt link through further RNAi studies, generating a double-headed Smed as a result. This was followed up in Zebrafish (Danio rerio), because there happened to be another building about 100m away that specialised in Zebrafish development.

The second study was/is a bit more similar to your current situation. We looked at sorted dendritic cell expression in mice after the application of two Th2-inducing treatment conditions (injected Nippostrongylus brasiliensis; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with DESeq2, and developed a Shiny App for result exploration after I got tired of getting asked to do a whole bunch of mechanical follow-up tasks (mostly generating graphs, lists, and heatmaps for different gene subsets). We ended up with thousands of differentially expressed genes, so it wasn't feasible to look individually at every gene for research evidence and pathways -- one scientist spent a couple of weeks looking at one interesting-looking candidate gene. The differentially-expressed genes were fed into the oh-so-expensive IPA (which gives results, but not necessarily good ones) and DAVID, among a few others. These helped a little bit, and led to what was the central idea of our paper (i.e. that there were at least two different Th2 responses), but we still had the problem of too many "significant" pathways, and lots of differentially-expressed genes that were members of a whole bunch of other pathways. We tried WGCNA; we did a bit of GSEA; and we compared Limma and DESeq2, but still couldn't find anything concrete or obvious to demonstrate what was triggering the Th2 response. There were experimental studies as well: blocking and mutant studies to investigate the two types of Th2 response (and experimentally confirm the differential expression). Our experiments and exploration are still ongoing, but we felt it was necessary to get one paper out at the intermediate stage, so just ended up taking the top 25 differentially-expressed genes in each experiment (by fold change), and finding some way of grouping some of them together (via expert knowledge and literature searches). It kind-of-sort-of worked, but we're left with a bit of frustration in that we (and the tools we used) were basically stumped by having too many biologically-relevant results.

• We should move to a chat,... but study 2 looks like a study I did: how did you make your enrichment analysis (outside IPA), specially of GO? How did you try WGCNA? Did you do a GSEA of the genes up-regulated in the comparison of a treatment vs control in the comparison of treatment1 vs treatment2?
– llrs
Jun 8, 2017 at 15:32
• WCGNA: We had about 80 different biological replicates comprising about 25 different treatment conditions, so there was enough expression variation to be able to pull out a few shared gene sets
– gringer
Jun 8, 2017 at 19:53
• GO enrichment analysis wasn't done; I don't know of any GOEA tools that use the hierarchical structure of the ontology map, and without that it's very prone to selection bias. For GSEA, we were mostly looking at treatment vs control (for various treatments and cell populations), and I think using the gene sets in the database that seemed to be immune-related. I'm not completely sure on that; it was other scientists that did the GSEA.
– gringer
Jun 8, 2017 at 19:58
• For the WGCNA, I first did a VST transform of the expression data in DESeq, normalised by the length of the longest transcript for that gene, then determined the 5,000 genes that were most variable (by MAD) across all our samples. This matrix of "lowFiltered" genes was fed into WGCNA by another scientist, and we spent a fair amount of time trying to understand the process and get the results looking right.
– gringer
Jun 8, 2017 at 20:02
• topGO package of Bioconductor uses the DAG structure of GO to find the enriched terms. WGCNA requires homogeneous samples, so if you mixed the 25 treatment conditions in those 80 samples to build a network, it would not detect the network for each treatment. You would need to do a consensus analysis (with multiset.* functions) to retrieve the network shared between each treatment, but WGCNA requires above 15 samples per condition to effectively work.
– llrs
Jun 9, 2017 at 6:56