Radioactive skin injury single cell analysis method and application thereof

Abstract

The invention provides a radioactive skin injury single cell analysis method and application thereof. Comprising the following steps: sample preparation and single cell sequencing, transcriptome data processing, transcript downstream analysis, cell circulation analysis: calculating circulation fractions by using a Cell CycleScoring function and distributing circulation stages, quasi-timing analysis including subgroup analysis, quasi-timing analysis and difference analysis, RNA velocity analysis including RNA velocity determination and velocity field projection, and ligand and receptor interaction analysis. The method comprises the following steps: carrying out intercellular communication analysis, displaying a result through netVisualbub and rankNet functions, and carrying out enrichment analysis. The invention provides a reliable analysis method for researching the change of cell types and mechanisms at different time points after radiotherapy, so that technicians in the field can more deeply understand the characteristics of different cells after radiotherapy-induced skin injury, which is beneficial to exploring new biomarkers and researching and developing new treatment strategies, and has a wide application prospect. And skin injury caused by radiotherapy is reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a single-cell analysis method for radiation-induced skin damage and an application thereof. Background Art

[0002] Radiation-induced skin injury (RISI) is a notable complication of radiotherapy, mainly affecting sensitive areas such as the face, neck, trunk, and extremities. The severity of RISI is variable, ranging from mild erythema to severe wet desquamation and ulceration, influenced by radiation type and dose factors. Notably, approximately 95% of individuals receiving tumor radiotherapy exhibit varying degrees of radiation dermatitis. Radiation injury can severely affect patients' quality of life and health, potentially leading to treatment interruptions and inadequate follow-up treatment. Although advances in IMRT technology have improved dose distribution, limited understanding of the inflammatory mechanisms and pathogenic triggers of radiation dermatitis has hampered progress in mitigating these adverse side effects. In fact, multiple factors, cells, microbiota, and signaling pathways are involved in the initiation of inflammatory, apoptotic, and oxidative responses during the development of radiation dermatitis. However, the exact role of skin cells, their disrupted mutual communication, and the cellular mediators that drive this inflammatory susceptibility have not been fully characterized. Therefore, a universally accepted standard of care for radiation injury is still lacking.

[0003] There is no single-cell analysis process and application for radiation-induced skin damage in the existing technology, especially the interactive analysis technology for cell subpopulations; therefore, we propose a single-cell analysis method for radiation-induced skin damage and its application. Summary of the invention

[0004] The purpose of the present invention is to provide a single-cell analysis process for radiation-induced skin damage and its application to solve the problems existing in the prior art. The present invention provides a reliable analysis method for studying the changes in cell types and mechanisms at different time points after radiotherapy, so that technicians in this field can have a deeper understanding of the characteristics of different cells after radiotherapy-induced skin damage, which is helpful to explore new biomarkers, research and develop new treatment strategies, and reduce the occurrence of skin damage caused by radiotherapy. In order to explore the cell composition in the skin after radiotherapy, the key analysis mechanisms of fibroblasts and keratinocytes, and determine the key mediators of inflammatory changes associated with the progression of radiation dermatitis, it is very necessary to develop the present invention. The following technical scheme: A single-cell analysis method for radiation-induced skin damage, comprising the following analysis steps:

[0005] Step 1: Sample preparation and single-cell sequencing: The sample preparation and single-cell sequencing process includes sample processing, sample collection, and organization;

[0006] Step 2, transcriptome data processing: the transcriptome data processing includes data mapping and data integration and quality control;

[0007] Step 3, transcript downstream analysis: transcript downstream analysis includes data standardization and batch effect removal as well as dimensionality reduction analysis and cell cluster identification. The cell cluster identification: using the JackStraw function in Seurat to identify significant PCA;

[0008] Step 4: Cell cycle analysis: The cycle score of each cell is calculated using the CellCycleScoring function, and each cycle stage is assigned to metadata.

[0009] Step 5: Pseudo-time series analysis, which includes subgroup analysis, reverse time series analysis, and difference analysis;

[0010] Step 6, RNA velocity analysis, wherein the RNA velocity analysis includes RNA velocity measurement and velocity field projection;

[0011] Step 7, ligand-receptor interaction analysis, including cell-cell communication analysis and displaying the results through netVisual_bubble and rankNet functions;

[0012] Step 8, enrichment analysis, the enrichment analysis includes biological differential expression analysis, and compares the expression levels of pathways between two samples by using the GSEA function; and by single sample enrichment scoring. As a preferred technical solution of the present invention, step 1, sample preparation and single cell sequencing, specifically includes the following steps:

[0013] Step 11, sample treatment: Female C57 / BL / 6 samples received a single dose of 20 Gy of radiation therapy to the right dorsal region, above the right forelimb;

[0014] Step 12, sample collection: on the 7th and 30th days after treatment, collect skin samples; Step 13, tissue dissociation: use the sample liver dissociation kit to process 300 mg of skin tissue and dissect it into small pieces (<1mm 3 ) and transferred to a gentleMACS tube containing an enzyme mix; the tissue was processed using a gentleMACS dissociator, diluted with DMEM and filtered through a 70 μm membrane;

[0015] Step 14, preparation of cell suspension: centrifuge the cell suspension at 300 × g, 4°C for 5 minutes and resuspend in PBS;

[0016] Step 15, Single-cell sequencing: Load the single-cell suspension onto the Chromium Single-cell Controller and generate single-cell gel beads (GEMs) using the Single-cell 3' Library and Gel Bead Kit V3.1; add 10,000 cells to each channel and capture up to 8,000 cells per library; perform reverse transcription reactions to generate barcoded cDNA and RNA, and then sequence each library using the NovaSeq platform.

[0017] As a preferred technical solution of the present invention, step 2, transcriptome data processing: specifically comprises the following steps:

[0018] Step 21. Data mapping: Use Cell Ranger (5.0.0) to map transcripts to the mm10 reference genome and generate a raw gene expression matrix for each scRNA-seq sample.

[0019] Step 22, data integration and quality control: The combined gene expression matrix was converted into a Seurat object using the R software package Seurat (4.3.0); abnormal cells expressing less than 200 or more than 10,000 genes were excluded.

[0020] As a preferred technical solution of the present invention, step 3, downstream analysis of transcripts: comprises the following steps: step 31, data normalization: column normalization and logarithmic transformation of the gene expression matrix; step 32, batch effect removal: using the HarmonyR package to remove batch effects;

[0021] Step 33: Dimensionality reduction analysis: perform principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimensionality reduction.

[0022] As a preferred technical solution of the present invention, step 5, pseudo-time series analysis, includes the following steps: step 51, subpopulation analysis: performing fibroblast subpopulation analysis based on fibroblast-specific markers; step 52, reverse time series analysis: transposing the subpopulation into a reverse time series analysis object, and using Monocle 2 (2.26.0) to simulate the developmental trajectory of the cell;

[0023] Step 53, differential analysis: using DDRTree's dimensionality reduction algorithm to identify pseudo-time dependent genes. As a preferred technical solution of the present invention, step 6, RNA velocity analysis, comprises the following steps: Step 61, RNA velocity determination: determining RNA velocity based on the sheared and unsheared transcript reads in the CellRanger output folder;

[0024] Step 62: Velocity field projection: Calculate spliced ​​and unspliced ​​reads based on the scVelo Python package, apply a stochastic model to estimate RNA velocity, and project the velocity field into UMAP space.

[0025] As a preferred technical solution of the present invention, step 7, ligand-receptor interaction analysis, comprises the following steps:

[0026] Step 71, Cell-cell communication analysis: Use the Cellchat (1.1.3) R package to infer ligand and receptor interactions between the population or subpopulations at different time points.

[0027] As a preferred technical solution of the present invention, step 8, enrichment analysis, includes the following steps: step 81, biological differential expression analysis: using the ClusterProfiler (4.9.0) R package to explore biological themes related to differentially expressed genes at different time points or between different cell types; step 82, single sample enrichment scoring: using the GSVA (1.46.0) R software package to determine the single sample enrichment score, and using the H.ARI.V7.4 gene library for enrichment scoring.

[0028] As a preferred technical solution of the present invention, step 8, enrichment analysis, retains biological processes with p values ​​less than 0.05 in the GO database; and uses the GSEA function to compare the expression levels of pathways between two samples. Application of a single-cell analysis method for radiation-induced skin damage, molecular analysis of radiation-induced skin damage in various radiotherapy stages of mouse data sets.

[0029] Compared with the prior art, the beneficial effects of the present invention are as follows: In the scheme of the present invention: The present invention provides a reliable analysis method for studying the changes in cell types and mechanisms at different time points after radiotherapy, so that technicians in this field can have a deeper understanding of the characteristics of different cells after radiotherapy-induced skin damage, which is helpful to explore new biomarkers, research and develop new treatment strategies, and reduce the occurrence of skin damage caused by radiotherapy. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 A schematic diagram of tissue acquisition provided by the present invention;

[0031] Figure 2 A schematic diagram of the cell types provided by the present invention;

[0032] Figure 3 A schematic diagram of the cell types provided by the present invention;

[0033] Figure 4 A schematic diagram of the cell types provided by the present invention;

[0034] Figure 5 A schematic diagram of analysis data provided by the present invention;

[0035] Figure 6 A schematic diagram of the distribution structure of the samples provided by the present invention;

[0036] Figure 7 This is a schematic diagram of the CellChat analysis and display structure provided by the present invention. DETAILED DESCRIPTION

[0037] To make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be described clearly and completely in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them.

[0038] Therefore, the following detailed description of the embodiments of the present invention is not intended to limit the scope of the claimed invention, but merely represents some embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. It should be noted that the embodiments of the present invention and the features and technical solutions in the embodiments can be combined with each other without conflict. It should be noted that similar numbers and letters represent similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0039] Example 1: Please refer to Figure 1-Figure 7 , a single cell analysis method for radiation-induced skin damage, comprising the following analysis steps:

[0040] Step 1: Sample preparation and single-cell sequencing: The sample preparation and single-cell sequencing process includes sample processing, sample collection, and organization;

[0041] Step 2, transcriptome data processing: the transcriptome data processing includes data mapping and data integration and quality control;

[0042] Step 3, transcript downstream analysis: transcript downstream analysis includes data standardization and batch effect removal as well as dimensionality reduction analysis and cell cluster identification. The cell cluster identification: using the JackStraw function in Seurat to identify significant PCA;

[0043] Step 4: Cell cycle analysis: The cycle score of each cell is calculated using the CellCycleScoring function, and each cycle stage is assigned to metadata.

[0044] Step 5: Pseudo-time series analysis, which includes subgroup analysis, reverse time series analysis, and difference analysis;

[0045] Step 6, RNA velocity analysis, wherein the RNA velocity analysis includes RNA velocity measurement and velocity field projection;

[0046] Step 7, ligand-receptor interaction analysis, including cell-cell communication analysis and displaying the results through netVisual_bubble and rankNet functions;

[0047] Step 8: Enrichment analysis, which includes biological differential expression analysis and compares the expression levels of pathways between two samples by using the GSEA function; and single sample enrichment scoring. Step 1: Sample preparation and single cell sequencing, specifically including the following steps:

[0048] Step 11, sample treatment: Female C57 / BL / 6 samples received a single dose of 20 Gy of radiation therapy to the right dorsal region, above the right forelimb;

[0049] Step 12, sample collection: Collect skin samples on the 7th and 30th days after treatment;

[0050] Step 13, tissue dissociation: 300 mg of skin tissue was processed using a sample liver dissociation kit and dissected into small pieces (<1 mm 3 ) and transferred to a gentleMACS tube containing an enzyme mix; the tissue was processed using a gentleMACS dissociator, diluted with DMEM and filtered through a 70 μm membrane;

[0051] Step 14, preparation of cell suspension: centrifuge the cell suspension at 300 × g, 4°C for 5 minutes and resuspend in PBS;

[0052] Step 15, Single-cell sequencing: Load the single-cell suspension onto the Chromium Single-cell Controller and generate single-cell gel beads (GEMs) using the Single-cell 3' Library and Gel Bead Kit V3.1; add 10,000 cells to each channel and capture up to 8,000 cells per library; perform reverse transcription reactions to generate barcoded cDNA and RNA, and then sequence each library using the NovaSeq platform.

[0053] Step 2: Transcriptome data processing: specifically including the following steps:

[0054] Step 21. Data mapping: Use Cell Ranger (5.0.0) to map transcripts to the mm10 reference genome and generate a raw gene expression matrix for each scRNA-seq sample.

[0055] Step 22, data integration and quality control: The combined gene expression matrix was converted into a Seurat object using the R software package Seurat (4.3.0); abnormal cells expressing less than 200 or more than 10,000 genes were excluded.

[0056] Step 3: Downstream analysis of transcripts: including the following steps:

[0057] Step 31, data normalization: perform column normalization and logarithmic transformation on the gene expression matrix; Step 32, batch effect removal: use the Harmony R package to remove batch effects;

[0058] Step 33: Dimensionality reduction analysis: perform principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimensionality reduction.

[0059] Step 5: Quasi-timing analysis, including the following steps:

[0060] Step 51, subpopulation analysis: fibroblast subpopulation analysis was performed based on fibroblast-specific markers; Step 52, reverse time series analysis: the subpopulation was transposed into a reverse time series analysis object, and Monocle 2 (2.26.0) was used to simulate the developmental trajectory of the cells;

[0061] Step 53: Differential analysis: Use DDRTree's dimensionality reduction algorithm to identify pseudo-time dependent genes. Step 6: RNA velocity analysis, including the following steps:

[0062] Step 61, RNA velocity determination: RNA velocity is determined based on the sheared and unsheared transcript reads in the Cell Ranger output folder;

[0063] Step 62: Velocity field projection: Calculate spliced ​​and unspliced ​​reads based on the scVelo Python package, apply a stochastic model to estimate RNA velocity, and project the velocity field into UMAP space.

[0064] Step 7, ligand-receptor interaction analysis, including the following steps:

[0065] Step 71, Cell-cell communication analysis: Use the Cellchat (1.1.3) R package to infer ligand and receptor interactions between the population or subpopulations at different time points.

[0066] Step 8, enrichment analysis, including the following steps:

[0067] Step 81, biological differential expression analysis: Use ClusterProfiler (4.9.0) R package to explore biological themes related to differentially expressed genes at different time points or between different cell types;

[0068] Step 82, single sample enrichment score: Use the GSVA (1.46.0) R software package to determine the single sample enrichment score, and use the H.ARI.V7.4 gene library for enrichment scoring.

[0069] Step 8: Enrichment analysis, retaining biological processes with p-values ​​less than 0.05 in the GO database; using the GSEA function to compare the expression levels of pathways between the two samples.

[0070] Application of a single-cell analysis method for radiation-induced skin lesions Molecular analysis of radiation-induced skin lesions at various stages of radiotherapy in a mouse dataset.

[0071] Example 2: The present invention provides a single cell analysis method for radiation skin damage. All codes are edited and organized in R language and Python, as follows:

[0072] 1. Single-cell sequencing of mouse skin samples

[0073] Female C57 / BL / 6 mice received a single (20 Gy) radiation treatment targeting the right dorsal region, specifically above the right forelimb. After radiotherapy, mouse skin samples were collected 7 and 30 days after treatment and analyzed using scRNA-seq and multiplex immunofluorescence staining. Approximately 300 mg of skin tissue per group was processed using a mouse liver dissociation kit (130-105-807, Miltenyi Biotec Inc., Auburn, CA, USA). Tissue samples were dissected into small pieces (<1 mm3) and transferred to a gentleMACS tube containing an enzyme mixture. Subsequently, the tissue was processed using a gentleMACS dissociator (Miltenyi Biotec Inc., Auburn, CA, USA). Subsequently, 2 ml of DMEM was added for dilution, and the tissue lysate was filtered with a 70 μm membrane. The cell suspension was centrifuged at 300 × g and 4 ° C for 5 minutes, the supernatant was fully aspirated, the cells were washed, and resuspended with PBS. The single-cell suspension was loaded onto the Chromium single-cell controller. Single cell gel beads (GEMs) were generated using the Single Cell 3' Library and Gel Beads Kit V3.1 (10xGenomics, 60318) and the Chromium Single Cell A Chip Kit (10xGenomics, 60318) using existing protocols. Approximately 10,000 cells were added to each channel, and up to 8,000 cells were captured per library. After cell capture and cleavage, a reverse transcription reaction was performed to generate barcoded cDNA. Subsequently, barcoded RNA was generated by reverse transcription in each individual GEMs. Each library was then sequenced using the NovaSeq platform (Illumina).

[0074] 2. Transcriptome data processing

[0075] Transcripts were mapped to the mm10 reference genome using Cell Ranger (5.0.0) to generate a raw gene expression matrix for each scRNA-seq sample. The combined gene expression matrix was converted into a Seurat object using the R package Seurat (4.3.0). Cells expressing less than 200 or more than 10,000 genes were considered aberrant and excluded.

[0076] Code:

[0077]

[0078] 3. Downstream analysis of transcripts The gene expression matrix was column normalized and logarithmically transformed. To identify cell clusters, a set of highly variable genes was first identified by packing the average expression of all genes into groups of uniform size and calculating the median dispersion in each group. The integrated single-cell object was processed using the NormalizeData function, and the normalization method LogNormalize was selected.

[0079] The R package Harmony was used to remove batch effects. Principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimensionality reduction were then performed. The JackStraw function in Seurat was used to identify significant PCAs.

[0080] Code:

[0081] seuratSet<-NormalizeData(seuratSet,normalization.method="LogNormalize",scale.factor=10000)

[0082] head(seuratSet@assays$RNA@data)

[0083] seuratSet<-FindVariableFeatures(seuratSet,selection.method="vst",nfeatures=3000)

[0084] library(harmony)

[0085] sce.all.int<-seuratSet

[0086] sce.all.int<-RunHarmony(sce.all,c("orig.ident"),

[0087] plot_convergence = T)

[0088] harmony_embeddings<-Embeddings(sce.all.int,'harmony')sce.all.int=RunTSNE(sce.all.int,reduction="harmony",dims=1:20)sce.all.int=RunUMAP(sce.all.int,reduction="harmony",dims=1:20)names(sce.all.int@reductions)

[0089] names(sce.all@reductions)

[0090] p1<-DimPlot(sce.all.int,reduction="umap",group.by="orig.ident",label=T,label.box=T)+NoAxes()p2<-DimPlot(sce.all.int,reduction="umap",group.by="RNA_snn_res.0.5",label=T,label.box=T)+NoAxes()plot(p1+p2)

[0091] seuratSet.harmony<-sce.all.int

[0092] seuratSet.harmony=FindNeighbors(seuratSet.harmony,reduction="harmony",dims=1:20)

[0093] for(res in c(0.1,0.2,0.3,0.5,0.8,1)){

[0094] seuratSet.harmony=FindClusters(seuratSet.harmony,resolution=res,algorithm=1)}

[0095] p1<-DimPlot(seuratSet,reduction="umap",group.by="orig.ident",label=T,label.box=T)+NoAxes()

[0096] p2<-DimPlot(seuratSet,reduction="umap",group.by=

[0097] "RNA_snn_res.0.5",label=T,label.box=T)+NoAxes()plot(p1+p2)

[0098] saveRDS(seuratSet.harmony,'seuratSet.harmony.rds')

[0099] 4. Cell cycle analysis

[0100] The cycling score for each cell was calculated using the CellCycleScoring function, and each cycle stage was assigned to metadata by setting reference genes associated with the cell cycle, including 42 G1 / S and 52 G2M core genes.

[0101] Code:

[0102] mouse_cell_cycle_genes<-readRDS("~ / mouse_cell_cycle_genes.rds")seuratSet.harmony<-CellCycleScoring(seuratSet.harmony,s.features=mouse_cell_cycle_genes$s.genes,g2m.features=

[0103] mouse_cell_cycle_genes$g2m.genes,set.ident=TRUE)

[0104] 5. Quasi-series analysis

[0105] We performed fibroblast subpopulation analysis based on fibroblast-specific markers. Afterwards, the subpopulations were transposed into reverse-chronological analysis objects. Monocle 2 (2.26.0) was used for simulation.

[0106] The developmental trajectories of these cells were analyzed using the dimensionality reduction algorithm of DDRTree, with other parameters set to default. Considering that all genes with q values ​​below 0.1 are pseudo-time-dependent genes, differential analysis found significantly changed pseudo-time-dependent genes.

[0107] Code:

[0108]

[0109]

[0110] 6. RNA velocity analysis

[0111] RNA velocity was determined based on the number of spliced ​​and unspliced ​​transcript reads in the Cell Ranger output folder. Based on the scVelo Python package, spliced ​​and unspliced ​​reads were counted, basic preprocessing was performed, and a stochastic model was applied to estimate RNA velocity. Subsequently, the velocity field was projected into

[0112] The other parameters use the default parameters in the velocyto R package.

[0113] Code:

[0114] #Integrate multiple loom files

[0115] import loompy

[0116] files=[' / CON1.loom',' / CON2.loom',' / CON3.loom',' / XRT1.loom',' / XRT2.loom',

[0117] ' / XRT_6.loom']

[0118] output_filename=' / combined.loom'

[0119] loompy.combine(files,output_filename,key="Accession")

[0120] #Integrate loom files and meta-data data

[0121] import scvelo as scv

[0122] import pandas as pd

[0123] import numpy as np

[0124] import os

[0125] loom_data=scv.read(' / combined.loom',cache=False)

[0126] loom_data.obs.head()

[0127] #Read meta information in seuratSet

[0128] sample_obs=pd.read_csv(os.path.join(meta_path,"cellID_obs.csv"))

[0129] cell_umap=pd.read_csv(os.path.join(meta_path,"cell_embeddings.csv"),header=0,names=["CellID","UMAP_1","UMAP_2"])

[0130] cell_tsne=pd.read_csv(os.path.join(meta_path,"cell_embeddings_tsne.csv"),header=0,names=["CellID","TSNE_1","TSNE_2"])

[0131] cell_clusters=pd.read_csv(os.path.join(meta_path,"cell_clusters.csv"),header=0,names=["CellID","cluster"])

[0132] cell_celltype=pd.read_csv(os.path.join(meta_path,"cell_celltype.csv"),header=0,names=["CellID","celltype"])

[0133] # Take the intersection of the cell file and the RNA shear rate file to retain the cell type of interest

[0134] sample_one=loom_data[np.isin(loom_data.obs.index,sample_obs)]sample_one.obs.head()

[0135] #Build umap coordinates, cluster, celltype information data frame

[0136] sample_one_index=pd.DataFrame(sample_one.obs.index)

[0137] sample_one_index = sample_one_index.rename(columns={0: 'Cell ID'})

[0138] umap_ordered = sample_one_index.merge(cell_umap, on="CellID") umap_ordered.head()

[0139] tsne_ordered = sample_one_index.merge(cell_tsne, on="CellID") tsne_ordered.head()

[0140] celltype_ordered = sample_one_index.merge(cell_celltype, on="CellID")

[0141] celltype_ordered.head()

[0142] clusters_ordered = sample_one_index.merge(cell_clusters, on="CellID")

[0143] clusters_ordered.head()

[0144] # Add umap coordinates and cluster information to sample_one

[0145] umap_ordered = umap_ordered.iloc[:, 1:]

[0146] tsne_ordered = tsne_ordered.iloc[:, 1:]

[0147] clusters_ordered = clusters_ordered.iloc[:, 1:]

[0148] celltype_ordered = celltype_ordered.iloc[:, 1:]

[0149] sample_one.obsm['X_umap'] = umap_ordered.values

[0150] sample_one.obsm['X_tsne']=tsne_ordered.values

[0151] sample_one.uns['clusters']=clusters_ordered.values

[0152] sample_one.obs['celltype']=celltype_ordered.values

[0153] adata=sample_one

[0154] adata.var_names_make_unique()

[0155] #save model to file

[0156] adata.write('Allcelltype_dynamicModel_Fib.h5ad',compression='gzip')

[0157] #read h5ad file

[0158] adata=scv.read('Allcelltype_dynamicModel_Fib.h5ad')

[0159] #Running RNA Velocity

[0160] scv.pp.filter_and_normalize(adata,min_shared_counts=30,n_top_genes=2000)

[0161] scv.pp.moments(adata,n_pcs=30,n_neighbors=30)

[0162] scv.tl.velocity(adata,mode="stochastic")

[0163] scv.tl.velocity_graph(adata)

[0164] scv.pl.velocity_embedding(adata,basis='X_umap',arrow_size=5)

[0165] adata.write('Allcelltype_dynamicModel_Fib_1.h5ad',compression='gzip')

[0166] scv.pl.velocity_embedding_stream(adata,basis='X_tsne',color="celltype")

[0167] scv.pl.velocity_embedding_stream(adata,basis='X_tsne',color="celltype",size=20,alpha=0.8)

[0168] 7. Analysis of ligand-receptor interaction

[0169] To understand intercellular communication, the Cellchat (1.1.3) R package was used to infer the interaction of ligands and receptors based on the total population or subpopulations at different time points. In this study, Cellchat analysis was performed using default parameters to infer the expression changes of specific pathways at different time points. The annotation library references CellChatDB.mouse. The results are mainly displayed by netVisual_bubble and rankNet functions.

[0170] Code:

[0171] cellchat<-createCellChat(object=seuratSet.harmony,group.by="celltype")

[0172] levels(cellchat@idents)

[0173] groupSize<-as.numeric(table(cellchat@idents))groupSize

[0174] #Set up the ligand-receptor interaction database CellChatDB<-CellChatDB.mouse

[0175] showDatabaseCategory(CellChatDB)

[0176] dplyr::glimpse(CellChatDB$interaction)

[0177] CellChatDB.use<-CellChatDB

[0178] cellchat@DB<-CellChatDB.use

[0179] #Preprocess expression data for cell communication analysis cellchat<-subsetData(cellchat)

[0180] future::plan("multiprocess",workers=4)

[0181] cellchat<-identifyOverExpressedGenes(cellchat)cellchat<-identifyOverExpressedInteractions(cellchat)cellchat<-projectData(cellchat,PPI.mouse)

[0182] #Calculate the communication probability and infer the cellchat network cellchat<-computeCommunProb(cellchat)

[0183] cellchat<-filterCommunication(cellchat,min.cells=10)#Extract the inferred cellchat network as a data frame df.net<-subsetCommunication(cellchat)

[0184] # Infer cell-cell communication at the pathway level cellchat <- computeCommunProbPathway(cellchat) cellchat <- aggregateNet(cellchat)

[0185] groupSize<-as.numeric(table(cellchat@idents))netVisual_circle(cellchat@net$count,vertex.weight=groupSize,weight.scale=T,label.edge=F,title.name="Number of interactions")netVisual_circle(cellchat@net$weight,vertex.weight=groupSize,weight.scale=T,label.edge=F,title.name="Interactionweights / strength")

[0186] 8. Enrichment analysis

[0187] ClusterProfiler (4.9.0) R package was used to explore biological themes related to differentially expressed genes at different time points or between different cell types. Clusterprofiler is widely used to identify the biological functions of genes, including GO and KEGG. In this study, biological processes with p-values ​​less than 0.05 in the GO database were mainly retained. The GSEA function was used to compare the expression levels of pathways between two samples with default parameters. The GSVA (1.46.0) R package was used to determine the single sample enrichment score, and the H.ARI.V7.4 gene library was used for enrichment scoring.

[0188] Code:

[0189]

[0190]

[0191] Experimental example:

[0192] Molecular analysis of radiation-induced skin damage in mouse datasets at various radiotherapy stages 1. Female C57 / BL / 6 mice received a single (20 Gy) radiation therapy targeting the right dorsal region, specifically above the right forelimb. After radiotherapy, mouse skin samples were collected 7 and 30 days after treatment and analyzed using scRNA-seq and multiplex immunofluorescence staining. The study was divided into three groups: day 0, day 7, and day 30 after radiotherapy.

[0193] 2. Through transcript data processing and downstream dimensionality reduction clustering analysis, a data set containing 71,412 cells was finally formed. Using UMAP for visualization, 10 major cell types were identified ( Figure 2C). All 10 clusters included cells from skin biopsies of control and irradiated mice (days 7 and 30). Figure 2 E). Cell composition analysis showed that the proportion of keratinocytes (KCs) and T cells increased at day 7 after radiation and further amplified at day 30 ( Figure 3 F). This suggests a potential pathogenic role for these cell types in RISI. Differential gene expression analysis, identification of cell type specific marker genes and visualization of the top 50 ranked genes for each cell type in a heatmap ( Figure 4 H); 3. The role of circulating keratinocytes in the BMP signaling pathway in wound healing and fibrosis in RISI. GO functional enrichment analysis revealed the upregulation of WNT and BMP-related signaling pathways ( Figure 5 B). Given that BMPs belong to the transforming growth factor β (TGF-β) superfamily, and TGF-β activates fibrosis through the Smad gene, gene expression analysis showed upregulation of Smad in KCs after radiation ( Figure 5 C). In addition, Smad expression and trajectory analysis revealed an increase in Smad levels over time ( Figure 5 D and 5E). These findings suggest that fibrosis develops gradually after RT.

[0194] 4. RISI terminally differentiated secretory papillary fibroblasts with immune cell recruitment ability. Gene differential analysis showed that in addition to genes related to ribosomes and mitochondria (such as Rpl35, Rps28 and mt-Co1), a large number of genes related to metabolism, inflammation and immune response (such as Nr1d1 and Bgn) were upregulated after radiotherapy ( Figure 6 A). This indicates that radiation exposure not only affects the metabolism of fibroblasts, but also affects immune subsets, ultimately triggering an inflammatory response within fibroblasts. Compared with day 7, the expression of collagen-related genes, such as COL15A1 and COL16A1, increased in samples on day 30.

[0195] Fibroblasts isolated from aggregated cell populations were classified into six different subpopulations ( Figure 6 B), the distribution of samples is shown in the UMAP diagram ( Figure 6 C). Through classic marker annotation, 6 subpopulations were annotated as different fibroblast subpopulations ( Figure 6 D and 6E), in addition, the composition of fibroblasts undergoes rapid and dramatic changes after radiotherapy ( Figure 6 F).

[0196] Enrichment analysis was performed to explore the role of fibroblast subsets in the pathogenesis of RISI.

[0197] The PI3K / AKT / mTOR signaling pathway is a hallmark of inflammatory mediators and is significantly associated with secretory-papillary fibroblasts ( Figure 6 G). Analysis showed that secretory papillary fibroblasts have a significant immune activation function, as evidenced by upregulated genes ( Figure 6 H) is proven.

[0198] 5. Development of regulatory networks between keratinocytes (KCs), fibroblasts (Fibs), and endothelial cell (ECs) subsets and immune cells in RISI.

[0199] Considering the key role of immune cells in RISI, analysis of intercellular communication is necessary. CellChat analysis showed that as radiation dermatitis progressed, there was an increase in ligand-receptor pair interactions between circulating keratinocytes, secretory papillary fibroblasts, lymphatic endothelial cells, and immune cells ( Figure 7 A). Extract specific cell subsets and delve into their complex interaction signals with immune cells. Upregulation of the Wnt pathway in circulating keratinocytes after radiotherapy is primarily driven by Wnt4 and Wnt10a ( Figure 7 B). In contrast, after an initial increase in BMP signaling in other subsets, interactions with specific immune cells, including blood Tcm, anti-inflammatory macrophages, and Trm2, were reduced, with a particularly pronounced deficiency observed in lymphatic endothelial cells ( Figure 7 C). In addition, the results showed that Gas6, an upregulated gene associated with aging, can bind to ligands of various immune cells ( Figure 7 D). In the context of skin diseases, studies have shown that the chemokine CXCL8 produced by M2 macrophages can interact with endothelial cells expressing ACKR1. Analysis results showed that after radiation therapy, there were enhanced interactions between KCs, Fibs, ECs subsets and immune cells involving chemokines (including CXCL1, CXCL2, CXCL10, CCL7 and CCL2) and Ackr1. These strong interactions suggest the important involvement of these subsets in recruiting leukocytes to the skin during RISI ( Figure 7 E). Overall, these specific receptor-ligand communications further emphasize the critical contributions of circulating keratinocytes, secretory papillary fibroblasts, and lymphatic endothelial cells in the development of RISI.

[0200] The above embodiments are only used to illustrate the present invention and are not intended to limit the technical solutions described in the present invention. Although the present invention has been described in detail with reference to the above embodiments, the present invention is not limited to the above specific implementation methods. Therefore, any modification or replacement of the present invention; and all technical solutions and improvements thereof that do not depart from the spirit and scope of the invention are included in the scope of the claims of the present invention.

Claims

1. A single cell analysis method for radiation-induced skin damage, characterized in that: The analysis steps include: Step 1: Sample preparation and single-cell sequencing: The sample preparation and single-cell sequencing process includes sample processing, sample collection, and organization; Step 2, transcriptome data processing: the transcriptome data processing includes data mapping and data integration and quality control; Step 3, transcript downstream analysis: transcript downstream analysis includes data standardization and batch effect removal as well as dimensionality reduction analysis and cell cluster identification. The cell cluster identification: using the JackStraw function in Seurat to identify significant PCA; Step 4: Cell cycle analysis: The cycle score of each cell is calculated using the CellCycleScoring function, and each cycle stage is assigned to metadata. Step 5: Pseudo-time series analysis, which includes subgroup analysis, reverse time series analysis, and difference analysis; Step 6, RNA velocity analysis, wherein the RNA velocity analysis includes RNA velocity measurement and velocity field projection; Step 7, ligand-receptor interaction analysis, including cell-cell communication analysis and displaying the results through netVisual_bubble and rankNet functions; Step 8: Enrichment analysis, which includes biological differential expression analysis and compares the expression levels of pathways between two samples by using the GSEA function; and by single sample enrichment scoring.

2. A single cell analysis method for radiation skin damage according to claim 1, characterized in that: Step 1: Sample preparation and single-cell sequencing, specifically including the following steps: Step 11, sample treatment: Female C57 / BL / 6 samples received a single dose of 20 Gy of radiation therapy to the right dorsal region, above the right forelimb; Step 12, sample collection: Collect skin samples on the 7th and 30th days after treatment; Step 13, tissue dissociation: 300 mg of skin tissue was processed using a sample liver dissociation kit and dissected into small pieces (<1 mm 3 ) and transferred to a gentleMACS tube containing an enzyme mix; the tissue was processed using a gentleMACS dissociator, diluted with DMEM and filtered through a 70 μm membrane; Step 14, preparation of cell suspension: centrifuge the cell suspension at 300 × g, 4°C for 5 minutes and resuspend in PBS; Step 15, Single-cell sequencing: Load the single-cell suspension onto the Chromium Single-cell Controller and generate single-cell gel beads (GEMs) using the Single-cell 3' Library and Gel Bead Kit V3.1; add 10,000 cells to each channel and capture up to 8,000 cells per library; perform reverse transcription reactions to generate barcoded cDNA and RNA, and then sequence each library using the NovaSeq platform.

3. A single cell analysis method for radiation skin damage according to claim 2, characterized in that: Step 2: Transcriptome data processing: specifically including the following steps: Step 21. Data mapping: Use Cell Ranger (5.0.0) to map transcripts to the mm10 reference genome and generate a raw gene expression matrix for each scRNA-seq sample. Step 22, data integration and quality control: The combined gene expression matrix was converted into a Seurat object using the R software package Seurat (4.3.0); abnormal cells expressing less than 200 or more than 10,000 genes were excluded.

4. A single cell analysis method for radiation skin damage according to claim 3, characterized in that: Step 3: Downstream analysis of transcripts: including the following steps: Step 31, data standardization: column normalization and logarithmic transformation of the gene expression matrix; Step 32: Batch effect removal: Use the Harmony R package to remove batch effects; Step 33: Dimensionality reduction analysis: perform principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimensionality reduction.

5. A single cell analysis method for radiation skin damage according to claim 4, characterized in that: Step 5: Quasi-timing analysis, including the following steps: Step 51, subpopulation analysis: performing fibroblast subpopulation analysis based on fibroblast-specific markers; Step 52, reverse time analysis: transpose the subpopulation into a reverse time analysis object, and use Monocle 2 (2.26.0) to simulate the developmental trajectory of cells; Step 53, differential analysis: Use DDRTree's dimensionality reduction algorithm to identify pseudo-time-dependent genes.

6. A single cell analysis method for radiation skin damage according to claim 5, characterized in that: Step 6, RNA velocity analysis, includes the following steps: Step 61, RNA velocity determination: RNA velocity is determined based on the sheared and unsheared transcript reads in the Cell Ranger output folder; Step 62: Velocity field projection: Calculate spliced ​​and unspliced ​​reads based on the scVelo Python package, apply a stochastic model to estimate RNA velocity, and project the velocity field into UMAP space.

7. A single cell analysis method for radiation-induced skin damage according to claim 6, characterized in that: Step 7, ligand-receptor interaction analysis, including the following steps: Step 71, Cell-cell communication analysis: Use the Cellchat (1.1.3) R package to infer ligand and receptor interactions between the population or subpopulations at different time points.

8. A single cell analysis method for radiation skin damage according to claim 7, characterized in that: Step 8: Enrichment analysis, including the following steps: Step 81, biological differential expression analysis: Use ClusterProfiler (4.9.0) R package to explore biological themes related to differentially expressed genes at different time points or between different cell types; Step 82, single sample enrichment score: Use the GSVA (1.46.0) R software package to determine the single sample enrichment score, and use the H.ARI.V7.4 gene library for enrichment scoring.

9. A single cell analysis method for radiation-induced skin damage according to claim 8, characterized in that: Step 8: Enrichment analysis, retaining biological processes with p-values ​​less than 0.05 in the GO database; using the GSEA function to compare the expression levels of pathways between the two samples.

10. An application of a single cell analysis method for radiation skin damage, characterized in that: Molecular analysis of radiation-induced skin damage across various radiotherapy phases in a mouse dataset.

Images