Bioinformatic Practices: How to conduct Enrichment Analysis with TCGA-BRCA RNA-seq (FPKM)?

Tissue:

Processed TCGA-BRCA RNAseq data in FPKM (2 excel files) have been given (normal and primary tumour; brca-rsem-fpkm-tcga.xlsx and brca-rsem-fpkm-tcga-t.xlsx, respectively). How to perform analysis to find out the the most enriched and downregulated genes and gene sets in primary breast tumours?

The two most popular methods in statistical approaches and software for RNA-seq data analysis are DESeq2 and edgeR. Both are built on the core assumption that sequencing data follow a negative binomial distribution, and they apply their own logic to normalize the raw read counts. If samples lack replicates, the negative binomial distribution may be noisy, in which case one might consider DEGseq based on the Poisson distribution. However, an underlying assumption that must be accepted for software like DESeq2 and edgeR to function properly is that there are no absolute differences between two groups of samples, or in other words, the expression of the vast majority of genes is the same across both groups.

A recent study[1] has indicated that once the analysis scales to large cohorts in population-level RNA-seq studies, with sample sizes ranging from several dozen to hundreds or thousands, the assumption of no differential gene expression across the majority of genes no longer holds. In population-level studies that have a larger sample size (at least several dozen), parameter assumptions often become unnecessary. Furthermore, with a larger number of samples, outliers are more likely to occur, violating assumptions and leading to unstable p-values, which in turn can invalidate the false discovery rate (FDR).

To explore how sample size affects the performance of six analytical methods, researchers conducted subsampling on each dataset, resulting in datasets with sample sizes ranging from 2 to 100. From the results shown above, it is evident that at a 1% FDR threshold, as long as the sample size per condition exceeds eight, the Wilcoxon rank-sum test performs comparably or better than the three parametric methods (DESeq2, edgeR, and limma-voom) and other non-parametric methods.

Therefore, the Wilcoxon rank-sum test should be employed for differential expression analysis in this case.

Unlike DESeq2, edgeR, and limma, as a non-regression-based method, the Wilcoxon rank-sum test cannot adjust for various potential confounding factors, such as differences in sequencing depth. Therefore, to utilize the Wilcoxon rank-sum test, RNA-seq samples should first be normalized to eliminate batch effects before analysis. Given that the provided type of scRNAseq is in FPKM format, which stands for fragments per kilobase of exon per million fragments mapped, FPKM cannot be used directly for differential expression analysis and should first be converted into another form. Here, we choose to transform it into TPM (Transcripts Per Million) because we will use the Wilcoxon Test to conduct a differential expression analysis.

Therefore, I have constructed a function that converts the expression matrix from FPKM to TPM.

Thus, we will obtain a table listing the results of the differential expression analysis.

2 Enrichment Analysis and Visualization

Next, let's conduct enrichment analysis, including GO, KEGG, and GSEA analysis.

Load in the differential expression data.

First, we need to specify the species databases for enrichment analysis: the species database for GO analysis is 'org.Hs.eg.db', and for KEGG analysis, it is 'hsa'.

The GO terms are grouped into different biological processes, including gland development, small GTPase mediated signal transduction, embryonic organ development, mononuclear cell differentiation, and others. The presence of these terms suggests these biological processes are potentially involved in the differences between normal and tumor samples in breast cancer. For instance, GO terms such as "gland development" and "embryonic organ development" with a high number of counts and very low p-values (dark red bars) suggest that these processes are highly relevant in the context of the tumor's biology and are significantly different compared to normal samples.

Similarly, terms like "DNA-binding transcription factor activity" and "RNA polymerase II-specific DNA-binding transcription factor binding" have a high count, which implies a significant change in gene expression related to transcription regulation between normal and tumor tissues. The blue bar at the bottom represents "GTPase regulator activity," which has a substantial number of genes involved but with a less significant p-value (as indicated by the lighter shade of red).

Due to this chart, the PI3K-Akt signaling pathway has the highest count of differentially expressed genes, suggesting that this pathway plays a significant role in the molecular events of BRCA tumor development or progression.

Pathways such as MAPK signaling, Cytokine-cytokine receptor interaction, and Calcium signaling also have high counts, underscoring their potential involvement in breast cancer pathophysiology. The Breast cancer pathway itself is highlighted, which is expected since the dataset specifically concerns breast cancer.

The color coding shows that pathways related to cancer, immune response, and cell signaling (such as Hepatocellular carcinoma, Human T-cell leukemia virus 1 infection, and Chemokine signaling pathway) are highly significant, as indicated by the deeper red colors.

Interestingly, pathways like Neuroactive ligand-receptor interaction and the Neurotrophin signaling pathway are also represented, which could suggest a role for these neuro-related pathways in breast cancer, a connection that has been explored in recent research.

The presence of non-cancer-specific pathways (such as Salmonella infection or AMPK signaling pathway) may reflect broader cellular processes affected by cancer or might be incidental findings due to the wide-ranging impact of the disease on cellular metabolism.

Processes like "small GTPase mediated signal transduction" and "mononuclear cell differentiation" have a relatively high Gene Ratio and a significant number of genes involved (larger red dots), indicating these are likely important processes in the tumor samples.

"Gland development" and "embryonic organ development" also show significant enrichment, which is expected in the context of breast cancer, as these processes are related to tissue and organ formation, and their dysregulation might play a role in tumorigenesis.

Some terms like "DNA-binding transcription factor activity" and "actin binding" appear less significant (indicated by blue dots) and involve fewer genes (smaller dots), which may suggest a lesser role in the specific aspects of breast cancer pathology being studied.

Overall, the dot plot visualizes which biological processes are most significantly associated with the genes that are differentially expressed between normal and breast cancer tumor samples.

Due to this dotplot, the Neuroactive ligand-receptor interaction and PI3K-Akt signaling pathway appear to be the most significantly enriched pathways, as indicated by their position towards the top of the plot and their red color.

MAPK signaling and Cytokine-cytokine receptor interaction pathways also show significant enrichment, suggesting these signaling pathways are potentially involved in the pathophysiology of breast cancer.

The presence of pathways such as Human papillomavirus infection may reflect underlying processes related to viral mechanisms that could be implicated in cancer or might be associated with the wider effects of tumor biology on cellular function.

The Calcium signaling pathway and cAMP signaling pathway are less significantly enriched but still noteworthy, as indicated by their purple and blue colors.

The size of the dots, particularly for pathways like Neuroactive ligand-receptor interaction and PI3K-Akt signaling, suggests a substantial number of genes are involved, highlighting the importance of these pathways in the context of BRCA tumor biology.

Now, let's discover the correlation between GO pathways.

The network plot visualizes the relationships and overlap between different GO biological process categories. Each node (circle) represents a GO term, and the size of the node reflects the number of genes associated with that term—larger nodes indicate more genes are involved. The color of the nodes corresponds to the adjusted p-value (p.adjust) of the enrichment, with the gradient from blue to red indicating less to more significant terms, respectively.

Central nodes, like "small GTPase mediated signal transduction" and "cell junction assembly," suggest they are pivotal biological processes with a larger number of genes involved, making them potentially crucial in the context of the pathology of breast cancer. Highly connected nodes indicate GO terms that share a significant number of genes, suggesting related or overlapping biological functions.

Red nodes, such as "response to metal ion" and "synapse organization," are among the most significantly enriched processes, which might be particularly relevant to the pathophysiology of breast cancer.

Peripheral nodes with fewer connections and smaller size may represent more specific processes that are altered between normal and tumor samples but are not as broadly represented across the differentially expressed genes.

The plot provides a comprehensive overview of the biological processes that are most impacted in breast cancer tumor tissue compared to normal breast tissue.

Furthermore, we can make a directed acyclic graph (DAG) representing a hierarchy of Biological Process (BP) Gene Ontology (GO) terms derived from a GO enrichment analysis.

And let's discover the correlation between KEGG pathways.

Central and larger nodes like those representing the Cell cycle and MAPK signaling pathway are significant and involve a high number of genes, suggesting these pathways play a central role in the biology of breast cancer presented in the data. Pathways such as Estrogen signaling, PI3K-Akt signaling, and Insulin signaling are also highlighted. These signaling pathways are known to be crucial in breast cancer development and progression, affecting cell growth, metabolism, and survival.

The presence of infection-related pathways like Human papillomavirus infection, Salmonella infection, and notably Coronavirus disease – COVID-19 in the context of breast cancer analysis might reflect the impact of these infections on cellular processes or could be due to the interplay between infectious agents and cancer pathways.

Peripheral nodes such as Peroxisome, Endocytosis, and Lysosome might indicate ancillary pathways that are less central but still impacted by the differential gene expression between normal and tumor samples.

The network layout with densely interconnected nodes in the center and less connected nodes on the periphery helps identify the core biological processes and secondary pathways that may contribute to the phenotype observed.

The central node colored in red, "GO:0042127 - regulation of cell proliferation", is a critical biological process often implicated in cancer, indicating that regulation of cell proliferation is significantly affected in BRCA tumor samples.

The nodes in orange may represent processes that are important but have a slightly lower significance level compared to the red node. They include processes such as "GO:0051270 - regulation of cell motion" and "GO:0030334 - regulation of cell migration", which are also relevant to cancer progression as they involve the movement and spread of cells.

Yellow nodes, such as "GO:0007155 - cell adhesion", could be indicative of processes that are important in the context of the pathology, like the adhesion of cancer cells within the tumor microenvironment or to distant tissues during metastasis.

The gradient from blue to red indicates increasing significance, with red being the most statistically significant (lowest q-value). Positive NES (to the right) indicates gene sets that are overrepresented at the top of the ranked list, suggesting a potential upregulation or activation in the condition being studied. And negative NES (to the left) indicates gene sets that are overrepresented at the bottom of the ranked list, suggesting potential downregulation or deactivation.

From this plot, we can infer which biological processes or pathways might be relevant in the condition or phenotype being studied. For example, gene sets related to "Antigen processing and presentation," "Nucleocytoplasmic transport", and "Motor proteins" show a positive NES, suggesting they might be more active. On the other hand, gene sets such as "ABC transporters" and "Insulin resistance" have a negative NES, indicating they might be less active or downregulated in the condition under study.

Pathways on the far right with positive NES values and red bars, such as Antigen processing and presentation and Nucleocytoplasmic transport, are significantly enriched in the tumor samples. This suggests these pathways are upregulated in the primary tumor samples compared to the normal samples and could be critical in the cancer’s biology.

The PI3K-Akt signaling pathway, often implicated in cancer, shows a moderate NES and a highly significant q-value. This pathway's role in cell survival, growth, and metabolism is well-documented in oncology, and its enrichment here supports its potential involvement in breast cancer tumorigenesis.

Pathways with negative NES values, depicted by blue bars, such as ABC transporters and Insulin resistance, are potentially downregulated in the tumor samples. These pathways might be less active in tumor cells or could be associated with the tumor microenvironment’s metabolic changes. The q-value color coding indicates that all the depicted pathways are significantly enriched, but those with the reddest bars, such as Nucleocytoplasmic transport and Antigen processing and presentation, are the most statistically significant.

References:

[1] Li, Yumei, Xinzhou Ge, Fanglue Peng, Wei Li, and Jingyi Jessica Li. “Exaggerated False Positives by Popular Differential Expression Methods When Analyzing Human Population Samples.” Genome Biology 23, no. 1 (March 15, 2022): 79. https://doi.org/10.1186/s13059-022-02648-4.