2
$\begingroup$

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

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ 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... $\endgroup$
    – llrs
    Commented Feb 8, 2019 at 15:41
  • $\begingroup$ @llrs I'm interested in co-expression in tumor, may be I will do spearately. For the second, ofcourse I go based on pvalue. $\endgroup$
    – stack_learner
    Commented Feb 8, 2019 at 16:01
  • 1
    $\begingroup$ 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 $\endgroup$
    – llrs
    Commented Feb 8, 2019 at 16:05
  • $\begingroup$ 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? $\endgroup$
    – stack_learner
    Commented Feb 8, 2019 at 16:08
  • $\begingroup$ 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) $\endgroup$
    – llrs
    Commented Feb 8, 2019 at 16:18

3 Answers 3

Reset to default
4
$\begingroup$

You have several options to approach this with 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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – stack_learner
    Commented Feb 8, 2019 at 18:12
  • 1
    $\begingroup$ 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. $\endgroup$
    – Peter Langfelder
    Commented Feb 10, 2019 at 1:02
  • 1
    $\begingroup$ 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. $\endgroup$
    – Peter Langfelder
    Commented Feb 10, 2019 at 1:02
  • $\begingroup$ 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]; $\endgroup$
    – stack_learner
    Commented Feb 10, 2019 at 16:49
  • 1
    $\begingroup$ As far as I can tell (the formatting is messed up) it should correctly give you the top hubs in the brown module. $\endgroup$
    – Peter Langfelder
    Commented Feb 13, 2019 at 23:14
2
$\begingroup$

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
Improve this answer
$\endgroup$
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.