Answers.Body: <p>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.</p>
  
  <p>In the <a href="http://www.sciencedirect.com/science/article/pii/S221112471401047X" rel="noreferrer">first study</a>, I did a <em>de-novo</em> transcriptome assembly of <em>Schmidtea mediterranea (Smed)</em> from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a <em>β-catenin</em> 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 <em>teashirt (tsh)</em> 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 <em>teashirt</em> link through further RNAi studies, generating a double-headed <em>Smed</em> as a result. This was followed up in Zebrafish <em>(Danio rerio)</em>, because there happened to be another building about 100m away that specialised in Zebrafish development.</p>
  
  <p>The <a href="https://www.growkudos.com/publications/10.1084%252Fjem.20160470" rel="noreferrer">second study</a> 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 <em>Nippostrongylus brasiliensis</em>; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with <a href="http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html" rel="noreferrer">DESeq2</a>, and developed a <a href="http://shiny.rstudio.com/gallery/" rel="noreferrer">Shiny App</a> 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 <a href="https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/" rel="noreferrer">IPA</a> (which gives results, but not necessarily good ones) and <a href="https://david.ncifcrf.gov/" rel="noreferrer">DAVID</a>, 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 <a href="https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/Rpackages/WGCNA/" rel="noreferrer">WGCNA</a>; we did a bit of <a href="http://software.broadinstitute.org/gsea/index.jsp" rel="noreferrer">GSEA</a>; and we compared <a href="https://bioconductor.org/packages/release/bioc/html/limma.html" rel="noreferrer">Limma</a> 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.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T15:32:38.703Z
Answers.Text: 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?
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T19:53:53.583Z
Answers.Text: 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
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T19:58:40.423Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T20:02:34.047Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T06:56:47.940Z
Answers.Text: 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.
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T07:49:57.360Z
Answers.Text: Thanks, I might give topGO a try once I've sorted out the proportionality calculations
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T12:38:57.067Z
Answers.LastActivityDate: 2017-06-08T14:37:56.240Z
Answers.LastEditDate: 2017-06-08T14:37:56.240Z
Answers.UserDisplayName: gringer
Body: <p>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.</p>
  
  <p>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.</p>
  
  <p>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.</p>
  
  <p>What are the best practices to process this list to answer some of the above?</p>
  
  <p>I've found PPI analysis on <a href="http://networkanalyst.ca" rel="noreferrer">NetworkAnalyst</a>, and  GO and pathway analysis on <a href="http://innatedb.com" rel="noreferrer" title="InnateDB">InnateDB</a> useful.</p>
  
  <p>What other options are available?</p>
  
  <p>I'm particularly interested in approaches using <code>R</code>.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T12:11:01.780Z
Comments.Text: 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.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T12:11:23.883Z
Comments.Text: 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T12:13:21.860Z
Comments.Text: 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?
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T12:53:41.173Z
Comments.Text: 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.
Comments.UserDisplayName: Peter
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T14:33:46.177Z
Comments.Text: 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" ?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-01-11T04:16:22.693Z
Comments.Text: It's best to stick to questions that are as specific as possible, and let other people decide on their own generalisations about the answers. There are so many different things that can be done in bioinformatics even for a single task that it's no useful to speculate about what other similar things could be done.
Comments.UserDisplayName: gringer
Tags: <rna-seq><r><deseq2><networks>
Title: What are the ways to process a list of differentially expressed genes?
ThreadText: # What are the ways to process a list of differentially expressed genes?
  
  (Asked by: username_0 on 2017-06-08)
  
  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](http://networkanalyst.ca), and GO and pathway analysis on [InnateDB](http://innatedb.com "InnateDB") useful.
  
  What other options are available?
  
  I'm particularly interested in approaches using `R`.
  
  **Tags:** rna-seq, r, deseq2, networks
  
  > **Comments:**
  > **username_1**: 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.  
  > **username_2**: 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.  
  > **username_3**: Hi username_0, 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?  
  > **username_0**: 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.  
  > **username_4**: 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" ?  
  > **username_3**: It's best to stick to questions that are as specific as possible, and let other people decide on their own generalisations about the answers. There are so many different things that can be done in bioinformatics even for a single task that it's no useful to speculate about what other similar things could be done.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_3 on 2017-06-08. Score: 7 | Accepted answer
  
  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](http://www.sciencedirect.com/science/article/pii/S221112471401047X), 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](https://www.growkudos.com/publications/10.1084%252Fjem.20160470) 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](http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html), and developed a [Shiny App](http://shiny.rstudio.com/gallery/) 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](https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/) (which gives results, but not necessarily good ones) and [DAVID](https://david.ncifcrf.gov/), 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](https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/Rpackages/WGCNA/); we did a bit of [GSEA](http://software.broadinstitute.org/gsea/index.jsp); and we compared [Limma](https://bioconductor.org/packages/release/bioc/html/limma.html) 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.
  
  > **Comments:**
  > **username_4**: 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?  
  > **username_3**: 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  
  > **username_3**: 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.  
  > **username_3**: 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.  
  > **username_4**: 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.  
  > **username_3**: Thanks, I might give topGO a try once I've sorted out the proportionality calculations  
  
  ---
  
  