Clinical specimens and patient features
Clinical specimens collected and analyzed in this study were harvested through informed consent. The study protocol was approved by the ethics committees of Peking University People’s Hospital and Beijing Children's Hospital of Capital Medical University. The clinical characteristics of enrolled patients are listed in Table 1.
Table 1
Demographic and clinical characteristics of the enrolled patients
No.
|
Age/Sex
|
Diagnosis
|
Site
|
Stage
|
Tumor thrombus site
|
Chemotherapy
|
Radiotherapy
|
Operation Time (mins)
|
Estimated blood loss (ml)
|
Follow-up time
|
Outcome
|
1
|
26/M
|
Osteosarcoma
|
Left ilium
|
III
|
Iliac vein
|
Yes
|
No
|
160
|
5000
|
11
|
DOD
|
2
|
26/M
|
Osteosarcoma
|
Left humerus
|
III
|
Axillary vein
|
Yes
|
No
|
215
|
900
|
9
|
DOD
|
3
|
17/F
|
Osteosarcoma
|
Pelvis
|
III
|
Inferior vena cava
|
Yes
|
No
|
230
|
1500
|
22
|
DOD
|
4
|
20/M
|
Osteosarcoma
|
Pelvis
|
III
|
Inferior vena cava
|
Yes
|
No
|
220
|
5100
|
20
|
DOD
|
5
|
24/M
|
Osteosarcoma
|
Left ilium
|
III
|
Iliac vein
|
Yes
|
No
|
225
|
3500
|
7
|
DOD
|
6
|
9/F
|
Osteosarcoma
|
Right femur
|
III
|
Femoral vein
|
Yes
|
No
|
210
|
500
|
13
|
DOD
|
7
|
18/M
|
Osteosarcoma
|
Pelvis
|
III
|
Iliac vein
|
Yes
|
No
|
195
|
1500
|
10
|
DOD
|
DOD, dead of disease. |
Bulk RNA-seq data analysis
We obtained clean data with the adapter removed and low-quality reads filtered. The paired-end reads were mapped to the human genome (UCSC hg38) by Hisat2 (version 2.1.0). Gene annotation and gene expression reads-count were performed by HTSeq (version 0.11.2). Then we selected the protein-coding genes and calculated the normalized expression value using the DESeq2 package (version 1.32.0) on the R platform (version 4.1.0).
Gene-set enrichment analysis (GSEA)
The normalized expression matrix was uploaded to the gene set enrichment analysis (GSEA) software (version 4.1.0). the hallmark gene sets provided by the Molecular Signatures Database (MSigDB) were selected as a reference database. The GSEA program was run with 1,000 permutations for statistical significance estimation, and the default signal-to-noise metric between the two phenotypes was used to rank all genes.
Ingenuity Pathway Analysis (IPA)
Differentially expressed genes (DEGs) were analyzed by the DESeq2 package. DEGs were judged with the threshold adjusted p-value less than 0.05 and absolute log2 fold-change greater than 1. Then these DEGs were uploaded to the Ingenuity Pathway Analysis (IPA) software. Canonical pathways, causal networks, and upstream regulators were enriched with default parameters in this software.
Sample collection and preparation of single-cell suspensions
Obtained during surgery, fresh tumor and thrombus tissues were stored in Dulbecco’s modified Eagle’s medium (DMEM, Gibco) with 10% fetal bovine serum (FBS, Gibco) and processed on ice. Subsequently, the tissues were rinsed with cold phosphate buffer saline (PBS, Gibco) three times and minced into ~1mm3 pieces. The tissue pieces were digested into single-cell suspensions using collagenase-2 (1 mg/mL) at 37℃ for 45min and filtered through a 70-μm cell strainer. Upon digestion, the single-cell suspensions were centrifugated at 400g for 5min and the supernatant was discarded. Adding red blood cell lysis buffer (Solarbio) to remove red blood cells, the suspensions were then passed through a 40-μm filter and centrifugated at 500g for 5min, and resuspended in PBS. Cell viability was validated under the phase-contrast light microscope (Leica) after staining with trypan blue (ACMEC).
Single-cell transcriptome library preparation and sequencing
After resection, tissue specimens were rapidly processed for single-cell RNA sequencing. Single-cell suspensions were prepared according to the protocol of Chromium™ Single Cell 3' Solution (V2 chemistry). All specimens were washed twice with cold 1× PBS. A hemocytometer (Thermo) was used to evaluate cell viability rates. Then, we used Countess (Thermo) to determine the concentration of the single-cell suspension and adjust the concentration to 1000 cells/μl. Samples with a cell concentration lower than that defined in the user guide (i.e., <400 cells/µl) were pelleted and resuspended in a reduced volume, and the concentration of the new solution was determined again. Finally, the cells in each sample were loaded, and libraries were constructed using a Chromium single-cell kit. Single-cell libraries were subjected to 150-bp paired-end sequencing on the Illumina NavoSeq6000 platform.
Single-cell RNA-seq data preprocessing and quality control
After obtaining the paired-end raw reads, we used CellRanger (10x Genomics, v3.1.0) to preprocess the single-cell RNA-seq data. The cell barcodes and unique molecular identifiers (UMIs) of the library were extracted from read1. Then, the reads were split according to their cell (barcode) IDs, and the UMI sequences from read2 were simultaneously recorded for each cell. Quality control of these raw readings was subsequently performed to eliminate adapter contamination, duplicates, and low-quality bases. After filtering barcodes and low-quality readings that were not related to cells, we used STAR (version 2.5.1b) to map the clean reads to the human genome (hg19), and we retained the uniquely mapped reads for UMI counts. Next, we estimated the molecular counts and generated a UMI count matrix for each cell by counting the UMIs for each sample. Finally, we generated a gene-barcode matrix that showed the barcoded cells and gene expression counts. Based on the number of total reads, the number of detected gene features, and the percentage of mitochondrial genes, we performed quality control filtering through Seurat to discard low-quality cells (v3.1.5). Briefly, the proportion of mitochondrial genes inside one cell was calculated to be lower than 10%, and the number of total reads in one cell was below 50000. Additionally, the cells were further filtered according to the following criteria: each sample with no more than 5000 genes detected, respectively, and at least 200 genes detected per cell in any sample. Low-quality cells and outliers were discarded, and the single cells that passed the QC criteria were used for downstream analyses. Doublets were predicted by DoubletFinder (v2.0).
Clustering analysis and cell phenotype recognition
The Seurat software package was used to perform cell clustering analysis to identify major cell types. All Seurat objects constructed from the filtered UMI-based gene expression matrixes of given samples were merged. We first applied the “SCTransform” function to implement normalization, variance stabilization, and feature selection through a regularized negative binomial model. Then, we reduced dimensionality through principal component analysis (PCA). According to standard steps implemented in Seurat, highly variable numbers of principal components (PCs) 1 to 20 were selected and used for clustering using the Uniform Manifold Approximation and Projection (UMAP). We identified the cell types of these groups based on the expression of canonical cell type markers. Finally, the 5 groups of T cells and 2 groups of ADSCs were respectively clustered for downstream analysis.
DEG analysis
The cell-type-specific genes were identified by running Seurat with the ‘FindAllMarkers’ function on a log-transformed expression matrix with the following parameter settings: min.pct=0.25, logfc.threshold=0.25 (that is, there is at least a 0.25 log-scale fold change between the cells inside and outside a cluster), and only.pos =TRUE (that is, only positive markers are returned). For heatmap and violin plots, the SCT-transformed data from the Seurat pipeline were used. Using the Seurat ‘FindMarkers’ function, we found differentially expressed genes (DEGs) between the two cell types. We also used the R package clusterProfiler with the default parameters to identify gene sets that exhibited significant and consistent differences between two given biological states.
Multiplexed immunofluorescence (mIF) staining
To validate the immunostimulatory microenvironment in the tumor thrombus of OS, mIF staining was performed using PANO 7-plex IHC kit (Cat# 0004100100, Panovue). CD3 (Cat# ab135372, Abcam), CD4 (Cat# ZM0418, ZSGB-BIO), CD8A (Cat# CST70306, CST), CD68 (Cat# CST76437, CST), CCL4 (Cat# ab235961, Abcam) and DAPI (Cat# D9542, Sigma-Aldrich) antibodies were sequentially applied, followed by horseradish peroxidase-conjugated secondary antibody incubation and tyramide signal amplification. Next, the slides were microwave heat-treated after each tyramide signal amplification operation. Nuclei were stained with DAPI after all of the antigens above had been labeled. To obtain multispectral images, the stained slides were scanned using the Mantra System (PerkinElmer), which captures the fluorescent spectra at 20-nm wavelength intervals from 420 to 720 nm with identical exposure time. 5 fields of immune cell enriched tumoral area for each slide were selected for image capture. The selected field were scanned to obtain multispectral images using the Mantra System, which captures the fluorescent spectra at 20-nm wavelength intervals from 420 to 720 nm with identical exposure time.
Quantification and statistical analysis
The specific tests used to analyze each set of experiments are indicated in the figure legends. Comparisons between two groups after IHC staining were performed using a two-tailed Student’s t-test, and comparisons among three or more groups were performed using one-way ANOVA. The significance of numerical values of clinical features was determined using unpaired multiple T-test (discovery determined using the two-stage linear step-up procedure of Benjamini, Krieger, and Yekutieli, with Q = 1%). Each row was analysed individually, without assuming a consistent SD. Numerical values were corrected for multiple comparisons using the Holm-Sidak method (alpha= 0.05, presented with padj). All statistical calculations were performed using GraphPad Prism software (ver. 7.0, GraphPad, USA) or R software (https://www.r-project.org/).