Differential Gene Expression Analysis Advanced?
Ishani Ray
2025-10-19
Introduction
Research Question
This analysis investigates the differential gene expression patterns in the SRP033432 dataset, examining how gene expression changes across different cellular conditions. Understanding these transcriptomic signatures is crucial for understanding cellular differentiation and development.
Specific aims: - Identify genes significantly differentially expressed between conditions - Characterize the biological pathways driving cellular changes - Understand the molecular signatures underlying these cellular states
Methods
Data Processing and Quality Control
library(recount3)
library(DESeq2)
library(tidyverse)
library(ggplot2)
library(pheatmap)
library(org.Hs.eg.db)
library(clusterProfiler)
library(EnhancedVolcano)
library(DT)
library(kableExtra)
library(ggrepel)
library(scales)
library(dplyr)
library(knitr)
library(SummarizedExperiment)
library(RColorBrewer)
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(biomaRt)# Load the RNA-seq data from recount3
human_projects <- available_projects()
srp_project <- subset(human_projects, project == "SRP033432")
# Create RangedSummarizedExperiment object
rse_gene <- create_rse(srp_project)
cat("RSE Object Summary:\n")## RSE Object Summary:cat("Dimensions (genes x samples):", dim(rse_gene), "\n\n")## Dimensions (genes x samples): 63856 9# Check sample information
cat("Available metadata columns:\n")## Available metadata columns:print(colnames(colData(rse_gene)))## [1] "rail_id"
## [2] "external_id"
## [3] "study"
## [4] "sra.sample_acc.x"
## [5] "sra.experiment_acc"
## [6] "sra.submission_acc"
## [7] "sra.submission_center"
## [8] "sra.submission_lab"
## [9] "sra.study_title"
## [10] "sra.study_abstract"
## [11] "sra.study_description"
## [12] "sra.experiment_title"
## [13] "sra.design_description"
## [14] "sra.sample_description"
## [15] "sra.library_name"
## [16] "sra.library_strategy"
## [17] "sra.library_source"
## [18] "sra.library_selection"
## [19] "sra.library_layout"
## [20] "sra.paired_nominal_length"
## [21] "sra.paired_nominal_stdev"
## [22] "sra.library_construction_protocol"
## [23] "sra.platform_model"
## [24] "sra.sample_attributes"
## [25] "sra.experiment_attributes"
## [26] "sra.spot_length"
## [27] "sra.sample_name"
## [28] "sra.sample_title"
## [29] "sra.sample_bases"
## [30] "sra.sample_spots"
## [31] "sra.run_published"
## [32] "sra.size"
## [33] "sra.run_total_bases"
## [34] "sra.run_total_spots"
## [35] "sra.num_reads"
## [36] "sra.num_spots"
## [37] "sra.read_info"
## [38] "sra.run_alias"
## [39] "sra.run_center_name"
## [40] "sra.run_broker_name"
## [41] "sra.run_center"
## [42] "recount_project.project"
## [43] "recount_project.organism"
## [44] "recount_project.file_source"
## [45] "recount_project.metadata_source"
## [46] "recount_project.date_processed"
## [47] "recount_qc.aligned_reads%.chrm"
## [48] "recount_qc.aligned_reads%.chrx"
## [49] "recount_qc.aligned_reads%.chry"
## [50] "recount_qc.bc_auc.all_reads_all_bases"
## [51] "recount_qc.bc_auc.all_reads_annotated_bases"
## [52] "recount_qc.bc_auc.unique_reads_all_bases"
## [53] "recount_qc.bc_auc.unique_reads_annotated_bases"
## [54] "recount_qc.bc_auc.all_%"
## [55] "recount_qc.bc_auc.unique_%"
## [56] "recount_qc.bc_frag.count"
## [57] "recount_qc.bc_frag.kallisto_count"
## [58] "recount_qc.bc_frag.kallisto_mean_length"
## [59] "recount_qc.bc_frag.mean_length"
## [60] "recount_qc.bc_frag.mode_length"
## [61] "recount_qc.bc_frag.mode_length_count"
## [62] "recount_qc.exon_fc.all_%"
## [63] "recount_qc.exon_fc.unique_%"
## [64] "recount_qc.exon_fc_count_all.total"
## [65] "recount_qc.exon_fc_count_all.assigned"
## [66] "recount_qc.exon_fc_count_unique.total"
## [67] "recount_qc.exon_fc_count_unique.assigned"
## [68] "recount_qc.gene_fc.all_%"
## [69] "recount_qc.gene_fc.unique_%"
## [70] "recount_qc.gene_fc_count_all.total"
## [71] "recount_qc.gene_fc_count_all.assigned"
## [72] "recount_qc.gene_fc_count_unique.total"
## [73] "recount_qc.gene_fc_count_unique.assigned"
## [74] "recount_qc.intron_sum"
## [75] "recount_qc.intron_sum_%"
## [76] "recount_qc.star.%_of_chimeric_reads"
## [77] "recount_qc.star.%_of_chimeric_reads2"
## [78] "recount_qc.star.%_of_reads_mapped_to_multiple_loci"
## [79] "recount_qc.star.%_of_reads_mapped_to_multiple_loci2"
## [80] "recount_qc.star.%_of_reads_mapped_to_too_many_loci"
## [81] "recount_qc.star.%_of_reads_mapped_to_too_many_loci2"
## [82] "recount_qc.star.%_of_reads_unmapped:_other"
## [83] "recount_qc.star.%_of_reads_unmapped:_other2"
## [84] "recount_qc.star.%_of_reads_unmapped:_too_many_mismatches"
## [85] "recount_qc.star.%_of_reads_unmapped:_too_many_mismatches2"
## [86] "recount_qc.star.%_of_reads_unmapped:_too_short"
## [87] "recount_qc.star.%_of_reads_unmapped:_too_short2"
## [88] "recount_qc.star.all_mapped_reads"
## [89] "recount_qc.star.all_mapped_reads2"
## [90] "recount_qc.star.average_input_read_length"
## [91] "recount_qc.star.average_input_read_length2"
## [92] "recount_qc.star.average_mapped_length"
## [93] "recount_qc.star.average_mapped_length2"
## [94] "recount_qc.star.deletion_average_length"
## [95] "recount_qc.star.deletion_average_length2"
## [96] "recount_qc.star.deletion_rate_per_base"
## [97] "recount_qc.star.deletion_rate_per_base2"
## [98] "recount_qc.star.insertion_average_length"
## [99] "recount_qc.star.insertion_average_length2"
## [100] "recount_qc.star.insertion_rate_per_base"
## [101] "recount_qc.star.insertion_rate_per_base2"
## [102] "recount_qc.star.mapping_speed,_million_of_reads_per_hour"
## [103] "recount_qc.star.mapping_speed,_million_of_reads_per_hour2"
## [104] "recount_qc.star.mismatch_rate_per_base,_%"
## [105] "recount_qc.star.mismatch_rate_per_base,_%2"
## [106] "recount_qc.star.number_of_chimeric_reads"
## [107] "recount_qc.star.number_of_chimeric_reads2"
## [108] "recount_qc.star.number_of_input_reads"
## [109] "recount_qc.star.number_of_input_reads2"
## [110] "recount_qc.star.number_of_reads_mapped_to_multiple_loci"
## [111] "recount_qc.star.number_of_reads_mapped_to_multiple_loci2"
## [112] "recount_qc.star.number_of_reads_mapped_to_too_many_loci"
## [113] "recount_qc.star.number_of_reads_mapped_to_too_many_loci2"
## [114] "recount_qc.star.number_of_reads_unmapped:_other"
## [115] "recount_qc.star.number_of_reads_unmapped:_other2"
## [116] "recount_qc.star.number_of_reads_unmapped:_too_many_mismatches"
## [117] "recount_qc.star.number_of_reads_unmapped:_too_many_mismatches2"
## [118] "recount_qc.star.number_of_reads_unmapped:_too_short"
## [119] "recount_qc.star.number_of_reads_unmapped:_too_short2"
## [120] "recount_qc.star.number_of_splices:_at/ac"
## [121] "recount_qc.star.number_of_splices:_at/ac2"
## [122] "recount_qc.star.number_of_splices:_annotated_(sjdb)"
## [123] "recount_qc.star.number_of_splices:_annotated_(sjdb)2"
## [124] "recount_qc.star.number_of_splices:_gc/ag"
## [125] "recount_qc.star.number_of_splices:_gc/ag2"
## [126] "recount_qc.star.number_of_splices:_gt/ag"
## [127] "recount_qc.star.number_of_splices:_gt/ag2"
## [128] "recount_qc.star.number_of_splices:_non-canonical"
## [129] "recount_qc.star.number_of_splices:_non-canonical2"
## [130] "recount_qc.star.number_of_splices:_total"
## [131] "recount_qc.star.number_of_splices:_total2"
## [132] "recount_qc.star.uniquely_mapped_reads_%"
## [133] "recount_qc.star.uniquely_mapped_reads_%2"
## [134] "recount_qc.star.uniquely_mapped_reads_number"
## [135] "recount_qc.star.uniquely_mapped_reads_number2"
## [136] "recount_qc.junction_count"
## [137] "recount_qc.junction_coverage"
## [138] "recount_qc.junction_avg_coverage"
## [139] "recount_qc.star.number_of_input_reads_both"
## [140] "recount_qc.star.all_mapped_reads_both"
## [141] "recount_qc.star.number_of_chimeric_reads_both"
## [142] "recount_qc.star.number_of_reads_mapped_to_multiple_loci_both"
## [143] "recount_qc.star.number_of_reads_mapped_to_too_many_loci_both"
## [144] "recount_qc.star.number_of_reads_unmapped:_other_both"
## [145] "recount_qc.star.number_of_reads_unmapped:_too_many_mismatches_both"
## [146] "recount_qc.star.number_of_reads_unmapped:_too_short_both"
## [147] "recount_qc.star.uniquely_mapped_reads_number_both"
## [148] "recount_qc.star.%_mapped_reads_both"
## [149] "recount_qc.star.%_chimeric_reads_both"
## [150] "recount_qc.star.%_reads_mapped_to_multiple_loci_both"
## [151] "recount_qc.star.%_reads_mapped_to_too_many_loci_both"
## [152] "recount_qc.star.%_reads_unmapped:_other_both"
## [153] "recount_qc.star.%_reads_unmapped:_too_many_mismatches_both"
## [154] "recount_qc.star.%_reads_unmapped:_too_short_both"
## [155] "recount_qc.star.uniquely_mapped_reads_%_both"
## [156] "recount_seq_qc.min_len"
## [157] "recount_seq_qc.max_len"
## [158] "recount_seq_qc.avg_len"
## [159] "recount_seq_qc.#distinct_quality_values"
## [160] "recount_seq_qc.#bases"
## [161] "recount_seq_qc.%a"
## [162] "recount_seq_qc.%c"
## [163] "recount_seq_qc.%g"
## [164] "recount_seq_qc.%t"
## [165] "recount_seq_qc.%n"
## [166] "recount_seq_qc.avgq"
## [167] "recount_seq_qc.errq"
## [168] "recount_pred.sample_acc.y"
## [169] "recount_pred.curated.type"
## [170] "recount_pred.curated.tissue"
## [171] "recount_pred.pattern.predict.type"
## [172] "recount_pred.pred.type"
## [173] "recount_pred.curated.cell_type"
## [174] "recount_pred.curated.cell_line"
## [175] "BigWigURL"# Check sample information
coldata_raw <- colData(rse_gene) %>%
as.data.frame() %>%
dplyr::select(starts_with("sra")) %>%
head() %>%
kable()
coldata_raw| sra.sample_acc.x | sra.experiment_acc | sra.submission_acc | sra.submission_center | sra.submission_lab | sra.study_title | sra.study_abstract | sra.study_description | sra.experiment_title | sra.design_description | sra.sample_description | sra.library_name | sra.library_strategy | sra.library_source | sra.library_selection | sra.library_layout | sra.paired_nominal_length | sra.paired_nominal_stdev | sra.library_construction_protocol | sra.platform_model | sra.sample_attributes | sra.experiment_attributes | sra.spot_length | sra.sample_name | sra.sample_title | sra.sample_bases | sra.sample_spots | sra.run_published | sra.size | sra.run_total_bases | sra.run_total_spots | sra.num_reads | sra.num_spots | sra.read_info | sra.run_alias | sra.run_center_name | sra.run_broker_name | sra.run_center | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| SRR1041710 | SRS510053 | SRX386044 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277968: Sample1_hpSC_Undiff; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;Human Parthenogenic Stem Cells|source_name;;Human Parthenogenic Stem Cells | GEO Accession;;GSM1277968 | NA | GSM1277968 | Sample1_hpSC_Undiff | 4521789200 | 22608946 | 2015-12-24 13:33:03 | 3282304956 | 4521789200 | 22608946 | 2 | 22608946 | index:0,count:22608946,average:100,stdev:0|index:1,count:22608946,average:100,stdev:0 | GSM1277968_r1 | GEO | NA | NA |
| SRR1041711 | SRS510054 | SRX386045 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277969: Sample2_hpSC_Undiff; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;Human Parthenogenic Stem Cells|source_name;;Human Parthenogenic Stem Cells | GEO Accession;;GSM1277969 | NA | GSM1277969 | Sample2_hpSC_Undiff | 3825829200 | 19129146 | 2015-12-24 13:33:03 | 2774244721 | 3825829200 | 19129146 | 2 | 19129146 | index:0,count:19129146,average:100,stdev:0|index:1,count:19129146,average:100,stdev:0 | GSM1277969_r1 | GEO | NA | NA |
| SRR1041712 | SRS510055 | SRX386046 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277970: Sample3_hpSC_Undiff; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;Human Parthenogenic Stem Cells|source_name;;Human Parthenogenic Stem Cells | GEO Accession;;GSM1277970 | NA | GSM1277970 | Sample3_hpSC_Undiff | 4377262800 | 21886314 | 2015-12-24 13:33:03 | 3180460261 | 4377262800 | 21886314 | 2 | 21886314 | index:0,count:21886314,average:100,stdev:0|index:1,count:21886314,average:100,stdev:0 | GSM1277970_r1 | GEO | NA | NA |
| SRR1041713 | SRS510056 | SRX386047 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277971: Sample4_hpSC_NSC; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;hpSC derived Neural Stem Cells|source_name;;Neural Stem Cells | GEO Accession;;GSM1277971 | NA | GSM1277971 | Sample4_hpSC_NSC | 10811876400 | 54059382 | 2015-12-24 13:33:03 | 7811347465 | 10811876400 | 54059382 | 2 | 54059382 | index:0,count:54059382,average:100,stdev:0|index:1,count:54059382,average:100,stdev:0 | GSM1277971_r1 | GEO | NA | NA |
| SRR1041714 | SRS510057 | SRX386048 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277972: Sample5_hpSC_NSC; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;hpSC derived Neural Stem Cells|source_name;;Neural Stem Cells | GEO Accession;;GSM1277972 | NA | GSM1277972 | Sample5_hpSC_NSC | 7978232200 | 39891161 | 2015-12-24 13:33:03 | 5768106176 | 7978232200 | 39891161 | 2 | 39891161 | index:0,count:39891161,average:100,stdev:0|index:1,count:39891161,average:100,stdev:0 | GSM1277972_r1 | GEO | NA | NA |
| SRR1041715 | SRS510059 | SRX386049 | SRA115150 | GEO | NA | Differential expression of human parthenogenic stem cells, neural stem cells and DA progenitors. | Neural stem cells (NSC) derived from human parthenogenic stem cells (hpSC) have been observed to show stronger positive functional effects than hpSC-derived dopaminergic neuron precursors (DAP) in treatment of induced Parkinson Disease in animal models. RNAseq of the two types of cells were normalized and analyzed to compare gene expression profiles. Overall design: cDNA library of hpsC, NSC and DAP triplicates were sequenced using Illumina HiSeq 2000. The sequence reads were mapped to hg19 reference genome and hits that passed quality filters were analyzed for differential expression. | NA | GSM1277973: Sample6_hpSC_NSC; Homo sapiens; RNA-Seq | NA | NA | NA | RNA-Seq | TRANSCRIPTOMIC | cDNA | paired | NA | NA | RNA was isolated from cells (RNeasy Plus Mini kit, Qiagen), quantified (Qubit RNA Assay Kit, Life Technologies, Inc.), and quality controlled (RNA6000 Nano Kit and BioAnalyzer 2100, Agilent). Approximately 600 ng was used as input for the Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) and sequencing libraries were created according to the manufacturer’s protocol. Briefly, poly-A containing mRNA molecules are purified using a reagent containing poly-T oligo-attached magnetic beads. After purification, the mRNA is fragmented and first strand cDNA is produced using random primers and reverse transcriptase. Second strand cDNA synthesis is then prepared using DNA Polymerase I and RNase H. The cDNA is then ligated to indexed Illumina adapters and enriched with PCR to create the cDNA library. | Illumina HiSeq 2000 | cell line;;LLC12PH|cell type;;hpSC derived Neural Stem Cells|source_name;;Neural Stem Cells | GEO Accession;;GSM1277973 | NA | GSM1277973 | Sample6_hpSC_NSC | 6085372000 | 30426860 | 2015-12-24 13:33:03 | 4407595990 | 6085372000 | 30426860 | 2 | 30426860 | index:0,count:30426860,average:100,stdev:0|index:1,count:30426860,average:100,stdev:0 | GSM1277973_r1 | GEO | NA | NA |
# Extract phenotype information - adjust based on actual column names
coldata <- colData(rse_gene) %>%
as.data.frame() %>%
mutate(
sample_id = rownames(.),
# Adjust this based on your actual sample names/metadata
cell_type = case_when(
str_detect(sra.sample_title, "Undiff") ~ "Undiff",
str_detect(sra.sample_title, "NSC") ~ "NSC",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(cell_type))
# Update RSE with filtered samples
rse_filtered <- rse_gene[, coldata$sample_id]
colData(rse_filtered) <- DataFrame(coldata)
cat("Sample counts per cell type:\n")## Sample counts per cell type:print(table(coldata$cell_type))##
## NSC Undiff
## 3 3# Clean Ensembl IDs (remove version numbers)
ensembl_ids <- rownames(rse_filtered)
ensembl_ids_clean <- gsub("\\..*", "", ensembl_ids)
# Step 1: Quick annotation with org.Hs.eg.db
gene_symbols_orgdb <- mapIds(org.Hs.eg.db,
keys = ensembl_ids_clean,
column = "SYMBOL",
keytype = "ENSEMBL",
multiVals = "first")
orgdb_annotated <- sum(!is.na(gene_symbols_orgdb))
cat(" ✓", orgdb_annotated, "genes annotated with org.Hs.eg.db\n\n")## ✓ 38709 genes annotated with org.Hs.eg.db# Step 2: Use biomaRt to fill gaps and get additional info
ensembl <- useEnsembl(biomart = "genes",
dataset = "hsapiens_gene_ensembl",
mirror = "useast")
# Query in chunks to avoid timeout
chunk_size <- 5000
n_chunks <- ceiling(length(ensembl_ids_clean) / chunk_size)
biomart_results <- list()
for (i in 1:n_chunks) {
start_idx <- (i - 1) * chunk_size + 1
end_idx <- min(i * chunk_size, length(ensembl_ids_clean))
chunk_ids <- ensembl_ids_clean[start_idx:end_idx]
cat(" Chunk", i, "/", n_chunks, "...")
tryCatch({
chunk_result <- getBM(
attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype'),
filters = 'ensembl_gene_id',
values = chunk_ids,
mart = ensembl
)
biomart_results[[i]] <- chunk_result
cat(" ✓\n")
Sys.sleep(0.5)
}, error = function(e) {
cat(" ERROR\n")
biomart_results[[i]] <- NULL
})
}## Chunk 1 / 13 ... ✓
## Chunk 2 / 13 ... ✓
## Chunk 3 / 13 ... ✓
## Chunk 4 / 13 ... ✓
## Chunk 5 / 13 ... ✓
## Chunk 6 / 13 ... ✓
## Chunk 7 / 13 ... ✓
## Chunk 8 / 13 ... ✓
## Chunk 9 / 13 ... ✓
## Chunk 10 / 13 ... ✓
## Chunk 11 / 13 ... ✓
## Chunk 12 / 13 ... ✓
## Chunk 13 / 13 ... ✓biomart_annotation <- bind_rows(biomart_results)
cat("\n ✓", nrow(biomart_annotation), "genes retrieved from biomaRt\n\n")##
## ✓ 62016 genes retrieved from biomaRt# Merge annotations
gene_annotation <- tibble(
ensembl_id = ensembl_ids,
ensembl_id_clean = ensembl_ids_clean,
gene_symbol_orgdb = gene_symbols_orgdb
) %>%
left_join(biomart_annotation, by = c("ensembl_id_clean" = "ensembl_gene_id")) %>%
mutate(
gene_symbol = coalesce(gene_symbol_orgdb, external_gene_name),
display_name = coalesce(gene_symbol, ensembl_id_clean)
) %>%
dplyr::select(ensembl_id, ensembl_id_clean, gene_symbol, gene_biotype, display_name)
# Summary
total_annotated <- sum(!is.na(gene_annotation$gene_symbol))
pct_annotated <- round(100 * total_annotated / nrow(gene_annotation), 1)
improvement <- total_annotated - orgdb_annotated
cat("=== ANNOTATION SUMMARY ===\n")## === ANNOTATION SUMMARY ===cat(" Total genes:", nrow(gene_annotation), "\n")## Total genes: 63856cat(" Annotated:", total_annotated, "(", pct_annotated, "%)\n")## Annotated: 62215 ( 97.4 %)cat(" Improvement:", improvement, "additional genes from biomaRt\n")## Improvement: 23506 additional genes from biomaRtcat(" Unannotated:", nrow(gene_annotation) - total_annotated, "\n\n")## Unannotated: 1641# Create DESeq2 object
dds <- DESeqDataSet(
rse_filtered,
design = ~ cell_type
)
# Determine smallest group size
smallest_group_size <- min(table(coldata$cell_type))
cat("\nSmallest group size:", smallest_group_size, "\n")##
## Smallest group size: 3# Pre-filtering: remove genes with very low counts
# At least 10 counts in smallest group
keep <- rowSums(counts(dds) >= 10) >= smallest_group_size
dds <- dds[keep, ]
cat("Genes retained after filtering:", sum(keep), "\n")## Genes retained after filtering: 32935cat("Genes removed:", nrow(rse_filtered) - sum(keep), "\n\n")## Genes removed: 30921# Set reference level
dds$cell_type <- relevel(dds$cell_type, ref = "Undiff")
cat("Reference condition set to:", levels(dds$cell_type)[1], "\n")## Reference condition set to: UndiffDifferential Expression Analysis
# Run DESeq2
dds <- DESeq(dds)
conditions <- levels(dds$cell_type)
conditions## [1] "Undiff" "NSC"# Get results
res_nsc <- results(dds, contrast = c("cell_type", "NSC", "Undiff"))
# Create summary
cat("DESeq2 Results Summary - NSC vs Undiff:\n")## DESeq2 Results Summary - NSC vs Undiff:print(summary(res_nsc))##
## out of 32935 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 10816, 33%
## LFC < 0 (down) : 12028, 37%
## outliers [1] : 305, 0.93%
## low counts [2] : 0, 0%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## NULLResults
PCA Plot
# VST transformation for visualization
vst_data <- vst(dds, blind = FALSE)
# Extract PCA data
pca_data <- plotPCA(vst_data, intgroup = "cell_type", returnData = TRUE)
pca_vars <- attr(pca_data, "percentVar")
# Create plot
pca_plot <- ggplot(pca_data, aes(x = PC1, y = PC2, color = cell_type, shape = cell_type)) +
geom_point(size = 4, alpha = 0.8) +
xlab(paste0("PC1 (", round(pca_vars[1], 2), "%)")) +
ylab(paste0("PC2 (", round(pca_vars[2], 2), "%)")) +
theme_minimal() +
theme(
legend.position = "bottom",
text = element_text(size = 12, face = "bold"),
plot.title = element_text(hjust = 0.5)
) +
labs(title = "PCA: Sample Clustering by Cell Type", color = "Cell Type", shape = "Cell Type") +
scale_color_brewer(palette = "Set1")
print(pca_plot)
Differential Expression Results
# Function to create DEG table with gene annotations
create_deg_table <- function(res_obj, gene_annotation) {
res_df <- as.data.frame(res_obj) %>%
rownames_to_column("ensembl_id") %>%
left_join(gene_annotation, by = c("ensembl_id" = "ensembl_id")) %>%
filter(!is.na(padj)) %>%
arrange(padj) %>%
mutate(
significant = ifelse(padj < 0.05 & abs(log2FoldChange) > 1, "Yes", "No"),
direction = ifelse(log2FoldChange > 0, "Up-regulated", "Down-regulated")
) %>%
dplyr::select(
ensembl_id, ensembl_id_clean, gene_symbol, display_name, log2FoldChange, stat, pvalue, padj, significant, direction, gene_biotype
)
return(res_df)
}
# Get DEG tables
degs_nsc <- create_deg_table(res_nsc, gene_annotation)
# Count significant genes
sig_count <- sum(degs_nsc$significant == "Yes")
up_count <- sum(degs_nsc$significant == "Yes" & degs_nsc$direction == "Up-regulated")
down_count <- sum(degs_nsc$significant == "Yes" & degs_nsc$direction == "Down-regulated")
cat("Significantly differentially expressed genes (padj < 0.05, |log2FC| > 1):\n")## Significantly differentially expressed genes (padj < 0.05, |log2FC| > 1):cat("Total:", sig_count, "\n")## Total: 13348cat("Up-regulated:", up_count, "\n")## Up-regulated: 6025cat("Down-regulated:", down_count, "\n\n")## Down-regulated: 7323Interactive DEG Table: NSC vs hpSC
degs_nsc %>%
filter(significant == "Yes") %>%
dplyr::select(gene_symbol, display_name, log2FoldChange, padj, direction, gene_biotype) %>%
arrange(padj) %>%
datatable(
extensions = 'Buttons',
options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel'),
pageLength = 15,
scrollX = TRUE,
columnDefs = list(
list(
targets = which(names(degs_nsc) %in% c("padj", "log2FoldChange")) - 1,
render = JS("function(data, type, row) {return parseFloat(data).toExponential(2);}")
)
)
),
caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: left;',
paste("Significantly Differentially Expressed Genes:")
)
)Volcano Plot
# Prepare data for volcano plot
volcano_data <- degs_nsc %>%
mutate(
log10_padj = -log10(padj),
label = ifelse(significant == "Yes" & abs(log2FoldChange) > 1.5, gene_symbol, NA)
) %>%
filter(!is.na(log2FoldChange))
volcano_plot <- ggplot(volcano_data, aes(x = log2FoldChange, y = log10_padj)) +
geom_point(aes(color = significant), alpha = 0.6, size = 2) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray50") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray50") +
ggrepel::geom_text_repel(aes(label = label), max.overlaps = 10, size = 3) +
scale_color_manual(values = c("No" = "gray70", "Yes" = "#e74c3c")) +
xlab("log2(Fold Change)") +
ylab("-log10(Adjusted p-value)") +
theme_minimal() +
theme(
legend.position = "bottom",
text = element_text(size = 12, face = "bold"),
plot.title = element_text(hjust = 0.5)
) +
labs(title = "Volcano Plot: NSC vs hpSC", color = "Significant")
print(volcano_plot)
Heatmap of Top Differentially Expressed Genes
# Get top 10 upregulated and downregulated genes
top_up <- degs_nsc %>%
filter(direction == "Up-regulated", significant == "Yes") %>%
arrange(padj) %>%
pull(ensembl_id) %>%
head(10)
top_down <- degs_nsc %>%
filter(direction == "Down-regulated", significant == "Yes") %>%
arrange(padj) %>%
pull(ensembl_id) %>%
head(10)
top_genes <- c(top_up, top_down)
# Create gene labels
gene_labels <- degs_nsc %>%
filter(ensembl_id %in% top_genes) %>%
arrange(match(ensembl_id, top_genes)) %>%
pull(display_name)
# Extract normalized counts for plotting
vst_mat <- assay(vst_data)[top_genes, ]
rownames(vst_mat) <- gene_labels
annotation <- coldata %>% dplyr::select(cell_type)
rownames(annotation) <- coldata$sample_id
# Create heatmap
pheatmap(
vst_mat,
annotation_col = annotation,
annotation_colors = list(
cell_type = c(
Undiff = "#e74c3c",
NSC = "#3498db"
)
),
main = "Top 10 Up and Down-regulated Genes (NSC vs Undiff)",
fontsize_row = 9,
fontsize_col = 10,
color = colorRampPalette(c("blue", "white", "red"))(50),
scale = "row"
)
Gene Set Enrichment Analysis
# Prepare ranked gene list for GSEA
ranked_genes <- degs_nsc %>%
filter(!is.na(stat)) %>%
arrange(desc(stat)) %>%
pull(stat, name = gene_symbol)
# Remove NA gene symbols
ranked_genes <- ranked_genes[!is.na(names(ranked_genes))]
# Convert gene symbols to Entrez IDs
gene_entrez <- mapIds(
org.Hs.eg.db,
keys = names(ranked_genes),
column = "ENTREZID",
keytype = "SYMBOL"
)
# Remove duplicates and NAs
ranked_genes_entrez <- ranked_genes[!is.na(gene_entrez)]
names(ranked_genes_entrez) <- gene_entrez[!is.na(gene_entrez)]
ranked_genes_entrez <- ranked_genes_entrez[!duplicated(names(ranked_genes_entrez))]
cat("Genes mapped to Entrez IDs:", length(ranked_genes_entrez), "\n")## Genes mapped to Entrez IDs: 25936# Load pathways
gsea_results <- gseKEGG(
geneList = ranked_genes_entrez,
organism = "hsa",
pvalueCutoff = 0.1,
verbose = FALSE
)
colnames(as.data.frame(gsea_results))## [1] "ID" "Description" "setSize" "enrichmentScore"
## [5] "NES" "pvalue" "p.adjust" "qvalue"
## [9] "rank" "leading_edge" "core_enrichment"# Get top pathways
top_pathways <- gsea_results@result %>%
arrange(p.adjust) %>% # use qvalues instead of padj
head(5)
datatable(
top_pathways %>%
dplyr::select(ID, Description, NES, p.adjust, core_enrichment) %>%
mutate(
NES = round(NES, 3),
padj = format(p.adjust, scientific = TRUE, digits = 2)
),
caption = "Top 5 KEGG Pathways by Adjusted p-value"
)# Plot top 5 pathways
gsea_plot <- ggplot(
top_pathways,
aes(x = reorder(Description, NES), y = NES, fill = p.adjust)
) +
geom_col() +
coord_flip() +
scale_fill_gradient(low = "#e74c3c", high = "#3498db") +
xlab("KEGG Pathway") +
ylab("Normalized Enrichment Score (NES)") +
theme_minimal() +
theme(
text = element_text(size = 11, face = "bold"),
plot.title = element_text(hjust = 0.5),
axis.text.y = element_text(size = 10)
) +
labs(
title = "Top 5 Enriched KEGG Pathways",
fill = "Adj. p-value"
)
print(gsea_plot)
Discussion
Key Findings
Based on the differential expression analysis comparing NSC and Undiff, we identified 13348 significantly dysregulated genes (adjusted p-value < 0.05, |log2 fold-change| > 1). This includes 6025 up-regulated and 7323 down-regulated genes. The PCA plot demonstrates clear separation between conditions, indicating substantial transcriptomic remodeling.
Biological Interpretation
The genes identified as significantly up-regulated likely include key transcription factors, signaling molecules, and structural proteins essential for the differentiated state. Conversely, down-regulated genes probably reflect the silencing of undifferentiated cell-specific programs.
The GSEA results highlight enriched biological pathways, providing mechanistic insights into the molecular processes driving the observed transcriptomic changes.
Limitations and Considerations
- Gene annotation relies on current Ensembl and gene symbol databases
- Multiple testing correction (adjusted p-values) is applied to control false discovery rate
- Results are limited to protein-coding genes in the current analysis
- Biological validation through qRT-PCR or other methods is recommended
Future Directions
Validation Experiments: Perform qRT-PCR or Western blotting on top candidate genes.
Pathway Analysis: Conduct more detailed functional annotation using GO terms or other pathway databases.
Single-Cell Analysis: Investigate cell-type-specific subpopulations using scRNA-seq if available.
Temporal Dynamics: Examine time-course data to understand dynamic gene expression patterns.
Comparative Studies: Compare findings with related published datasets for conservation of key signatures.
Session Information
sessionInfo()## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] biomaRt_2.58.2 msigdbr_7.5.1
## [3] enrichplot_1.24.4 RColorBrewer_1.1-3
## [5] knitr_1.48 scales_1.4.0
## [7] kableExtra_1.4.0 DT_0.33
## [9] EnhancedVolcano_1.22.0 ggrepel_0.9.6
## [11] clusterProfiler_4.12.6 org.Hs.eg.db_3.19.1
## [13] AnnotationDbi_1.66.0 pheatmap_1.0.12
## [15] lubridate_1.9.3 forcats_1.0.0
## [17] stringr_1.5.1 dplyr_1.1.4
## [19] purrr_1.0.2 readr_2.1.5
## [21] tidyr_1.3.1 tibble_3.2.1
## [23] ggplot2_4.0.0 tidyverse_2.0.0
## [25] DESeq2_1.44.0 recount3_1.14.0
## [27] SummarizedExperiment_1.34.0 Biobase_2.64.0
## [29] GenomicRanges_1.56.1 GenomeInfoDb_1.40.1
## [31] IRanges_2.38.1 S4Vectors_0.42.1
## [33] BiocGenerics_0.50.0 MatrixGenerics_1.16.0
## [35] matrixStats_1.4.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.4.1 BiocIO_1.14.0 bitops_1.0-9
## [4] ggplotify_0.1.2 filelock_1.0.3 R.oo_1.26.0
## [7] polyclip_1.10-7 XML_3.99-0.17 lifecycle_1.0.4
## [10] httr2_1.0.5 lattice_0.22-6 MASS_7.3-61
## [13] crosstalk_1.2.1 magrittr_2.0.3 sass_0.4.9
## [16] rmarkdown_2.28 jquerylib_0.1.4 yaml_2.3.10
## [19] sessioninfo_1.2.2 cowplot_1.1.3 DBI_1.2.3
## [22] abind_1.4-8 zlibbioc_1.50.0 R.utils_2.12.3
## [25] ggraph_2.2.1 RCurl_1.98-1.16 yulab.utils_0.1.7
## [28] tweenr_2.0.3 rappdirs_0.3.3 GenomeInfoDbData_1.2.12
## [31] tidytree_0.4.6 svglite_2.1.3 codetools_0.2-20
## [34] DelayedArray_0.30.1 xml2_1.3.6 DOSE_3.30.5
## [37] ggforce_0.4.2 tidyselect_1.2.1 aplot_0.2.3
## [40] UCSC.utils_1.0.0 farver_2.1.2 viridis_0.6.5
## [43] BiocFileCache_2.12.0 GenomicAlignments_1.40.0 jsonlite_1.8.9
## [46] tidygraph_1.3.1 systemfonts_1.3.1 progress_1.2.3
## [49] tools_4.4.1 treeio_1.28.0 snow_0.4-4
## [52] Rcpp_1.1.0 glue_1.7.0 gridExtra_2.3
## [55] SparseArray_1.4.8 xfun_0.46 qvalue_2.36.0
## [58] withr_3.0.2 fastmap_1.2.0 digest_0.6.36
## [61] timechange_0.3.0 R6_2.5.1 gridGraphics_0.5-1
## [64] colorspace_2.1-1 GO.db_3.19.1 dichromat_2.0-0.1
## [67] RSQLite_2.3.7 R.methodsS3_1.8.2 generics_0.1.3
## [70] data.table_1.16.0 rtracklayer_1.64.0 prettyunits_1.2.0
## [73] graphlayouts_1.2.0 httr_1.4.7 htmlwidgets_1.6.4
## [76] S4Arrays_1.4.1 scatterpie_0.2.4 pkgconfig_2.0.3
## [79] gtable_0.3.6 blob_1.2.4 S7_0.2.0
## [82] XVector_0.44.0 shadowtext_0.1.4 htmltools_0.5.8.1
## [85] fgsea_1.30.0 png_0.1-8 ggfun_0.1.6
## [88] rstudioapi_0.16.0 tzdb_0.4.0 reshape2_1.4.4
## [91] rjson_0.2.23 nlme_3.1-166 curl_5.2.3
## [94] cachem_1.1.0 parallel_4.4.1 restfulr_0.0.15
## [97] pillar_1.10.1 grid_4.4.1 vctrs_0.6.5
## [100] dbplyr_2.5.0 evaluate_1.0.0 cli_3.6.3
## [103] locfit_1.5-9.10 compiler_4.4.1 Rsamtools_2.20.0
## [106] rlang_1.1.4 crayon_1.5.3 labeling_0.4.3
## [109] plyr_1.8.9 fs_1.6.4 stringi_1.8.4
## [112] viridisLite_0.4.2 BiocParallel_1.38.0 babelgene_22.9
## [115] Biostrings_2.72.1 lazyeval_0.2.2 GOSemSim_2.30.2
## [118] Matrix_1.7-0 hms_1.1.3 patchwork_1.3.1
## [121] bit64_4.5.2 KEGGREST_1.44.1 highr_0.11
## [124] igraph_2.0.3 memoise_2.0.1 bslib_0.8.0
## [127] ggtree_3.12.0 fastmatch_1.1-4 bit_4.5.0
## [130] ape_5.8-1 gson_0.1.0