I was pointed to a folder with data on our shared PetaLibrary. I found FASTQ files in 5 different subfolders. The RNA-seq modalities that were requested from the sequencing core were: pooled single-cell RNA-seq, dissociated bulk tissue polyA+ bulk sequenced, dissociated bulk tissue rRNA- bulk sequenced, and bulk tissue rRNA- bulk sequenced.
The subfolder and its contents with pooled scRNA-seq were clearly labelled. That leaves 4 different subfolders with unknown modalities. We label them by the dates that the data was generated: 230414, 230418, 230509, and 230626.
Scrutinizing the contents of all the subfolders, we found that folders 230414 and 230418 contained data from disjoint sets of samples (i.e. all different samples),
and likewise 230509 and 230626 contained data from disjoint sets of samples.
Comparing folders 230414 and 230418 with 230509 and 230626 reveals that the samples in 230414 and 230418 are a subset of the samples in 230509 and 230626 (i.e. all samples in 230414 and 230418 are represented in 230509 and 230626).
Folders 230414 and 230418 both contained README noting that the data in that folder was from dissociated cells that were poly-A captured and bulk sequenced.
All the dependencies are provided in a yml file data_detective.yml that can be
built, updated, and activated as follows
conda env create -f data_detective.yml
conda env update -f data_detective.yml
conda activate data_detective
The sequences were aligned using STAR and reference genome refdata-gex-GRCh38-2024-A that was downloaded from the 10x Genomics website.
00_create_STAR_index.shcreates an index for the reference genome01_align_samples.shaligns transcriptomic samples with a reference genome.
In deseq_analysis.ipynb, we're using DEA on the dataset to try and figure out whether the various subsets of data were created using poly A selection or rRNA depletion, and ideally whether the subsets contain data from dissociated bulk or bulk chunks.
As input there is a folder STAR (from the sequence alignment step) with subfolders reflecting the PL organizational structure and each of those subfolders containing folder named with the sample id.
E.g. STAR/230414/2018.
Each of the sample id folders contains (among other things) a ReadsPerGene.out.tab file.
A count matrix is prepared by going through all of the ReadsPerGene files, and a metadata table is created with the folder names.
From the reference genome previously used to align the sequences, we have a genes.gtf file that helps to convert Ensembl gene id to gene names.
Red blood cell (erythrocyte) and adipocyte marker genes are used to test our hypotheses. These can be found:
- Adipocyte marker genes: Emont et al. 2022 (SI Table 1)
- RBC marker genes: Zie et al. 2021 (SI)
We conclude that 230414 and 230418 contain poly A+ dissociated bulk and 230509 and 230626 contain rRNA- bulk.
From the pilot study paper, we found the data GEO: GSE217517. On GEO there are zipped .tsv files contain the information needed to build a gene expression matrix. This study has bulk chunk rRNA-, bulk dissociated rRNA-, and bulk dissociated poly A+ data for 8 samples each.
We originally tried to do DEA to directly compare our current data with the pilot study data, but the analysis was dominated by batch effects.
Instead, in deseq_pilot_study_comp.ipynb we're redoing the analysis from deseq_analysis.ipynb and comparing the results of DEA on just the current data with the results of DEA on just the pilot study data.
This leads us to conclude that the 230414 and 230418 subfolders contain data for poly A+ dissociated bulk tissue, and the 230509 and 230626 subfolders contain data for rRNA- bulk chunk tissue.
The conclusions for 230414 and 230418 match the README associated with those folders. The other folders do not have README's.
To compliment the DEA, we're looking at PCA plots in pca_normalization.ipynb to check for clustering and separation between the different modalities.
The data is normalized first to account for the varying number of reads in the samples, and then PCA is applied.
No clearly separated clusters emerged, but the datapoints from the two different RNA-seq modalities are generally grouped by modality, and it supports our conclusions.