Integrate bulk RNA and ATAC-seq genes

Question

I have up and down-regulated genes from bulk RNA result from DeSeq2. Genes that are differential in chromatin accessibility between two conditions from DiffBind. I would like to find the genes overlap between bulk RNA and ATAC-seq assay as these papers but don't know the code to get it. Would you please have a suggestion? Thank you so much 😃.

paper 1 paper 2

---

in silico DEG details I used intersect function in R which can do the same job. Which I am not sure is the experiment try to find transcription factors that cause the difference in phenotype between 2 conditions. DEG in bulk-RNA from DeSeq2 (log2 fold change >2.5 or log2 fold change < -2.5) intersect with genes differently in chromatin accessibility from DiffBind (log2 fold change >2.5 or log2 fold change < -2.5). I narrowed down to 142 genes.

For atac-seq result, after applying filter log2fold > 2.5 or log2fold < -2.5 and p value < 0.00001, I still got around 12k genes that are different between 2 conditions. Meanwhile bulk-RNA after apply filter has only around 300 genes. The number of genes in human and pig are similar so I am not sure about the number of gene I have in atac-seq

What do you think about this approach because I don't use MEME or TOMTOM?

---

I share with you what the author did.

" In order to screen new transcription factors by combining with star transcription factors. First, we used MEME software to predict the motif of Duroc pigs' peaks. Second, we compared the transcription factor database with TOMTOM to predict all transcription factors that may bind to the chromatin open regions . Third, we found all possible binding sites of MYOD, MYOG, MYF5, MYF6 and MEF2 families and did a union. A total of 99 sites were found. "Fourthly, we counted other transcription factors that may bind to these 99 sites and sorted them by the number of binding sites. All of the above statistics were manually counted."

Will read more to give you a good reply. Thank you so much 😃

Answer 1

You requested I make a response. Okay, I'm reluctant to do this firstly, all my eukaryotes genomics is out of date except for immunology stuff and secondly, this is a hardcore eukaryote gene expression question. I am really not sure about the pig-human side of this experiment.

What you have is an extremely good data set. What you appear to be doing is assessing the complexity of chromatin based gene expression, using very hardcore data. What we don't know is the variation e.g. in experiments leading to the RNA-seq. I would guess the RNA-seq is varying and being compared against atac-seq.

I have not read the papers, read the abstracts you posted. MEME and TOMTOM don't seem to be impressive in comparison to an absolute biological experiment RNA-seq vs. atac-seq. Predictions are never perfect (the very, very best predictions are 80%) and you have an absolute data set - why even be concerned about predictions based analysis?

I do know DEG and the approach you are using is very conservative. It's good to visualise this on a volcano plot. Personally I just use a probabilistic threshold and then assess the log2 for each gene group. You are using both -/+2.5 log2 and <1E-6 . Personally I would be concerned with the stringency and also look at results with e.g. just <1E-6 or <1E-6 and -/+ 0.5log2. Keep in mind I don't do eukaryotes.

The threshold criteria will be the number 1 issue that affects your results and that is the most important information I can give. So do your approach and also try a more relaxed threshold. Personally I wouldn't relax the p<1E-6, but I would relax the log2. The reasons for this are a bit complicated, but they are statistically sound.

Your intersect approach seems fine, viz.

intersect(object1, object2) 

@haci's method is cool, the above does the job too. The question is the genes vs thresholds and the experimental set ups and then the comparison between experimental setups, i.e. which genes are switching under chromatin control.

What I am not clear about is the experimental design and how you are performing DSeq2,

1. I assume you are using an internal control for RNAseq and a comparable one for atac-seq. I assume you are not comparing atac-seq verse RNAseq in the DEG.
2. The alternative is atac-seq vs. RNAseq DEG and the comparison is pig versus human. I honestly dunno.

Anway the caveat I would add is an intersect SHOULD NOT ignore whether it is +log2 or -log2 ... you need to do four intersects

• +ve verse +ve
• -v verse -ve
• +ve verse -ve
• -ve verse +ve

Again its four separate intersects.

Good luck looks a really interesting experiment and analytics look simple but powerful in context.

Answer 2

If I am getting your question right, all you need is https://bioinfogp.cnb.csic.es/tools/venny/. I assume that your two lists corresponding to RNA-seq and ATAC-seq use the same set of gene identifiers (gene symbol, Entrez ID, Ensembl ID etc).

The Results section is populated with entries that belongs to the corresponding "region" upon clicking on the regions of the Venn diagram.

Personally I have never used it but my wet-lab colleagues use it all the time.