# BulkRNAseq analyses

This repository contains the key scripts used to generate the final tables and figures featured in the accompanying manuscript:  
**Zonari et al**, _Expansion of Hematopoietic Stem and Progenitor Cells from Human Mobilized Peripheral Blood for Gene Therapy_, XXX, 2025.  

## Abstract 
Ex vivo expansion of mobilized peripheral blood (mPB) hematopoietic stem cells (HSCs)
represents a promising approach to advance cell and gene therapy strategies yet is hampered by
loss of stem cell function when applying commonly used culture protocols. We performed an indepth
characterization of mPB expansion cultures by single cell RNA sequencing, which
highlighted differentiation trajectories with preservation of lineage fidelity in committed
progenitors. Defining a putative HSC cluster allowed an estimation of transduction efficiency in
ex vivo cultures, which correlated with long-term gene marking in xenografts and patients enrolled
in a gene therapy study. We then developed a clinically translatable, GMP-compliant process to
expand lentivirus (LV)-transduced HSCs from mPB of pediatric patients and adult donors, by
biologically informed protocol improvements of cytokine supplementation, media choice, timing
of LV transduction and combinations of small molecules preventing the activation of
differentiation programs. Our optimized process outperforms validated state-of-the-art cord blood
expansion protocols when applied to mPB. LV integration site analysis and genomic barcodebased
clonal tracking provided definitive proof for symmetric HSC self-renewal divisions
occurring during ex vivo culture. These results warrant clinical testing of this HSC
transduction/expansion process in an upcoming clinical gene therapy trial for autosomal recessive
osteopetrosis (EU CT 2024-518972-30), addressing the need of obtaining a high number of
functioning osteoclast progenitors for rapid bone marrow niche remodeling, where gene-corrected
HSC can engraft and allow long-term disease correction.

GEO: [GSE292605](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE292605). 

Below a brief description with Bulk and scRNAseq workflow adopted in this work.

**Bulk RNAseq** analysis was performed using a standard pipeline that includes the follwing steps:
1. Quality control by [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
2. Trimming of bad quality reads with [TrimGalore](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/)<details><summary>Running command</summary>trim_galore --quality 20 --fastqc --length 25 --output_dir {outdir} --paired {input.r1} {inout.r2}</details>
3. Alignment with [STAR](https://github.com/alexdobin/STAR)
    <details><summary>Running command</summary>
            "STAR " +
            "--runThreadN {threads} " +
            "--genomeDir {input.genome} " +
            "--readFilesIn {params.trim_seq} " +
            "--outSAMstrandField intronMotif " +
            "--outFileNamePrefix {params.aln_seq_prefix} " +
            "--outSAMtype BAM SortedByCoordinate " +
            "--outSAMmultNmax 1 " +
            "--outFilterMismatchNmax 10 " +
            "--outReadsUnmapped Fastx " +
            "--readFilesCommand zcat "
    </details>
4. Gene expression quantification with [FeatureCounts](https://academic.oup.com/bioinformatics/article/30/7/923/232889)
    <details><summary>Running command</summary>
            "featureCounts " +
            "-a {input.annot} " +
            "-o {output.fcount} " +
            "-g gene_name " +
            "-p -B -C " +
            "-s {params.strand} " +
            "--minOverlap 10 " +
            "-T {threads} " +
            "{input.bams} "
    </details>

5. DGE analysis with [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html):

We employed a data analysis workflow that relied on edgeR DEGs identification setting custom biological variation coefficients (BCV) to 0.1-0.3 for the different 1vs1 comparisons performed (72h vs 24h – in each batch separately). Counts were normalized using TMM method and differential test was performed with exactTest function provided by the edgeR package (v3.32). Genes with an adjusted p.value (Benjiamini Hochberg FDR method) < 0.05 were considered differentially expressed.  

Details are provided in the script: [DGE.R](http://www.bioinfotiget.it/gitlab/custom/zonari_mpbhscexp_2025/zonari_mpbhscexp_2025_bulkrnaseq/-/blob/main/DGE/DGE.R)

6. Dowstream functional Analysis with [ClusterProfiler](https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html).  
    In order to retrieve functional annotation from DE analysis, we performed **O**ver **R**epresentation **A**nalysis and **G**ene **S**et        **E**nrichment **A**nalysis by using the functions provided by the package.  
    **ORA** analysis was performed in particular on Gene Ontology Terms (Biological Process sub-group) retrieved by the function.  
    **GSEA** analysis was performed against MSigDB gene sets C1 to C6 + Hallmarks and Senescence signatures (*.gmt files)

Details are provided in the script: [Fig5F_processing_and_plot.R](http://www.bioinfotiget.it/gitlab/custom/zonari_mpbhscexp_2025/zonari_mpbhscexp_2025_bulkrnaseq/-/blob/main/Fig5F_processing_and_plot.R)