<p>A cursory search to <a href="https://bedtools.readthedocs.io/en/latest/content/tools/closest.html" rel="nofollow noreferrer">bedtools documentation</a> will reveal the <code>bedtools closest</code> feature - which might be exactly what you are looking in your third question. You can download the list of protein-coding genes from Ensembl biomart or UCSC table browser in bed format, and also convert the list of your lncRNAs into a bed compatible file. Rest is easy:</p>

<pre><code>bedtools closest -a lncRNAs.bed -b proteinCoding.bed [-iu/-id] -D ref
</code></pre>

CC BY-SA 4.0
2019-02-11T09:38:07.903Z
2019-02-11T09:38:07.903Z
Siddharth
<blockquote>
  <p>2) If not WGCNA, if I use cor function in R with method Pearson, on
  what cutoff I should select the target genes?</p>
</blockquote>

<p>Good call, this is a nice approach. It is a good question because the central issue is that a purists approach of 0.7 results in insufficient information to understand the behaviour of the data.</p>

<p>The guys here were using 0.4, in my experience co-variance of 0.5 and above is good evidence, whilst 0.7+ is strong evidence. </p>

<p>I've done this (not on RNAseq) and I used 0.25 to 0.45 and 0.45 and above then used a scatterplot with one half of the graph showing low co-variance and the other half strong co-variance. So in one sense I ignored my own advice (0.5+). In my defence its easy to reassign the boundaries. Each scatterplot (1000s) was then subject to tSNE analysis to assess the group behaviour. Its a different context but gives an idea.</p>

<p>Networks are cool for visualisation but they are not particularly  quantitive for most (not all) network theory.</p>

<hr>

<p>Just to explain to everyone else, what happens if you get boundary effects with the data and you are assigning a discrete boundary on a continous distribution (its a negative binomial distribution in my case). I observed a vast amount of results  occurring between 0.2 to 0.25 - but 0.2 isn't good evidence, so that data has to be lost. These sort of boundaries were in days before genomic throughput. Thus in my context I didn't consider a hard boundary particularly informative.</p>

CC BY-SA 4.0
2019-02-08T15:43:57.853Z
2019-02-08T16:01:58.460Z
2019-02-08T16:01:58.460Z
M__
<p>You have several options to approach this with <a href="https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/" rel="nofollow noreferrer">WGCNA</a> (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.</p>

<p>Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).  </p>

<p>I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples. </p>

<p>I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.</p>

CC BY-SA 4.0
2019-02-08T18:12:59.773Z
Thanks a lot for the answer Peter. With WGCNA, I detected the modules and from that modules having mRNAs and lncRNAs how to detect hub genes inside modules?
stack_learner
CC BY-SA 4.0
2019-02-10T01:02:45.817Z
Hub genes are the ones with high intramodular or eigengene based connectivity. Use the function signedKME to calculate kME values (it's really just a correlation of module eigengenes with individual expression profiles) and choose the top genes.
Peter Langfelder
CC BY-SA 4.0
2019-02-10T01:02:59.680Z
There is no hard and fast criterion for what constitutes a hub gene; I'd probably look at the first maybe 20 genes for a smaller module (less than say 200 genes) and maybe a few more (say up to 50) for larger modules; if you have a module of more than a thousand gene, looking at the top 100 may be useful. You can also go by kME (i.e., correlation); anything higher than 0.7 is pretty strong, and sometimes you can go down to say 0.6 to still find useful answers.
Peter Langfelder
CC BY-SA 4.0
2019-02-10T16:49:14.957Z
Thanks peter. After module detection, I used this code. Do you think this is right?         MEList = moduleEigengenes(datExpr, colors = dynamicColors);
MEs = MEList$eigengenes;

datKME=signedKME(datExpr, MEs);
head(datKME);

    FilterGenes= abs(datKME$kMEbrown)>.8;
table(FilterGenes);
dimnames(data.frame(datExpr))[[2]][FilterGenes];
stack_learner
CC BY-SA 4.0
2019-02-11T12:16:30.570Z
Hi Peter, Could you please check the above code whether that is right way to get hub genes?
stack_learner
CC BY-SA 4.0
2019-02-13T23:14:22.243Z
As far as I can tell (the formatting is messed up) it should correctly give you the top hubs in the brown module.
Peter Langfelder
CC BY-SA 4.0
2019-02-08T17:42:18.030Z
2019-02-10T00:56:02.353Z
2019-02-10T00:56:02.353Z
Peter Langfelder
<p>I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)</p>

<p>I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.</p>

<p>Few doubts:</p>

<p>1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?</p>

<p>2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?</p>

<p>3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.</p>

<p>Thank you</p>

CC BY-SA 4.0
2019-02-08T15:41:44.760Z
I think you should ask a separate question for your third question (remember to show what you have tried/searched...). For the question 1: which is your biological question? Are you interested in correlation between lncRNAs and mRNAs regardless if there is a tumor or not? Or are you interested in knowing if the co-expression in normal samples are different from the tumor? For question 2: You should be aware that you can find lots of false positive on those correlations. You should at least restrict yourself to those significantly correlated. But I wouldn't recommend to follow this path...
llrs
CC BY-SA 4.0
2019-02-08T16:01:35.347Z
@llrs I'm interested in co-expression in tumor, may be I will do spearately. For the second, ofcourse I go based on pvalue.
stack_learner
CC BY-SA 4.0
2019-02-08T16:05:49.977Z
In this case it is better if you only use the tumor samples. If you go based on p-value consider correcting the p-value to reduce the amount of false positives
llrs
CC BY-SA 4.0
2019-02-08T16:08:04.080Z
And for WGCNA, after module detection, there are some important genes Im interested in were kept in the grey module which is unassigned genes. What to do in this situation?
stack_learner
CC BY-SA 4.0
2019-02-08T16:18:58.147Z
You can then say that these genes are not co-expressed following a clear pattern. Either they don't follow a pattern, or are the ones regulating the expression of the other modules, or you don't have enough power/samples or it is too noisy to find a signal... The important genes might be for other reasons not for co-expression between mRNA and lncRNA (or they still are but you can't detect it)
llrs
<r><rna-seq><bioconductor><rna><wgcna>
Question about Co-expression analysis and finding targets for lncRNAs
# Question about Co-expression analysis and finding targets for lncRNAs

(Asked by: username_0 on 2019-02-08)

I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)

I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.

Few doubts:

1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?

2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?

3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.

Thank you

**Tags:** r, rna-seq, bioconductor, rna, wgcna
> **Comments:**
> **username_1**: I think you should ask a separate question for your third question (remember to show what you have tried/searched...). For the question 1: which is your biological question? Are you interested in correlation between lncRNAs and mRNAs regardless if there is a tumor or not? Or are you interested in knowing if the co-expression in normal samples are different from the tumor? For question 2: You should be aware that you can find lots of false positive on those correlations. You should at least restrict yourself to those significantly correlated. But I wouldn't recommend to follow this path...  
> **username_0**: @llrs I'm interested in co-expression in tumor, may be I will do spearately. For the second, ofcourse I go based on pvalue.  
> **username_1**: In this case it is better if you only use the tumor samples. If you go based on p-value consider correcting the p-value to reduce the amount of false positives  
> **username_0**: And for WGCNA, after module detection, there are some important genes Im interested in were kept in the grey module which is unassigned genes. What to do in this situation?  
> **username_1**: You can then say that these genes are not co-expressed following a clear pattern. Either they don't follow a pattern, or are the ones regulating the expression of the other modules, or you don't have enough power/samples or it is too noisy to find a signal... The important genes might be for other reasons not for co-expression between mRNA and lncRNA (or they still are but you can't detect it)  

---

## 3 answers

### Answer by username_4 on 2019-02-08. Score: 4 | Accepted answer

You have several options to approach this with [WGCNA](https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/) (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.

Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).

I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.

I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.

> **Comments:**
> **username_0**: Thanks a lot for the answer Peter. With WGCNA, I detected the modules and from that modules having mRNAs and lncRNAs how to detect hub genes inside modules?  
> **username_4**: Hub genes are the ones with high intramodular or eigengene based connectivity. Use the function signedKME to calculate kME values (it's really just a correlation of module eigengenes with individual expression profiles) and choose the top genes.  
> **username_4**: There is no hard and fast criterion for what constitutes a hub gene; I'd probably look at the first maybe 20 genes for a smaller module (less than say 200 genes) and maybe a few more (say up to 50) for larger modules; if you have a module of more than a thousand gene, looking at the top 100 may be useful. You can also go by kME (i.e., correlation); anything higher than 0.7 is pretty strong, and sometimes you can go down to say 0.6 to still find useful answers.  
> **username_0**: Thanks peter. After module detection, I used this code. Do you think this is right?         MEList = moduleEigengenes(datExpr, colors = dynamicColors); MEs = MEList$eigengenes;  datKME=signedKME(datExpr, MEs); head(datKME);      FilterGenes= abs(datKME$kMEbrown)>.8; table(FilterGenes); dimnames(data.frame(datExpr))[[2]][FilterGenes];  
> **username_0**: Hi Peter, Could you please check the above code whether that is right way to get hub genes?  
> **username_4**: As far as I can tell (the formatting is messed up) it should correctly give you the top hubs in the brown module.  

---

### Answer by username_2 on 2019-02-11. Score: 2

A cursory search to bedtools documentation will reveal the bedtools closest feature - which might be exactly what you are looking in your third question. You can download the list of protein-coding genes from Ensembl biomart or UCSC table browser in bed format, and also convert the list of your lncRNAs into a bed compatible file. Rest is easy:
bedtools closest -a lncRNAs.bed -b proteinCoding.bed [-iu/-id] -D ref


---

### Answer by username_3 on 2019-02-08. Score: 1

> 2) If not WGCNA, if I use cor function in R with method Pearson, on
> what cutoff I should select the target genes?

Good call, this is a nice approach. It is a good question because the central issue is that a purists approach of 0.7 results in insufficient information to understand the behaviour of the data.

The guys here were using 0.4, in my experience co-variance of 0.5 and above is good evidence, whilst 0.7+ is strong evidence.

I've done this (not on RNAseq) and I used 0.25 to 0.45 and 0.45 and above then used a scatterplot with one half of the graph showing low co-variance and the other half strong co-variance. So in one sense I ignored my own advice (0.5+). In my defence its easy to reassign the boundaries. Each scatterplot (1000s) was then subject to tSNE analysis to assess the group behaviour. Its a different context but gives an idea.

Networks are cool for visualisation but they are not particularly quantitive for most (not all) network theory.

---

Just to explain to everyone else, what happens if you get boundary effects with the data and you are assigning a discrete boundary on a continous distribution (its a negative binomial distribution in my case). I observed a vast amount of results occurring between 0.2 to 0.25 - but 0.2 isn't good evidence, so that data has to be lost. These sort of boundaries were in days before genomic throughput. Thus in my context I didn't consider a hard boundary particularly informative.


---

