I downloaded a new version of cellranger along with the various reference files provided by 10x for the VD(J) references etc. I got a bit distracted by the pipeline language implemented by 10x called ‘martian’. I have the feeling that it might prove a good thing to play with.
Here are the commands I ran to separate the samples and perform the alignments. There are 4 sample names and each was done with one run of the ‘normal’ GEX scRNASeq method and one of the (new to me) V(D)J library.
April kindly sent some information from 10x which shows that I should have used the multi pipline when preprocessing the data.
Intra-Muscular vs. Nasal
I wrote 4 separate configuration csv files using the templates I downloaded and following a little reading. It seemed to me that I should be able to process them all as a single csv file, but when I attempted that, cellranger did not react well. It also took a few tries before I got the various reference/library options correct.
Note that once cellranger successfully ran for the samples I moved them to the multi/ directory so that I can compare the outputs to when I simply did the ‘count’ operation.
The following invocations of cellranger all appear to work without any problems. Ideally I would like them to be done in a single run, though.
My attempts so far to use the csv configuration to concatenate multiple vdj libraries have not worked, so I chose to do it the stupid way, which is what I should have just done to begin with. Caveat, it works fine for the gex libraries to do it the way the documentation suggests.
cd preprocessing
for i in R1 R2; do
for j in Control Mock_Mex09 IM_Mex09 IN_Mex09; do
A_file=$(/bin/ls A_${j}_VDJ*_${i}_001.fastq.gz)
B_file=$(/bin/ls B_${j}_VDJ*_${i}_001.fastq.gz)
out_file="Concat_${j}_VDJ_${i}.fastq.gz"
cp_cmd="cp ${A_file} ${out_file}"
echo "Running: ${cp_cmd}."
eval $cp_cmd
cat_cmd="cat ${B_file} >> ${out_file}"
echo "Running: ${cat_cmd}."
eval $cat_cmd
done
done
module add cellranger
cellranger multi --id control --csv sample_sheets/multi_config_try05_control.csv
cellranger multi --id mock --csv sample_sheets/multi_config_try05_mock.csv
cellranger multi --id m --csv sample_sheets/multi_config_try05_m.csv
cellranger multi --id n --csv sample_sheets/multi_config_try05_n.csv
mv control mock m n 01multi_combined/
I wonder if I can put the gene annotations into the misc slot of the seurat data structure? And perhaps overload fData() to use it?
<- load_biomart_annotations()$annotation annotations
## The biomart annotations file already exists, loading from it.
<- unique(annotations[, c("hgnc_symbol", "description", "ensembl_gene_id")]) brief
#prefix <- "multi"
<- "01multi_combined" prefix
The following block is mostly a cut/paste of itself where I set the (over)simplified name of each sample. This then becomes the template for the path and parameters used to read the data, create a seurat object, and add the clonotype data from the vdj run.
I wrote a little function to make loading the Seurat data from a sample sheet easier. My intention is to have some of this code write back to that sample sheet.
In case I decide to play with the 4 samples separately, I immediately turned around and extracted all cells from each sample into separate data structures named control_scd, mock_scd, muscular_scd, and nasal_scd (I am using scd as a suffix to denote single-cell datastructure).
<- create_scd("sample_sheets/all_samples_202301.xlsx",
all_scd expression_column = "dataroot",
vdj_t_column = "vdjtcells")
## Did not find the batch column in the sample sheet.
## Filling it in as undefined.
## Warning in CheckDuplicateCellNames(object.list = objects): Some cell names are
## duplicated across objects provided. Renaming to enforce unique cell names.
## Warning in CheckDuplicateCellNames(object.list = objects): Some cell names are
## duplicated across objects provided. Renaming to enforce unique cell names.
I want to take a couple minutes to add some annotations to the seurat object, notably I want to state the identity relationships with some sort of name.
Thus I will make a vector of the the sample IDs and for each one make a category of self/not-self. Note that Seurat comes with a function ‘FindAllMarkers()’ which compares each self to all other samples, so this may be redundant; but it is kind of nice to be able to see the categories as a set of binary indexes.
<- add_binary_states(all_scd) all_scd
I recorded into the xlsx sample sheet the number of cells with defined clonotypes according to cellranger. We should be able to immediately recapitulate those numbers by asking how many cells have a clonotype_id.
<- !is.na(all_scd[["raw_clonotype_id"]])
all_clono summary(all_clono)
## raw_clonotype_id
## Mode :logical
## FALSE:39253
## TRUE :12341
<- !is.na(all_scd[["raw_clonotype_id"]]) &
control_clono "orig.ident"]] == "control"
all_scd[[summary(control_clono)
## raw_clonotype_id
## Mode :logical
## FALSE:49623
## TRUE :1971
<- !is.na(all_scd[["raw_clonotype_id"]]) &
mock_clono "orig.ident"]] == "mock"
all_scd[[summary(mock_clono)
## raw_clonotype_id
## Mode :logical
## FALSE:48504
## TRUE :3090
<- !is.na(all_scd[["raw_clonotype_id"]]) &
m_clono "orig.ident"]] == "m"
all_scd[[summary(m_clono)
## raw_clonotype_id
## Mode :logical
## FALSE:47377
## TRUE :4217
<- !is.na(all_scd[["raw_clonotype_id"]]) &
n_clono "orig.ident"]] == "n"
all_scd[[summary(n_clono)
## raw_clonotype_id
## Mode :logical
## FALSE:48531
## TRUE :3063
The numbers above correspond to the column ‘vdjt_cells’ in the excel sample sheet.
We have some simple functions which can record some potentially useful simple metrics into the metadata of the single cell datastructure. The record_seurat_samples appends to our input samplesheet extra columns containing metrics which might be of interest. In the following example it adds a column of how many cells Seurat thinks there are. This should match what we observed from cellranger as well as the binary sorter above.
Following that, we count up the number of genes observed followed by mean values for:
For a little extra fun, if the verbose parameter is left TRUE, this function uses the new(to me) skim function, which provides a pretty summary of the category of interest. It therefore also makes it easy, if we want, to add SD or quartiles in addition to mean values.
<- is.na(all_scd[["raw_clonotype_id"]])
not_vdj_idx "tcell"]] <- "Yes"
all_scd[[@meta.data[not_vdj_idx, "tcell"] <- "No" all_scd
<- record_seurat_samples(all_scd, type="num_cells") %>%
all_scd record_seurat_samples(type="nFeature_RNA") %>%
record_seurat_samples(type="nCount_RNA") %>%
record_seurat_samples(type="reads", column_name="clonotype_reads") %>%
record_seurat_samples(type="pct_mito", pattern="^mt-") %>%
record_seurat_samples(type="pct_ribo", pattern="^Rp[sl]")
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 51594
## Number of columns 46
## _______________________
## Column type frequency:
## numeric 1
## ________________________
## Group variables orig.ident
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable orig.ident n_missing complete_rate mean sd p0 p25 p50
## 1 nFeature_RNA control 0 1 1999. 1113. 20 1158 1700.
## 2 nFeature_RNA m 0 1 2066. 1275. 22 1101 1651
## 3 nFeature_RNA mock 0 1 2086. 1218. 18 1149 1734
## 4 nFeature_RNA n 0 1 2224. 1416. 27 1063. 1734.
## p75 p100 hist
## 1 2637. 7532 ▇▇▃▁▁
## 2 2778 7590 ▇▆▂▁▁
## 3 2785 7442 ▇▇▃▁▁
## 4 3228 8128 ▇▅▃▁▁
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 51594
## Number of columns 46
## _______________________
## Column type frequency:
## numeric 1
## ________________________
## Group variables orig.ident
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable orig.ident n_missing complete_rate mean sd p0 p25 p50
## 1 nCount_RNA control 0 1 5115. 4137. 500 2149 3618
## 2 nCount_RNA m 0 1 5729. 5383. 500 2155 3601
## 3 nCount_RNA mock 0 1 5539. 4862. 500 2144 3665
## 4 nCount_RNA n 0 1 6488. 6052. 500 2094. 3912
## p75 p100 hist
## 1 7057. 40853 ▇▂▁▁▁
## 2 7642 53707 ▇▁▁▁▁
## 3 7518 44494 ▇▂▁▁▁
## 4 9417. 53836 ▇▂▁▁▁
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 51594
## Number of columns 46
## _______________________
## Column type frequency:
## numeric 1
## ________________________
## Group variables orig.ident
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable orig.ident n_missing complete_rate mean sd p0 p25 p50
## 1 reads control 13209 0.130 7066. 6037. 41 3114. 5365
## 2 reads m 9664 0.304 3685. 4114. 24 1403 2565
## 3 reads mock 10835 0.222 4709. 4085. 42 1998. 3618
## 4 reads n 5545 0.356 3460. 3604. 25 1377 2475
## p75 p100 hist
## 1 9344. 53188 ▇▂▁▁▁
## 2 4469 48039 ▇▁▁▁▁
## 3 6198. 40634 ▇▁▁▁▁
## 4 4232. 36105 ▇▁▁▁▁
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 51594
## Number of columns 46
## _______________________
## Column type frequency:
## numeric 1
## ________________________
## Group variables orig.ident
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable orig.ident n_missing complete_rate mean sd p0 p25 p50 p75
## 1 pct_mito control 0 1 2.47 3.79 0 1.35 1.99 2.92
## 2 pct_mito m 0 1 2.32 3.42 0 1.35 1.97 2.75
## 3 pct_mito mock 0 1 2.10 2.53 0 1.19 1.77 2.55
## 4 pct_mito n 0 1 2.06 3.36 0 1.21 1.76 2.45
## p100 hist
## 1 98.7 ▇▁▁▁▁
## 2 98.4 ▇▁▁▁▁
## 3 98.9 ▇▁▁▁▁
## 4 97.3 ▇▁▁▁▁
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 51594
## Number of columns 46
## _______________________
## Column type frequency:
## numeric 1
## ________________________
## Group variables orig.ident
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable orig.ident n_missing complete_rate mean sd p0 p25 p50 p75
## 1 pct_ribo control 0 1 10.0 7.39 0 5.48 7.69 10.9
## 2 pct_ribo m 0 1 13.3 10.7 0 5.54 8.34 18.9
## 3 pct_ribo mock 0 1 9.08 7.41 0 4.44 6.10 10.6
## 4 pct_ribo n 0 1 13.7 10.5 0 5.82 8.31 20.8
## p100 hist
## 1 51.4 ▇▂▁▁▁
## 2 50.1 ▇▂▂▂▁
## 3 45.0 ▇▂▁▁▁
## 4 48.8 ▇▂▂▂▁
Let us start filtering the data, leading off with a definition of the minimum number of RNAs, minimum amount of rRNA, and maximum mitochondrial. In addition, let us print how much of each are observed before filtering. Before we can print/filter these attributes, we must use the PercentageFeatureSet() to get the numbers…
<- 200
min_num_rna <- 5
min_pct_ribo <- 20 max_pct_mito
Here is what skim looks like when given the entire dataset.
Genes observed:
skim(all_scd[["nFeature_RNA"]])
Name | all_scd[[“nFeature_RNA”]] |
Number of rows | 51594 |
Number of columns | 1 |
_______________________ | |
Column type frequency: | |
numeric | 1 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
---|---|---|---|---|---|---|---|---|---|---|
nFeature_RNA | 0 | 1 | 2078 | 1242 | 18 | 1124 | 1702 | 2798 | 8128 | ▇▆▂▁▁ |
Number counts:
skim(all_scd[["nCount_RNA"]])
Name | all_scd[[“nCount_RNA”]] |
Number of rows | 51594 |
Number of columns | 1 |
_______________________ | |
Column type frequency: | |
numeric | 1 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
---|---|---|---|---|---|---|---|---|---|---|
nCount_RNA | 0 | 1 | 5623 | 5052 | 500 | 2141 | 3659 | 7659 | 53836 | ▇▁▁▁▁ |
Percent ribosomal proteins:
skim(all_scd[["pct_mito"]])
Name | all_scd[[“pct_mito”]] |
Number of rows | 51594 |
Number of columns | 1 |
_______________________ | |
Column type frequency: | |
numeric | 1 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
---|---|---|---|---|---|---|---|---|---|---|
pct_mito | 0 | 1 | 2.26 | 3.32 | 0 | 1.28 | 1.88 | 2.69 | 98.86 | ▇▁▁▁▁ |
Percent mitochondrial RNA:
skim(all_scd[["pct_ribo"]])
Name | all_scd[[“pct_ribo”]] |
Number of rows | 51594 |
Number of columns | 1 |
_______________________ | |
Column type frequency: | |
numeric | 1 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
---|---|---|---|---|---|---|---|---|---|---|
pct_ribo | 0 | 1 | 11.26 | 9.15 | 0 | 5.19 | 7.46 | 13.84 | 51.42 | ▇▂▁▁▁ |
Number of clonotypes:
## Length and reads are for only those cells with clonotypes.
skim(all_scd[["reads"]])
Name | all_scd[[“reads”]] |
Number of rows | 51594 |
Number of columns | 1 |
_______________________ | |
Column type frequency: | |
numeric | 1 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
---|---|---|---|---|---|---|---|---|---|---|
reads | 39253 | 0.24 | 4425 | 4534 | 24 | 1639 | 3089 | 5505 | 53188 | ▇▁▁▁▁ |
And just for my curiosity, how many cells have some chain associated with them?
## How many cells have specific chains associated with them
sum(!is.na(all_scd$chain))
## [1] 12341
::kable(all_scd@misc[["sample_metadata"]]) knitr
sampleid | condition | dataroot | gexcells | gexmedianrpc | gexmediangpc | gextotalgenes | gexmedianupc | vdjtcells | vdjtvjspanning | mediantraupc | mediantrbpc | areads | breads | mappedgenomepct | mappedtxpct | vdjreads | vdjrpc | vdjmeanupc | batch | num_cells | mean_nFeature_RNA | mean_nCount_RNA | mean_clonotype_reads | mean_pct_mito | mean_pct_ribo | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
control | control | control | 01multi_combined/control | 15180 | 8398 | 1701 | 21761 | 3618 | 1971 | 1514 | 4 | 10 | 88882712 | 118449275 | 92.12 | 66.90 | 61009543 | 30954 | 15025 | undefined | 15180 | 1999 | 5115 | 7066 | 2.473 | 10.001 |
mock | mock | mock | 01multi_combined/mock | 13925 | 8295 | 1734 | 21536 | 3665 | 3090 | 2466 | 4 | 10 | 99158139 | 100331735 | 93.28 | 68.35 | 56928487 | 18423 | 10219 | undefined | 13925 | 2066 | 5729 | 3685 | 2.322 | 13.305 |
m | m | muscular | 01multi_combined/m | 13881 | 7736 | 1651 | 21977 | 3601 | 4217 | 3385 | 4 | 10 | 104613172 | 89191429 | 92.56 | 68.02 | 56660131 | 13436 | 8384 | undefined | 13881 | 2086 | 5539 | 4709 | 2.096 | 9.083 |
n | n | nasal | 01multi_combined/n | 8608 | 8336 | 1735 | 21183 | 3912 | 3063 | 2489 | 4 | 11 | 53224250 | 81908022 | 92.46 | 68.30 | 38838833 | 12680 | 8056 | undefined | 8608 | 2224 | 6488 | 3460 | 2.065 | 13.675 |
Seurat provides a violin plot for visualizing the distribution of some of these metrics easy. So lets look!
VlnPlot(all_scd, features="nFeature_RNA", pt.size=0)
VlnPlot(all_scd, features="pct_mito", pt.size=0)
VlnPlot(all_scd, features="pct_ribo", pt.size=0)
VlnPlot(all_scd, features="nCount_RNA", pt.size=0)
VlnPlot(all_scd, features="reads", pt.size=0)
## Warning: Removed 39253 rows containing non-finite values (`stat_ydensity()`).
## I am curious about the length of the clonotype sequences.
VlnPlot(all_scd, features="length", pt.size=0)
## Warning: Removed 39253 rows containing non-finite values (`stat_ydensity()`).
I wrote a little helper function which should make it easy to filter and observe what the data looks like before/after filtering.
<- filter_scd(all_scd, min_num_rna = 500) filtered
## All filters removed 12510 (24.247005%) cells, 10607 (32.854267%) genes, 55591307 (19.160320%) counts.
"pre_plots"]][["count_vs_genes"]] filtered[[
"post_plots"]][["count_vs_genes"]] filtered[[
"pre_plots"]][["count_vs_ribo"]] filtered[[
"post_plots"]][["count_vs_ribo"]] filtered[[
"pre_plots"]][["count_vs_mito"]] filtered[[
"post_plots"]][["count_vs_mito"]] filtered[[
"pre_plots"]][["ribo_vs_mito"]] filtered[[
"post_plots"]][["ribo_vs_mito"]] filtered[[
<- filtered[["scd"]] filt_scd
Now let us look at the various metrics following filtering as mean/sample.
::kable(filt_scd@misc[["sample_metadata"]]) knitr
sampleid | condition | dataroot | gexcells | gexmedianrpc | gexmediangpc | gextotalgenes | gexmedianupc | vdjtcells | vdjtvjspanning | mediantraupc | mediantrbpc | areads | breads | mappedgenomepct | mappedtxpct | vdjreads | vdjrpc | vdjmeanupc | batch | num_cells | mean_nFeature_RNA | mean_nCount_RNA | mean_clonotype_reads | mean_pct_mito | mean_pct_ribo | mean_filt_nfeature | mean_filt_ncount | mean_filt_pct_mito | mean_filt_pct_ribo | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
control | control | control | 01multi_combined/control | 15180 | 8398 | 1701 | 21761 | 3618 | 1971 | 1514 | 4 | 10 | 88882712 | 118449275 | 92.12 | 66.90 | 61009543 | 30954 | 15025 | undefined | 12016 | 1999 | 5115 | 7066 | 2.473 | 10.001 | 2107 | 5555 | 2.464 | 11.28 |
mock | mock | mock | 01multi_combined/mock | 13925 | 8295 | 1734 | 21536 | 3665 | 3090 | 2466 | 4 | 10 | 99158139 | 100331735 | 93.28 | 68.35 | 56928487 | 18423 | 10219 | undefined | 8920 | 2066 | 5729 | 3685 | 2.322 | 13.305 | 2143 | 6166 | 2.339 | 15.52 |
m | m | muscular | 01multi_combined/m | 13881 | 7736 | 1651 | 21977 | 3601 | 4217 | 3385 | 4 | 10 | 104613172 | 89191429 | 92.56 | 68.02 | 56660131 | 13436 | 8384 | undefined | 10999 | 2086 | 5539 | 4709 | 2.096 | 9.083 | 2088 | 5741 | 2.322 | 11.80 |
n | n | nasal | 01multi_combined/n | 8608 | 8336 | 1735 | 21183 | 3912 | 3063 | 2489 | 4 | 11 | 53224250 | 81908022 | 92.46 | 68.30 | 38838833 | 12680 | 8056 | undefined | 7149 | 2224 | 6488 | 3460 | 2.065 | 13.675 | 2282 | 6823 | 2.027 | 15.48 |
<- NormalizeData(object=all_scd) %>%
all_norm FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors() %>%
FindClusters() %>%
RunTSNE(check_duplicates=FALSE) %>%
RunUMAP(reduction="pca", dims=1:10)
## Centering and scaling data matrix
## PC_ 1
## Positive: Rac2, Coro1a, Arhgdib, Laptm5, Cd3e, Ms4a4b, Selplg, Cd52, Cd3g, Ptprc
## Cd3d, Lat, Gimap3, Lck, Ms4a6b, Thy1, Ltb, Ptprcap, Vps37b, Sept1
## Il7r, Skap1, Itgb7, Lsp1, Cd53, Cytip, Gramd3, Cd37, Cd2, Fyb
## Negative: Sparc, Mettl7a1, Ager, Cd63, Emp2, Cldn18, Selenbp1, Cystm1, Selenbp2, Alcam
## Gsn, Limch1, Aqp5, Clic3, Krt7, Cryab, Timp3, Npnt, Scnn1a, Igfbp6
## Krt18, Igfbp7, Myh14, Gprc5a, Scd2, Krt8, Fads3, Sec14l3, Mal2, Prnp
## PC_ 2
## Positive: Wfdc2, Atp1b1, Ctsh, Cldn3, Cbr2, Sftpd, Napsa, Sftpb, Cldn7, Cxcl15
## Chchd10, Slc34a2, Muc1, Sftpa1, Ppp1r14c, Chil1, Lamp3, Epcam, Ptprf, Krt18
## Irx1, Lgi3, Sftpc, Lyz2, Car8, Sdc1, Baiap2l1, Nkx2-1, Scd1, Krt8
## Negative: Sod3, Serping1, Prelp, Bgn, Olfml3, Plpp3, Loxl1, Spon1, Tcf21, Hsd11b1
## Rarres2, Plxdc2, Col1a2, Col1a1, Fn1, Pdgfra, Itga8, Cdo1, Ptgis, Dpep1
## Clec3b, Mfap4, Colec12, Gpx3, Lrp1, Pcolce2, Col3a1, C7, Atp1a2, Mxra8
## PC_ 3
## Positive: Rtkn2, Scnn1g, Spock2, Col4a3, Igfbp2, Pdpn, Scnn1b, Col4a4, Flrt3, Sema3e
## Cytl1, Lama3, Hopx, Lamb3, Itgb6, Ndnf, Tmem37, Clic5, Bdnf, Cyp2b10
## Crlf1, Ptpre, Sema3a, Tspan8, Krt7, Tacstd2, Cdkn2b, Cyp4b1, Abca5, Hck
## Negative: Cbr2, Cldn3, Mgst1, Nupr1, Ppp1r14c, Car8, Cxcl15, Sftpd, Chil1, Muc1
## Dram1, Sftpb, Slc34a2, Dynlrb2, Hp, Sftpa1, Lamp3, Cxcl17, Acot1, Lyz2
## Foxj1, Napsa, Scd1, Dmkn, Lbp, Fam183b, Cpm, Sftpc, Lgi3, Hdc
## PC_ 4
## Positive: Slc34a2, Lamp3, Napsa, Sftpb, Sftpa1, Sftpc, Cxcl15, Chil1, Lyz2, Dram1
## Lgi3, Sftpd, Hc, Scd1, Egfl6, S100g, Muc1, Pla2g1b, Etv5, Fabp5
## Ctsc, Lpcat1, Sfta2, Slc26a9, Mlc1, Chia1, Ppp1r14c, Tspan11, Fasn, Kcnj15
## Negative: Fam183b, Tmem212, Foxj1, Ccdc153, 1700007K13Rik, Cyp2s1, 1110017D15Rik, Gm19935, Arhgdig, Tctex1d4
## Spaca9, BC051019, 1700094D03Rik, Cfap126, 1700001C02Rik, 1700016K19Rik, Rsph1, Gm867, Ccdc113, Nme5
## Dnali1, Tekt4, AU040972, Pifo, Odf3b, Sntn, Ak7, 2410004P03Rik, Nme9, Mns1
## PC_ 5
## Positive: Lgals1, Coro1a, Rac2, Laptm5, Arhgdib, Cd3e, Cd52, Selplg, Cd3g, Ms4a4b
## Ptprc, Lsp1, Ptprcap, Lat, Thy1, Cytip, Sept1, Cd3d, Lck, Ltb
## Ms4a6b, Itgb7, Ifi27l2a, Gimap3, Skap1, Cd37, Icos, Il2rb, Cd2, Cd53
## Negative: Ly6a, Tm4sf1, Ly6c1, Sema3c, Kdr, Tmem252, Car4, H2-Ab1, Plaur, Tmcc2
## Id1, Efnb2, Cd74, Gja4, Lyve1, Tbx3, Sox17, Emcn, Sema7a, Plk2
## H2-Eb1, Scarb1, Nes, Hspb1, Cyp4b1, Prx, Ccnd1, Serpina3i, Id3, H2-Aa
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 51594
## Number of edges: 1601833
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9148
## Number of communities: 29
## Elapsed time: 17 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 09:44:01 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:44:01 Read 51594 rows and found 10 numeric columns
## 09:44:01 Using Annoy for neighbor search, n_neighbors = 30
## 09:44:01 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:44:10 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a3d3288c2
## 09:44:10 Searching Annoy index using 1 thread, search_k = 3000
## 09:44:39 Annoy recall = 100%
## 09:44:42 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:44:48 Initializing from normalized Laplacian + noise (using irlba)
## 09:44:59 Commencing optimization for 200 epochs, with 2133286 positive edges
## 09:45:26 Optimization finished
DimPlot(object=all_norm, reduction="tsne")
<- DimPlot(all_norm, reduction="umap", group.by="orig.ident", label=TRUE)
plotted plotted
<- NormalizeData(object=filt_scd) %>%
fnorm_scd FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors() %>%
FindClusters() %>%
RunTSNE() %>%
RunUMAP(reduction="pca", dims=1:10)
## Centering and scaling data matrix
## PC_ 1
## Positive: Sparc, Emp2, Cd63, Igfbp7, Mettl7a1, Ager, Gsn, Cldn18, Timp3, Limch1
## Cst3, Gpx3, Cystm1, Selenbp1, Cryab, Bgn, Npnt, Selenbp2, Bcam, Pmp22
## Clic3, Apoe, Hspb1, Krt7, Ces1d, Aqp5, Anxa3, Ifitm3, Alcam, Fhl1
## Negative: Rac2, Coro1a, Arhgdib, Laptm5, Cd3e, Ms4a4b, Selplg, Cd52, Cd3g, Ptprc
## Cd3d, Lat, Lck, Ms4a6b, Thy1, Ltb, Ptprcap, Sept1, Il7r, Itgb7
## Lsp1, Gramd3, Fyb, Cd2, Cd28, Icos, Cd69, Bin2, Il2rb, AW112010
## PC_ 2
## Positive: Bgn, Serping1, Sod3, Plpp3, Prelp, Rarres2, Col1a2, Sparcl1, Hsd11b1, Vim
## Olfml3, Col1a1, Loxl1, Tcf21, Clec3b, Ptgis, Mxra8, Spon1, Plxdc2, Fxyd1
## Gpx3, Fhl1, Ltbp4, Col3a1, Lrp1, Ppp1r14a, Cdo1, Pcolce2, Mfap4, Colec12
## Negative: Sftpd, Napsa, Wfdc2, Sftpb, Cxcl15, Slc34a2, Cbr2, Sftpa1, Cldn3, Atp1b1
## Chil1, Lamp3, Muc1, Ppp1r14c, Sftpc, Lgi3, Lyz2, Ctsh, Dram1, Scd1
## Hc, Car8, Egfl6, Lpcat1, S100g, Irx1, Pi4k2b, Abca3, Pla2g1b, Ank3
## PC_ 3
## Positive: Rtkn2, Scnn1g, Hopx, Cyp4b1, Igfbp2, Clic5, Spock2, Col4a3, Sema3e, Pdpn
## Cyp2b10, Scnn1b, Tspan8, Col4a4, Lama3, Flrt3, Cytl1, Krt7, Tacstd2, Sec14l3
## Cdkn2b, Lamb3, Itgb6, Krt19, Mab21l4, Pakap.1, Tmem37, Scnn1a, Lmo7, Akap5
## Negative: Mgst1, Col6a1, Nupr1, Loxl1, Col3a1, Col1a1, Sod3, Col1a2, C3, Dram1
## Ces1d, Col6a2, Apoe, Slc34a2, Mmp2, Plxdc2, Cxcl15, Pdgfra, Dpep1, Chil1
## Lamp3, Sftpb, Spon1, Prelp, Cd302, Sftpd, Olfml3, Lyz2, Atp1a2, Sftpa1
## PC_ 4
## Positive: Cldn5, Plvap, Hpgd, H2-Ab1, Cd93, Cd74, Tmem252, H2-Aa, H2-Eb1, Sema3c
## Pcdh1, Mcam, Gja4, Efnb2, Sox17, Dll4, Lyve1, Stmn2, Tm4sf1, Hspa1a
## Kdr, Ly6c1, Plk2, Hilpda, Tmcc2, Emcn, Hspa1b, Scn7a, Ly6a, Ifitm3
## Negative: S100a6, Lgals1, Birc5, Lgals3, Mki67, Anxa1, Pclaf, Prdx6, Tacstd2, Top2a
## Pdpn, Cyp2b10, Sema3e, Ccna2, Krt19, Scnn1g, Coro1a, Ube2c, Spock2, Rac2
## Scnn1b, S100a11, Rtkn2, Cdk1, Mab21l4, Arhgdig, Col4a4, Ccnb2, Rrm2, Col4a3
## PC_ 5
## Positive: Fam183b, Foxj1, Tmem212, 1700007K13Rik, 1110017D15Rik, Ccdc153, Gm19935, Cfap126, Spaca9, Tctex1d4
## BC051019, 1700094D03Rik, Odf3b, 1700001C02Rik, 1700016K19Rik, Tekt4, Dnali1, Rsph1, AU040972, Cyp2s1
## Ak7, Ccdc113, Sntn, Gm867, Nme5, Nme9, 2410004P03Rik, Dynlrb2, Pifo, Lrrc10b
## Negative: Col4a4, Aqp5, Scnn1g, Ager, Rtkn2, Ndnf, Col4a3, Slc39a8, Gprc5a, Spock2
## Tmem37, Pdpn, Scnn1b, Igfbp6, Igfbp2, Flrt3, Crlf1, Scd2, Lamb3, Fads3
## Cldn18, Rgcc, Itgb6, Mal2, Tmod1, Fam189a2, Rnase4, Hs2st1, Abca5, Lamc2
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 1210541
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9102
## Number of communities: 25
## Elapsed time: 10 seconds
## 09:51:06 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:51:06 Read 39084 rows and found 10 numeric columns
## 09:51:06 Using Annoy for neighbor search, n_neighbors = 30
## 09:51:06 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:51:13 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a6482ce89
## 09:51:13 Searching Annoy index using 1 thread, search_k = 3000
## 09:51:35 Annoy recall = 100%
## 09:51:38 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:51:43 Initializing from normalized Laplacian + noise (using irlba)
## 09:51:52 Commencing optimization for 200 epochs, with 1600638 positive edges
## 09:52:15 Optimization finished
DimPlot(object=fnorm_scd, reduction="tsne")
<- DimPlot(fnorm_scd, reduction="umap", group.by="orig.ident", label=TRUE)
plotted plotted
<- JackStraw(fnorm_scd, num.replicate=10)
fnorm_scd <- ScoreJackStraw(fnorm_scd)
fnorm_scd JackStrawPlot(fnorm_scd)
## Warning: Removed 7000 rows containing missing values (`geom_point()`).
ElbowPlot(fnorm_scd)
## So I am thinking maybe 4-10?
<- 6
wanted_dims
<- FindNeighbors(fnorm_scd, dims=1:wanted_dims) %>%
fnorm_scd FindClusters(resolution=0.5) %>%
StashIdent(save.name="res0p5_clusters")
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 1144964
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9257
## Number of communities: 19
## Elapsed time: 12 seconds
## With Seurat 3.X, stashing identity classes can be accomplished with the following:
## .[["res0p5_clusters"]] <- Idents(object = .)
RunUMAP(fnorm_scd, dims=1:9)
## 09:55:20 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:55:20 Read 39084 rows and found 9 numeric columns
## 09:55:20 Using Annoy for neighbor search, n_neighbors = 30
## 09:55:20 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:55:27 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a5d935ddb
## 09:55:27 Searching Annoy index using 1 thread, search_k = 3000
## 09:55:49 Annoy recall = 100%
## 09:55:52 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:55:58 Initializing from normalized Laplacian + noise (using irlba)
## 09:56:07 Commencing optimization for 200 epochs, with 1594748 positive edges
## 09:56:28 Optimization finished
## An object of class Seurat
## 21678 features across 39084 samples within 1 assay
## Active assay: RNA (21678 features, 2000 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
DimPlot(fnorm_scd, label=TRUE)
<- FindClusters(fnorm_scd, resolution=0.1) %>%
fnorm_scd FindNeighbors(k.param=6) %>%
StashIdent(save.name="res0p1_clusters")
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 1144964
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9712
## Number of communities: 7
## Elapsed time: 8 seconds
## Computing nearest neighbor graph
## Computing SNN
## With Seurat 3.X, stashing identity classes can be accomplished with the following:
## .[["res0p1_clusters"]] <- Idents(object = .)
RunUMAP(fnorm_scd, dims=1:9)
## 09:56:52 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:56:52 Read 39084 rows and found 9 numeric columns
## 09:56:52 Using Annoy for neighbor search, n_neighbors = 30
## 09:56:52 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:56:59 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a53a136ee
## 09:56:59 Searching Annoy index using 1 thread, search_k = 3000
## 09:57:21 Annoy recall = 100%
## 09:57:24 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:57:30 Initializing from normalized Laplacian + noise (using irlba)
## 09:57:39 Commencing optimization for 200 epochs, with 1594748 positive edges
## 09:58:01 Optimization finished
## An object of class Seurat
## 21678 features across 39084 samples within 1 assay
## Active assay: RNA (21678 features, 2000 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
DimPlot(fnorm_scd, label=TRUE)
<- FindClusters(fnorm_scd, resolution=0.2) %>%
fnorm_scd FindNeighbors(k.param=10) %>%
StashIdent(save.name="res0p2_clusters")
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 507533
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9586
## Number of communities: 14
## Elapsed time: 6 seconds
## Computing nearest neighbor graph
## Computing SNN
## With Seurat 3.X, stashing identity classes can be accomplished with the following:
## .[["res0p2_clusters"]] <- Idents(object = .)
RunUMAP(fnorm_scd, dims=1:10)
## 09:58:26 UMAP embedding parameters a = 0.9922 b = 1.112
## 09:58:26 Read 39084 rows and found 10 numeric columns
## 09:58:26 Using Annoy for neighbor search, n_neighbors = 30
## 09:58:26 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 09:58:33 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a79ca4c81
## 09:58:33 Searching Annoy index using 1 thread, search_k = 3000
## 09:58:54 Annoy recall = 100%
## 09:58:57 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 09:59:03 Initializing from normalized Laplacian + noise (using irlba)
## 09:59:11 Commencing optimization for 200 epochs, with 1600638 positive edges
## 09:59:33 Optimization finished
## An object of class Seurat
## 21678 features across 39084 samples within 1 assay
## Active assay: RNA (21678 features, 2000 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
DimPlot(fnorm_scd, label=TRUE)
Add into the metadata a concatenation of the sample ID and the cluster ID
<- fnorm_scd[["orig.ident"]][["orig.ident"]]
identity_vector <- as.character(fnorm_scd[["res0p1_clusters"]][["res0p1_clusters"]])
cluster_vector <- paste0(identity_vector, "_", cluster_vector)
concatenated_vector "cluster1_sample"]] <- concatenated_vector
fnorm_scd[[summary(as.factor(fnorm_scd@meta.data[["cluster1_sample"]]))
## control_0 control_1 control_2 control_3 control_4 control_5 control_6 m_0
## 4790 1847 2427 1513 928 462 49 2869
## m_1 m_2 m_3 m_4 m_5 m_6 mock_0 mock_1
## 3790 1955 1009 730 583 63 2765 2877
## mock_2 mock_3 mock_4 mock_5 mock_6 n_0 n_1 n_2
## 1765 760 459 269 25 1000 2893 1595
## n_3 n_4 n_5 n_6
## 657 669 277 58
<- as.character(fnorm_scd[["res0p2_clusters"]][["res0p2_clusters"]])
cluster_vector2 <- paste0(identity_vector, "_", cluster_vector2)
concatenated_vector2 "cluster2_sample"]] <- concatenated_vector2
fnorm_scd[[summary(as.factor(fnorm_scd@meta.data[["cluster2_sample"]]))
## control_0 control_1 control_10 control_11 control_12 control_13 control_2
## 3504 1056 211 429 130 49 768
## control_3 control_4 control_5 control_6 control_7 control_8 control_9
## 1516 1362 895 905 762 319 110
## m_0 m_1 m_10 m_11 m_12 m_13 m_2
## 1818 2437 354 170 66 63 1300
## m_3 m_4 m_5 m_6 m_7 m_8 m_9
## 1169 864 773 699 560 287 439
## mock_0 mock_1 mock_10 mock_11 mock_12 mock_13 mock_2
## 1828 1052 274 155 91 26 1797
## mock_3 mock_4 mock_5 mock_6 mock_7 mock_8 mock_9
## 1081 668 659 436 520 190 143
## n_0 n_1 n_10 n_11 n_12 n_13 n_2
## 571 1556 72 58 38 58 1322
## n_3 n_4 n_5 n_6 n_7 n_8 n_9
## 709 437 862 640 371 227 228
I am not yet certain of how Seurat handles (non)normalized data for the various FindMarkers functions. Thus, I am adding the clusters from the dimension reductions to the non-normalized data here.
"res0p1_clusters"]] <- fnorm_scd[["res0p1_clusters"]]
filt_scd[["res0p2_clusters"]] <- fnorm_scd[["res0p2_clusters"]]
filt_scd[["cluster1_sample"]] <- fnorm_scd[["cluster1_sample"]]
filt_scd[["cluster2_sample"]] <- fnorm_scd[["cluster2_sample"]] filt_scd[[
If we decide to limit the comparisons to the relatively small subset of genes with significant variance, these will be the winners; so we can decide now if that is a goal. Given that this dataset is not super-huge, I am not included to do it; but it is nice to see the set of highly-variable genes.
<- FindVariableFeatures(fnorm_scd)
var <- head(VariableFeatures(var), 30)
most_var <- VariableFeaturePlot(var)
variable_plot <- LabelPoints(plot=variable_plot, points=most_var, repel=TRUE) variable_plot
## When using repel, set xnudge and ynudge to 0 for optimal results
variable_plot
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: ggrepel: 15 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
I assume that Dr. Park would like to play with these various results in a xlsx file, so I will dump them all to a list and use that to define the sheets of this putative output file.
<- list()
de_lst <- 0.6 fc_threshold
The Clonotype containing cells are pink, the rest are blue in the final plot in the following block.
<- filt_scd
vdj_scd Idents(vdj_scd) <- vdj_scd[["tcell"]]
<- NormalizeData(vdj_scd) %>%
vdj_scd_viz FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors() %>%
FindClusters() %>%
RunTSNE() %>%
RunUMAP(reduction="pca", dims=1:10)
## Centering and scaling data matrix
## PC_ 1
## Positive: Sparc, Emp2, Cd63, Igfbp7, Mettl7a1, Ager, Gsn, Cldn18, Timp3, Limch1
## Cst3, Gpx3, Cystm1, Selenbp1, Cryab, Bgn, Npnt, Selenbp2, Bcam, Pmp22
## Clic3, Apoe, Hspb1, Krt7, Ces1d, Aqp5, Anxa3, Ifitm3, Alcam, Fhl1
## Negative: Rac2, Coro1a, Arhgdib, Laptm5, Cd3e, Ms4a4b, Selplg, Cd52, Cd3g, Ptprc
## Cd3d, Lat, Lck, Ms4a6b, Thy1, Ltb, Ptprcap, Sept1, Il7r, Itgb7
## Lsp1, Gramd3, Fyb, Cd2, Cd28, Icos, Cd69, Bin2, Il2rb, AW112010
## PC_ 2
## Positive: Bgn, Serping1, Sod3, Plpp3, Prelp, Rarres2, Col1a2, Sparcl1, Hsd11b1, Vim
## Olfml3, Col1a1, Loxl1, Tcf21, Clec3b, Ptgis, Mxra8, Spon1, Plxdc2, Fxyd1
## Gpx3, Fhl1, Ltbp4, Col3a1, Lrp1, Ppp1r14a, Cdo1, Pcolce2, Mfap4, Colec12
## Negative: Sftpd, Napsa, Wfdc2, Sftpb, Cxcl15, Slc34a2, Cbr2, Sftpa1, Cldn3, Atp1b1
## Chil1, Lamp3, Muc1, Ppp1r14c, Sftpc, Lgi3, Lyz2, Ctsh, Dram1, Scd1
## Hc, Car8, Egfl6, Lpcat1, S100g, Irx1, Pi4k2b, Abca3, Pla2g1b, Ank3
## PC_ 3
## Positive: Rtkn2, Scnn1g, Hopx, Cyp4b1, Igfbp2, Clic5, Spock2, Col4a3, Sema3e, Pdpn
## Cyp2b10, Scnn1b, Tspan8, Col4a4, Lama3, Flrt3, Cytl1, Krt7, Tacstd2, Sec14l3
## Cdkn2b, Lamb3, Itgb6, Krt19, Mab21l4, Pakap.1, Tmem37, Scnn1a, Lmo7, Akap5
## Negative: Mgst1, Col6a1, Nupr1, Loxl1, Col3a1, Col1a1, Sod3, Col1a2, C3, Dram1
## Ces1d, Col6a2, Apoe, Slc34a2, Mmp2, Plxdc2, Cxcl15, Pdgfra, Dpep1, Chil1
## Lamp3, Sftpb, Spon1, Prelp, Cd302, Sftpd, Olfml3, Lyz2, Atp1a2, Sftpa1
## PC_ 4
## Positive: Cldn5, Plvap, Hpgd, H2-Ab1, Cd93, Cd74, Tmem252, H2-Aa, H2-Eb1, Sema3c
## Pcdh1, Mcam, Gja4, Efnb2, Sox17, Dll4, Lyve1, Stmn2, Tm4sf1, Hspa1a
## Kdr, Ly6c1, Plk2, Hilpda, Tmcc2, Emcn, Hspa1b, Scn7a, Ly6a, Ifitm3
## Negative: S100a6, Lgals1, Birc5, Lgals3, Mki67, Anxa1, Pclaf, Prdx6, Tacstd2, Top2a
## Pdpn, Cyp2b10, Sema3e, Ccna2, Krt19, Scnn1g, Coro1a, Ube2c, Spock2, Rac2
## Scnn1b, S100a11, Rtkn2, Cdk1, Mab21l4, Arhgdig, Col4a4, Ccnb2, Rrm2, Col4a3
## PC_ 5
## Positive: Fam183b, Foxj1, Tmem212, 1700007K13Rik, 1110017D15Rik, Ccdc153, Gm19935, Cfap126, Spaca9, Tctex1d4
## BC051019, 1700094D03Rik, Odf3b, 1700001C02Rik, 1700016K19Rik, Tekt4, Dnali1, Rsph1, AU040972, Cyp2s1
## Ak7, Ccdc113, Sntn, Gm867, Nme5, Nme9, 2410004P03Rik, Dynlrb2, Pifo, Lrrc10b
## Negative: Col4a4, Aqp5, Scnn1g, Ager, Rtkn2, Ndnf, Col4a3, Slc39a8, Gprc5a, Spock2
## Tmem37, Pdpn, Scnn1b, Igfbp6, Igfbp2, Flrt3, Crlf1, Scd2, Lamb3, Fads3
## Cldn18, Rgcc, Itgb6, Mal2, Tmod1, Fam189a2, Rnase4, Hs2st1, Abca5, Lamc2
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 1210541
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9102
## Number of communities: 25
## Elapsed time: 10 seconds
## 10:05:15 UMAP embedding parameters a = 0.9922 b = 1.112
## 10:05:15 Read 39084 rows and found 10 numeric columns
## 10:05:15 Using Annoy for neighbor search, n_neighbors = 30
## 10:05:15 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 10:05:22 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a2033519b
## 10:05:22 Searching Annoy index using 1 thread, search_k = 3000
## 10:05:43 Annoy recall = 100%
## 10:05:46 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 10:05:52 Initializing from normalized Laplacian + noise (using irlba)
## 10:06:00 Commencing optimization for 200 epochs, with 1600638 positive edges
## 10:06:23 Optimization finished
DimPlot(vdj_scd_viz, reduction="tsne")
Idents(vdj_scd_viz) <- vdj_scd_viz[["tcell"]]
<- DimPlot(vdj_scd_viz, reduction="umap", group.by="ident", label=TRUE)
plotted plotted
<- FindAllMarkers(vdj_scd, logfc.threshold=fc_threshold) vdj_markers
## Calculating cluster Yes
## Calculating cluster No
head(vdj_markers)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## Ms4a4b 0 2.250 0.889 0.053 0 Yes Ms4a4b
## Ccl5 0 2.199 0.274 0.034 0 Yes Ccl5
## Rac2 0 2.054 0.932 0.113 0 Yes Rac2
## Coro1a 0 2.042 0.932 0.121 0 Yes Coro1a
## Cd3e 0 1.884 0.928 0.088 0 Yes Cd3e
## Arhgdib 0 1.837 0.918 0.162 0 Yes Arhgdib
<- as.data.frame(vdj_markers)
combined rownames(combined) <- toupper(rownames(combined))
<- merge(combined, brief, by.x="row.names", by.y="hgnc_symbol",
annotated_markers all.x=TRUE)
"vdj_vs_all"]] <- annotated_markers
de_lst[[head(annotated_markers)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 ABCA3 0 -0.8406 0.115 0.376 0 Yes Abca3
## 2 ABCA3.1 0 0.8406 0.376 0.115 0 No Abca3
## 3 ACE 0 -1.4942 0.130 0.648 0 Yes Ace
## 4 ACE.1 0 1.4942 0.648 0.130 0 No Ace
## 5 ACER2 0 -0.6631 0.092 0.406 0 Yes Acer2
## 6 ACER2.1 0 0.6631 0.406 0.092 0 No Acer2
## description
## 1 ATP binding cassette subfamily A member 3 [Source:HGNC Symbol;Acc:HGNC:33]
## 2 <NA>
## 3 angiotensin I converting enzyme [Source:HGNC Symbol;Acc:HGNC:2707]
## 4 <NA>
## 5 alkaline ceramidase 2 [Source:HGNC Symbol;Acc:HGNC:23675]
## 6 <NA>
## ensembl_gene_id
## 1 ENSG00000167972
## 2 <NA>
## 3 ENSG00000159640
## 4 <NA>
## 5 ENSG00000177076
## 6 <NA>
Question: Is it smart enough to use the raw data if I give FindAllMarkers the normalized data? For the moment I do not think I will risk it.
<- FindAllMarkers(filt_scd, logfc.threshold=fc_threshold) combined_markers
## Calculating cluster control
## Calculating cluster m
## Calculating cluster mock
## Calculating cluster n
head(combined_markers)
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## Tmem252 0 1.0119 0.499 0.315 0 control Tmem252
## Plvap 0 0.8076 0.605 0.415 0 control Plvap
## Ms4a6b 0 -0.6020 0.112 0.335 0 control Ms4a6b
## Phf11d 0 -0.6152 0.124 0.305 0 control Phf11d
## Samhd1 0 -0.6259 0.364 0.578 0 control Samhd1
## H2-Q7 0 -0.6674 0.468 0.682 0 control H2-Q7
<- as.data.frame(combined_markers)
combined rownames(combined) <- toupper(rownames(combined))
<- merge(combined, brief, by.x="row.names", by.y="hgnc_symbol",
annotated_markers all.x=TRUE)
"samples_vs_all"]] <- annotated_markers de_lst[[
Since I am not using the fnorm_scd data structure, I will need to pull the cluster information from the normalized copy…
<- filt_scd
cluster_scd Idents(cluster_scd) <- cluster_scd[["res0p2_clusters"]]
<- FindAllMarkers(cluster_scd, only.pos=TRUE, logfc.threshold=fc_threshold) cluster_markers
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
## Calculating cluster 12
## Calculating cluster 13
<- as.data.frame(cluster_markers)
cluster_genes rownames(cluster_genes) <- toupper(rownames(cluster_genes))
<- merge(cluster_genes, brief, by.x="row.names", by.y="hgnc_symbol",
annotated_clusters all.x=TRUE)
head(annotated_clusters)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj cluster
## 1 0610010K14RIK 3.960e-173 0.7199 0.611 0.246 8.585e-169 9
## 2 1110004F10RIK 0.000e+00 0.6234 0.666 0.329 0.000e+00 6
## 3 1110008P14RIK 0.000e+00 0.7656 0.707 0.286 0.000e+00 5
## 4 1110008P14RIK1 0.000e+00 0.7679 0.701 0.292 0.000e+00 6
## 5 1110017D15RIK 0.000e+00 2.9414 0.939 0.002 0.000e+00 13
## 6 1110032A03RIK 7.044e-193 1.1535 0.719 0.102 1.527e-188 13
## gene description ensembl_gene_id
## 1 0610010K14Rik <NA> <NA>
## 2 1110004F10Rik <NA> <NA>
## 3 1110008P14Rik <NA> <NA>
## 4 1110008P14Rik <NA> <NA>
## 5 1110017D15Rik <NA> <NA>
## 6 1110032A03Rik <NA> <NA>
"samplecluster_vs_all"]] <- annotated_clusters
de_lst[[
%>%
annotated_clusters group_by(cluster) %>%
::top_n(n=10, wt=avg_log2FC) %>%
dplyras.data.frame()
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 ACTA2 2.625e-163 2.0222 0.140 0.030 5.691e-159 7 Acta2
## 2 AGER1 0.000e+00 3.8490 1.000 0.389 0.000e+00 6 Ager
## 3 AQP51 0.000e+00 3.0677 0.991 0.234 0.000e+00 6 Aqp5
## 4 AU040972 0.000e+00 4.8516 0.898 0.001 0.000e+00 13 AU040972
## 5 BTG1 0.000e+00 0.9181 0.943 0.675 0.000e+00 1 Btg1
## 6 C32 1.187e-145 3.2906 0.514 0.240 2.574e-141 8 C3
## 7 CBR2 0.000e+00 3.2504 0.989 0.144 0.000e+00 3 Cbr2
## 8 CBR22 5.288e-179 3.7216 0.980 0.237 1.146e-174 13 Cbr2
## 9 CCL11 0.000e+00 2.9691 0.221 0.004 0.000e+00 8 Ccl11
## 10 CCL21A 0.000e+00 2.9437 0.523 0.003 0.000e+00 12 Ccl21a
## 11 CCL41 0.000e+00 3.0918 0.344 0.025 0.000e+00 11 Ccl4
## 12 CCL5 0.000e+00 2.7185 0.433 0.059 0.000e+00 2 Ccl5
## 13 CCL52 3.358e-58 1.0894 0.272 0.105 7.280e-54 10 Ccl5
## 14 CCL6 0.000e+00 3.9990 0.729 0.005 0.000e+00 11 Ccl6
## 15 CCR7 0.000e+00 1.7068 0.661 0.077 0.000e+00 1 Ccr7
## 16 CD24A1 8.143e-227 3.7231 0.985 0.217 1.765e-222 13 Cd24a
## 17 CD3E1 0.000e+00 1.3816 0.923 0.261 0.000e+00 2 Cd3e
## 18 CD3E3 4.396e-199 0.7934 0.832 0.338 9.530e-195 10 Cd3e
## 19 CD3G 0.000e+00 1.4120 0.886 0.213 0.000e+00 2 Cd3g
## 20 CD3G2 5.893e-190 0.7413 0.754 0.292 1.277e-185 10 Cd3g
## 21 CD52 0.000e+00 1.3965 0.891 0.239 0.000e+00 2 Cd52
## 22 CDKN1C1 0.000e+00 1.9037 0.888 0.140 0.000e+00 7 Cdkn1c
## 23 CHIL1 0.000e+00 3.0752 0.958 0.114 0.000e+00 3 Chil1
## 24 CHIL11 0.000e+00 2.3075 0.953 0.145 0.000e+00 5 Chil1
## 25 CHIL3 0.000e+00 2.8178 0.586 0.001 0.000e+00 11 Chil3
## 26 CLDN5 0.000e+00 0.9772 0.914 0.385 0.000e+00 0 Cldn5
## 27 CLEC14A 0.000e+00 0.9318 0.806 0.316 0.000e+00 0 Clec14a
## 28 CLIC31 0.000e+00 3.0069 0.996 0.267 0.000e+00 6 Clic3
## 29 CLU3 1.884e-07 1.6426 0.363 0.276 4.084e-03 12 Clu
## 30 COL1A12 0.000e+00 2.9670 0.737 0.170 0.000e+00 8 Col1a1
## 31 COL1A22 0.000e+00 2.8036 0.693 0.197 0.000e+00 8 Col1a2
## 32 COL3A11 0.000e+00 3.4811 0.859 0.157 0.000e+00 8 Col3a1
## 33 CORO1A1 0.000e+00 3.0619 0.788 0.363 0.000e+00 9 Coro1a
## 34 COX4I22 0.000e+00 2.8431 0.979 0.144 0.000e+00 7 Cox4i2
## 35 CTLA2A 0.000e+00 0.9380 0.896 0.491 0.000e+00 0 Ctla2a
## 36 CTSS 0.000e+00 2.6584 0.881 0.091 0.000e+00 11 Ctss
## 37 CXCL15 0.000e+00 3.4288 0.995 0.145 0.000e+00 3 Cxcl15
## 38 CXCL151 0.000e+00 2.3410 0.983 0.176 0.000e+00 5 Cxcl15
## 39 CXCL2 0.000e+00 3.0352 0.539 0.048 0.000e+00 11 Cxcl2
## 40 CYBB 0.000e+00 2.5947 0.823 0.007 0.000e+00 11 Cybb
## 41 CYP2B101 0.000e+00 2.9079 0.971 0.127 0.000e+00 6 Cyp2b10
## 42 CYP2F21 0.000e+00 4.6976 0.837 0.063 0.000e+00 13 Cyp2f2
## 43 CYTL12 2.662e-35 1.6723 0.372 0.151 5.770e-31 12 Cytl1
## 44 DCN 0.000e+00 4.4693 0.745 0.013 0.000e+00 8 Dcn
## 45 DUSP10 0.000e+00 0.8373 0.563 0.156 0.000e+00 1 Dusp10
## 46 EGFL7 0.000e+00 1.0153 0.957 0.421 0.000e+00 0 Egfl7
## 47 EMP21 0.000e+00 2.9943 1.000 0.521 0.000e+00 6 Emp2
## 48 FAM183B 0.000e+00 3.8059 0.985 0.001 0.000e+00 13 Fam183b
## 49 FCER1G 0.000e+00 2.8584 0.898 0.029 0.000e+00 11 Fcer1g
## 50 FHL1 0.000e+00 1.9510 0.945 0.313 0.000e+00 4 Fhl1
## 51 FN1 0.000e+00 1.9900 0.906 0.208 0.000e+00 4 Fn1
## 52 GPIHBP1 0.000e+00 1.0345 0.904 0.380 0.000e+00 0 Gpihbp1
## 53 GPX3 0.000e+00 2.2317 0.996 0.429 0.000e+00 4 Gpx3
## 54 GRAMD3 0.000e+00 0.9715 0.699 0.194 0.000e+00 1 Gramd3
## 55 GSN 0.000e+00 2.3819 0.991 0.431 0.000e+00 4 Gsn
## 56 GUCY1A12 0.000e+00 2.6652 0.977 0.134 0.000e+00 7 Gucy1a1
## 57 GUCY1B11 0.000e+00 2.2422 0.946 0.118 0.000e+00 7 Gucy1b1
## 58 GYG 0.000e+00 1.9484 0.945 0.319 0.000e+00 4 Gyg
## 59 GZMA 0.000e+00 1.3923 0.140 0.016 0.000e+00 2 Gzma
## 60 H2AFZ1 0.000e+00 3.3186 0.997 0.764 0.000e+00 9 H2afz
## 61 HMGB2 0.000e+00 4.0079 0.991 0.387 0.000e+00 9 Hmgb2
## 62 HOPX1 0.000e+00 3.5571 1.000 0.524 0.000e+00 6 Hopx
## 63 HP2 5.651e-158 4.0197 0.898 0.190 1.225e-153 13 Hp
## 64 HPGD 0.000e+00 0.9877 0.809 0.348 0.000e+00 0 Hpgd
## 65 IGFBP21 0.000e+00 2.8773 0.844 0.088 0.000e+00 6 Igfbp2
## 66 IGFBP41 5.928e-89 2.0522 0.889 0.501 1.285e-84 12 Igfbp4
## 67 IGFBP5 0.000e+00 3.4222 0.624 0.024 0.000e+00 8 Igfbp5
## 68 IL1B 0.000e+00 2.7539 0.518 0.005 0.000e+00 11 Il1b
## 69 IL7R 0.000e+00 1.0182 0.744 0.182 0.000e+00 1 Il7r
## 70 IL7R2 1.028e-203 0.8539 0.718 0.259 2.230e-199 10 Il7r
## 71 INMT 0.000e+00 2.3950 0.975 0.322 0.000e+00 4 Inmt
## 72 KRT191 0.000e+00 2.8640 0.978 0.187 0.000e+00 6 Krt19
## 73 KRT71 0.000e+00 3.2883 0.999 0.276 0.000e+00 6 Krt7
## 74 LAPTM5 0.000e+00 1.3035 0.947 0.264 0.000e+00 2 Laptm5
## 75 LCN21 0.000e+00 2.3795 0.839 0.182 0.000e+00 5 Lcn2
## 76 LEF1 0.000e+00 0.8895 0.546 0.110 0.000e+00 1 Lef1
## 77 LGALS11 0.000e+00 3.5272 0.882 0.458 0.000e+00 9 Lgals1
## 78 LGALS32 0.000e+00 2.5184 0.899 0.304 0.000e+00 11 Lgals3
## 79 LTB 0.000e+00 1.3167 0.785 0.212 0.000e+00 2 Ltb
## 80 LTB2 6.739e-220 0.8539 0.756 0.277 1.461e-215 10 Ltb
## 81 LYZ1 0.000e+00 3.4318 0.576 0.098 0.000e+00 3 Lyz1
## 82 LYZ11 0.000e+00 2.4354 0.592 0.113 0.000e+00 5 Lyz1
## 83 LYZ2 0.000e+00 3.5296 0.962 0.199 0.000e+00 3 Lyz2
## 84 LYZ21 0.000e+00 2.4025 0.961 0.227 0.000e+00 5 Lyz2
## 85 MFAP4 0.000e+00 2.6586 0.985 0.323 0.000e+00 4 Mfap4
## 86 MFGE82 0.000e+00 2.0515 0.984 0.397 0.000e+00 7 Mfge8
## 87 MGP2 0.000e+00 3.2015 0.951 0.283 0.000e+00 8 Mgp
## 88 MGP3 6.435e-13 3.7647 0.418 0.300 1.395e-08 12 Mgp
## 89 MKI67 0.000e+00 3.0267 0.877 0.012 0.000e+00 9 Mki67
## 90 MMRN1 0.000e+00 2.3739 0.520 0.003 0.000e+00 12 Mmrn1
## 91 MS4A4B 0.000e+00 0.8832 0.831 0.217 0.000e+00 1 Ms4a4b
## 92 MS4A4B1 0.000e+00 1.5024 0.834 0.233 0.000e+00 2 Ms4a4b
## 93 MS4A4B3 1.115e-169 0.7910 0.757 0.302 2.418e-165 10 Ms4a4b
## 94 NDUFA4L21 0.000e+00 1.9368 0.896 0.085 0.000e+00 7 Ndufa4l2
## 95 NKG7 0.000e+00 1.5423 0.550 0.091 0.000e+00 2 Nkg7
## 96 PCLAF 0.000e+00 2.9114 0.840 0.005 0.000e+00 9 Pclaf
## 97 PDE5A1 0.000e+00 1.9600 0.914 0.169 0.000e+00 7 Pde5a
## 98 PDZD22 0.000e+00 2.2018 0.934 0.210 0.000e+00 7 Pdzd2
## 99 PLTP 0.000e+00 0.9974 0.951 0.452 0.000e+00 0 Pltp
## 100 PLVAP 0.000e+00 1.1891 0.880 0.373 0.000e+00 0 Plvap
## 101 POSTN2 0.000e+00 2.8807 0.965 0.090 0.000e+00 7 Postn
## 102 PRDX61 0.000e+00 2.8521 1.000 0.681 0.000e+00 6 Prdx6
## 103 RAC2 0.000e+00 1.4305 0.961 0.277 0.000e+00 2 Rac2
## 104 RAC22 1.460e-234 0.7545 0.915 0.355 3.165e-230 10 Rac2
## 105 RAMP2 0.000e+00 1.0279 0.976 0.434 0.000e+00 0 Ramp2
## 106 RARRES23 0.000e+00 3.2696 0.966 0.295 0.000e+00 8 Rarres2
## 107 RETNLA3 8.896e-125 3.9550 0.398 0.045 1.929e-120 13 Retnla
## 108 RPL12 0.000e+00 0.8919 0.974 0.860 0.000e+00 1 Rpl12
## 109 S100A101 0.000e+00 2.8948 0.968 0.696 0.000e+00 9 S100a10
## 110 SCGB1A1 0.000e+00 3.1611 0.697 0.425 0.000e+00 3 Scgb1a1
## 111 SCGB1A12 5.123e-84 6.0978 0.888 0.454 1.110e-79 13 Scgb1a1
## 112 SCGB3A2 1.762e-245 3.6370 0.260 0.010 3.819e-241 13 Scgb3a2
## 113 SELP 1.004e-218 1.6324 0.305 0.024 2.175e-214 12 Selp
## 114 SELPLG2 1.603e-194 0.7257 0.805 0.315 3.476e-190 10 Selplg
## 115 SERPING1 0.000e+00 1.9592 0.962 0.298 0.000e+00 4 Serping1
## 116 SFTPA1 0.000e+00 3.5185 0.997 0.190 0.000e+00 3 Sftpa1
## 117 SFTPA11 0.000e+00 2.4925 0.988 0.220 0.000e+00 5 Sftpa1
## 118 SFTPB 0.000e+00 3.5137 0.997 0.149 0.000e+00 3 Sftpb
## 119 SFTPB1 0.000e+00 2.2938 0.988 0.181 0.000e+00 5 Sftpb
## 120 SFTPC 0.000e+00 3.6318 0.994 0.523 0.000e+00 3 Sftpc
## 121 SFTPC1 0.000e+00 2.5258 0.992 0.540 0.000e+00 5 Sftpc
## 122 SFTPD 0.000e+00 3.2333 0.995 0.131 0.000e+00 3 Sftpd
## 123 SFTPD1 0.000e+00 2.3517 0.991 0.162 0.000e+00 5 Sftpd
## 124 SLC34A21 0.000e+00 2.1997 0.971 0.142 0.000e+00 5 Slc34a2
## 125 SOD3 0.000e+00 2.2334 0.961 0.245 0.000e+00 4 Sod3
## 126 SRGN2 1.925e-161 0.8425 0.960 0.681 4.173e-157 10 Srgn
## 127 STK17B 0.000e+00 0.8169 0.741 0.275 0.000e+00 1 Stk17b
## 128 TCF21 0.000e+00 2.1114 0.898 0.215 0.000e+00 4 Tcf21
## 129 TIMP1 0.000e+00 3.5343 0.477 0.057 0.000e+00 8 Timp1
## 130 TM4SF12 1.293e-21 2.2322 0.665 0.447 2.802e-17 12 Tm4sf1
## 131 TOP2A 0.000e+00 2.9707 0.829 0.016 0.000e+00 9 Top2a
## 132 TPPP32 1.775e-228 3.7031 0.990 0.220 3.847e-224 13 Tppp3
## 133 TSPAN7 0.000e+00 1.1471 0.959 0.412 0.000e+00 0 Tspan7
## 134 TUBB5 0.000e+00 3.1306 0.989 0.716 0.000e+00 9 Tubb5
## 135 TYROBP 0.000e+00 3.5928 0.937 0.016 0.000e+00 11 Tyrobp
## 136 UBE2C 0.000e+00 3.1423 0.723 0.007 0.000e+00 9 Ube2c
## 137 VCAM12 1.645e-14 1.9137 0.280 0.147 3.565e-10 12 Vcam1
## 138 VPS37B 0.000e+00 0.9892 0.779 0.297 0.000e+00 1 Vps37b
## 139 VPS37B3 1.298e-150 0.7941 0.765 0.362 2.814e-146 10 Vps37b
## 140 VWF1 5.878e-51 2.2979 0.375 0.128 1.274e-46 12 Vwf
## description
## 1 actin alpha 2, smooth muscle [Source:HGNC Symbol;Acc:HGNC:130]
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 BTG anti-proliferation factor 1 [Source:HGNC Symbol;Acc:HGNC:1130]
## 6 <NA>
## 7 <NA>
## 8 <NA>
## 9 C-C motif chemokine ligand 11 [Source:HGNC Symbol;Acc:HGNC:10610]
## 10 <NA>
## 11 <NA>
## 12 C-C motif chemokine ligand 5 [Source:HGNC Symbol;Acc:HGNC:10632]
## 13 <NA>
## 14 <NA>
## 15 C-C motif chemokine receptor 7 [Source:HGNC Symbol;Acc:HGNC:1608]
## 16 <NA>
## 17 <NA>
## 18 <NA>
## 19 CD3g molecule [Source:HGNC Symbol;Acc:HGNC:1675]
## 20 <NA>
## 21 CD52 molecule [Source:HGNC Symbol;Acc:HGNC:1804]
## 22 <NA>
## 23 <NA>
## 24 <NA>
## 25 <NA>
## 26 claudin 5 [Source:HGNC Symbol;Acc:HGNC:2047]
## 27 C-type lectin domain containing 14A [Source:HGNC Symbol;Acc:HGNC:19832]
## 28 <NA>
## 29 <NA>
## 30 <NA>
## 31 <NA>
## 32 <NA>
## 33 <NA>
## 34 <NA>
## 35 <NA>
## 36 cathepsin S [Source:HGNC Symbol;Acc:HGNC:2545]
## 37 <NA>
## 38 <NA>
## 39 C-X-C motif chemokine ligand 2 [Source:HGNC Symbol;Acc:HGNC:4603]
## 40 cytochrome b-245 beta chain [Source:HGNC Symbol;Acc:HGNC:2578]
## 41 <NA>
## 42 <NA>
## 43 <NA>
## 44 decorin [Source:HGNC Symbol;Acc:HGNC:2705]
## 45 dual specificity phosphatase 10 [Source:HGNC Symbol;Acc:HGNC:3065]
## 46 EGF like domain multiple 7 [Source:HGNC Symbol;Acc:HGNC:20594]
## 47 <NA>
## 48 <NA>
## 49 Fc fragment of IgE receptor Ig [Source:HGNC Symbol;Acc:HGNC:3611]
## 50 four and a half LIM domains 1 [Source:HGNC Symbol;Acc:HGNC:3702]
## 51 fibronectin 1 [Source:HGNC Symbol;Acc:HGNC:3778]
## 52 glycosylphosphatidylinositol anchored high density lipoprotein binding protein 1 [Source:HGNC Symbol;Acc:HGNC:24945]
## 53 glutathione peroxidase 3 [Source:HGNC Symbol;Acc:HGNC:4555]
## 54 <NA>
## 55 gelsolin [Source:HGNC Symbol;Acc:HGNC:4620]
## 56 <NA>
## 57 <NA>
## 58 <NA>
## 59 granzyme A [Source:HGNC Symbol;Acc:HGNC:4708]
## 60 <NA>
## 61 high mobility group box 2 [Source:HGNC Symbol;Acc:HGNC:5000]
## 62 <NA>
## 63 <NA>
## 64 15-hydroxyprostaglandin dehydrogenase [Source:HGNC Symbol;Acc:HGNC:5154]
## 65 <NA>
## 66 <NA>
## 67 insulin like growth factor binding protein 5 [Source:HGNC Symbol;Acc:HGNC:5474]
## 68 interleukin 1 beta [Source:HGNC Symbol;Acc:HGNC:5992]
## 69 interleukin 7 receptor [Source:HGNC Symbol;Acc:HGNC:6024]
## 70 <NA>
## 71 indolethylamine N-methyltransferase [Source:HGNC Symbol;Acc:HGNC:6069]
## 72 <NA>
## 73 keratin 71 [Source:HGNC Symbol;Acc:HGNC:28927]
## 74 lysosomal protein transmembrane 5 [Source:HGNC Symbol;Acc:HGNC:29612]
## 75 <NA>
## 76 lymphoid enhancer binding factor 1 [Source:HGNC Symbol;Acc:HGNC:6551]
## 77 <NA>
## 78 <NA>
## 79 lymphotoxin beta [Source:HGNC Symbol;Acc:HGNC:6711]
## 80 <NA>
## 81 <NA>
## 82 <NA>
## 83 <NA>
## 84 <NA>
## 85 microfibril associated protein 4 [Source:HGNC Symbol;Acc:HGNC:7035]
## 86 <NA>
## 87 <NA>
## 88 <NA>
## 89 marker of proliferation Ki-67 [Source:HGNC Symbol;Acc:HGNC:7107]
## 90 multimerin 1 [Source:HGNC Symbol;Acc:HGNC:7178]
## 91 <NA>
## 92 <NA>
## 93 <NA>
## 94 <NA>
## 95 natural killer cell granule protein 7 [Source:HGNC Symbol;Acc:HGNC:7830]
## 96 PCNA clamp associated factor [Source:HGNC Symbol;Acc:HGNC:28961]
## 97 <NA>
## 98 <NA>
## 99 phospholipid transfer protein [Source:HGNC Symbol;Acc:HGNC:9093]
## 100 plasmalemma vesicle associated protein [Source:HGNC Symbol;Acc:HGNC:13635]
## 101 <NA>
## 102 <NA>
## 103 Rac family small GTPase 2 [Source:HGNC Symbol;Acc:HGNC:9802]
## 104 <NA>
## 105 receptor activity modifying protein 2 [Source:HGNC Symbol;Acc:HGNC:9844]
## 106 <NA>
## 107 <NA>
## 108 ribosomal protein L12 [Source:HGNC Symbol;Acc:HGNC:10302]
## 109 <NA>
## 110 secretoglobin family 1A member 1 [Source:HGNC Symbol;Acc:HGNC:12523]
## 111 <NA>
## 112 secretoglobin family 3A member 2 [Source:HGNC Symbol;Acc:HGNC:18391]
## 113 selectin P [Source:HGNC Symbol;Acc:HGNC:10721]
## 114 <NA>
## 115 serpin family G member 1 [Source:HGNC Symbol;Acc:HGNC:1228]
## 116 surfactant protein A1 [Source:HGNC Symbol;Acc:HGNC:10798]
## 117 <NA>
## 118 surfactant protein B [Source:HGNC Symbol;Acc:HGNC:10801]
## 119 <NA>
## 120 surfactant protein C [Source:HGNC Symbol;Acc:HGNC:10802]
## 121 <NA>
## 122 surfactant protein D [Source:HGNC Symbol;Acc:HGNC:10803]
## 123 <NA>
## 124 <NA>
## 125 superoxide dismutase 3 [Source:HGNC Symbol;Acc:HGNC:11181]
## 126 <NA>
## 127 serine/threonine kinase 17b [Source:HGNC Symbol;Acc:HGNC:11396]
## 128 transcription factor 21 [Source:HGNC Symbol;Acc:HGNC:11632]
## 129 TIMP metallopeptidase inhibitor 1 [Source:HGNC Symbol;Acc:HGNC:11820]
## 130 <NA>
## 131 DNA topoisomerase II alpha [Source:HGNC Symbol;Acc:HGNC:11989]
## 132 <NA>
## 133 tetraspanin 7 [Source:HGNC Symbol;Acc:HGNC:11854]
## 134 <NA>
## 135 transmembrane immune signaling adaptor TYROBP [Source:HGNC Symbol;Acc:HGNC:12449]
## 136 ubiquitin conjugating enzyme E2 C [Source:HGNC Symbol;Acc:HGNC:15937]
## 137 <NA>
## 138 VPS37B subunit of ESCRT-I [Source:HGNC Symbol;Acc:HGNC:25754]
## 139 <NA>
## 140 <NA>
## ensembl_gene_id
## 1 ENSG00000107796
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 ENSG00000133639
## 6 <NA>
## 7 <NA>
## 8 <NA>
## 9 ENSG00000172156
## 10 <NA>
## 11 <NA>
## 12 ENSG00000271503
## 13 <NA>
## 14 <NA>
## 15 ENSG00000126353
## 16 <NA>
## 17 <NA>
## 18 <NA>
## 19 ENSG00000160654
## 20 <NA>
## 21 ENSG00000169442
## 22 <NA>
## 23 <NA>
## 24 <NA>
## 25 <NA>
## 26 ENSG00000184113
## 27 ENSG00000176435
## 28 <NA>
## 29 <NA>
## 30 <NA>
## 31 <NA>
## 32 <NA>
## 33 <NA>
## 34 <NA>
## 35 <NA>
## 36 ENSG00000163131
## 37 <NA>
## 38 <NA>
## 39 ENSG00000081041
## 40 ENSG00000165168
## 41 <NA>
## 42 <NA>
## 43 <NA>
## 44 ENSG00000011465
## 45 ENSG00000143507
## 46 ENSG00000172889
## 47 <NA>
## 48 <NA>
## 49 ENSG00000158869
## 50 ENSG00000022267
## 51 ENSG00000115414
## 52 ENSG00000277494
## 53 ENSG00000211445
## 54 <NA>
## 55 ENSG00000148180
## 56 <NA>
## 57 <NA>
## 58 <NA>
## 59 ENSG00000145649
## 60 <NA>
## 61 ENSG00000164104
## 62 <NA>
## 63 <NA>
## 64 ENSG00000164120
## 65 <NA>
## 66 <NA>
## 67 ENSG00000115461
## 68 ENSG00000125538
## 69 ENSG00000168685
## 70 <NA>
## 71 ENSG00000241644
## 72 <NA>
## 73 ENSG00000139648
## 74 ENSG00000162511
## 75 <NA>
## 76 ENSG00000138795
## 77 <NA>
## 78 <NA>
## 79 ENSG00000227507
## 80 <NA>
## 81 <NA>
## 82 <NA>
## 83 <NA>
## 84 <NA>
## 85 ENSG00000166482
## 86 <NA>
## 87 <NA>
## 88 <NA>
## 89 ENSG00000148773
## 90 ENSG00000138722
## 91 <NA>
## 92 <NA>
## 93 <NA>
## 94 <NA>
## 95 ENSG00000105374
## 96 ENSG00000166803
## 97 <NA>
## 98 <NA>
## 99 ENSG00000100979
## 100 ENSG00000130300
## 101 <NA>
## 102 <NA>
## 103 ENSG00000128340
## 104 <NA>
## 105 ENSG00000131477
## 106 <NA>
## 107 <NA>
## 108 ENSG00000197958
## 109 <NA>
## 110 ENSG00000149021
## 111 <NA>
## 112 ENSG00000164265
## 113 ENSG00000174175
## 114 <NA>
## 115 ENSG00000149131
## 116 ENSG00000122852
## 117 <NA>
## 118 ENSG00000168878
## 119 <NA>
## 120 ENSG00000168484
## 121 <NA>
## 122 ENSG00000133661
## 123 <NA>
## 124 <NA>
## 125 ENSG00000109610
## 126 <NA>
## 127 ENSG00000081320
## 128 ENSG00000118526
## 129 ENSG00000102265
## 130 <NA>
## 131 ENSG00000131747
## 132 <NA>
## 133 ENSG00000156298
## 134 <NA>
## 135 ENSG00000011600
## 136 ENSG00000175063
## 137 <NA>
## 138 ENSG00000139722
## 139 <NA>
## 140 <NA>
<- count_clonotype_by_cluster(cluster_scd) cluster_scd
## This result may be found in the res0p2_clusters element of scd@misc.
@misc$res0p2_clusters cluster_scd
## cluster_name cells has_clono proportion
## 1 0 7721 73 0.00945473384276648
## 2 1 6101 5162 0.846090804786101
## 3 2 5187 4311 0.831116252168884
## 4 3 4475 240 0.0536312849162011
## 5 4 3331 277 0.0831582107475233
## 6 5 3189 193 0.0605205393540295
## 7 6 2680 251 0.0936567164179104
## 8 7 2213 143 0.0646181653863534
## 9 8 1023 59 0.0576735092864125
## 10 9 920 631 0.685869565217391
## 11 10 911 706 0.774972557628979
## 12 11 812 83 0.102216748768473
## 13 12 325 12 0.0369230769230769
## 14 13 196 25 0.127551020408163
## 15 cluster_sum 39084 12166 0.31127827243885
<- cluster_scd@misc$res0p2_clusters[["proportion"]] >= 0.75
high_vdj_cluster_idx <- cluster_scd@misc$res0p2_clusters[high_vdj_cluster_idx, "cluster_name"] high_vdj_clusters
Lets add a new annotation column ‘high_vdj_prop’; this is the set of all cells in clusters with a proportion of vdj’s >= 0.75.
<- cluster_scd@meta.data[["res0p2_clusters"]] %in% high_vdj_clusters
high_vdj_cell_idx summary(high_vdj_cell_idx)
## Mode FALSE TRUE
## logical 26885 12199
"high_vdj_cluster"]] <- "No"
cluster_scd[[@meta.data[high_vdj_cell_idx, "high_vdj_cluster"] <- "Yes"
cluster_scd
@meta.data[["sample_vdjcluster"]] <- paste0(cluster_scd@meta.data[["orig.ident"]], "_",
cluster_scd@meta.data[["high_vdj_cluster"]]) cluster_scd
This asks what is different from all cells in high VDJ clusters vs. all others. It is likely we want to concatenate this will sample ID in order to perform specific comparisons.
<- cluster_scd
highvdj_scd Idents(highvdj_scd) <- highvdj_scd[["high_vdj_cluster"]]
<- FindAllMarkers(highvdj_scd, only.pos=TRUE, logfc.threshold=fc_threshold) cluster_markers
## Calculating cluster Yes
## Calculating cluster No
<- as.data.frame(cluster_markers)
cluster_genes rownames(cluster_genes) <- toupper(rownames(cluster_genes))
<- merge(cluster_genes, brief, by.x="row.names", by.y="hgnc_symbol",
annotated_clusters all.x=TRUE)
head(annotated_clusters)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## 1 1810037I17RIK 0 0.9773 0.698 0.366 0 No 1810037I17Rik
## 2 4931406P16RIK 0 0.7634 0.452 0.078 0 No 4931406P16Rik
## 3 ABCA3 0 1.1392 0.403 0.058 0 No Abca3
## 4 ACAA1A 0 0.6541 0.431 0.117 0 No Acaa1a
## 5 ACAA2 0 0.6983 0.493 0.189 0 No Acaa2
## 6 ACADL 0 0.8614 0.566 0.162 0 No Acadl
## description
## 1 <NA>
## 2 <NA>
## 3 ATP binding cassette subfamily A member 3 [Source:HGNC Symbol;Acc:HGNC:33]
## 4 <NA>
## 5 acetyl-CoA acyltransferase 2 [Source:HGNC Symbol;Acc:HGNC:83]
## 6 acyl-CoA dehydrogenase long chain [Source:HGNC Symbol;Acc:HGNC:88]
## ensembl_gene_id
## 1 <NA>
## 2 <NA>
## 3 ENSG00000167972
## 4 <NA>
## 5 ENSG00000167315
## 6 ENSG00000115361
"highvdjcluster_vs_all"]] <- annotated_clusters de_lst[[
A few clusters have a large majority of the VDJ information. Previously I recorded which they are, but because of the vagaries of UMAP, this appears to change from invocation to invocation.
Among the cluster2_samples, cluster 1 has 85% VDJ’s; so let us look at mock vs. nasal and mock vs. muscular in it.
I do an explicit query of my clonotype percentages in order to extract the cluster with the most VDJ markers. I will then perform an explicit pairwise comparison of what I perceive to be the most interesting categories therein: muscular/mock, nasal/mock, and muscular/nasal.
<- order(cluster_scd@misc[["res0p2_clusters"]][["proportion"]], decreasing=TRUE)
highest_cluster_idx <- cluster_scd@misc[["res0p2_clusters"]][highest_cluster_idx[1], "cluster_name"]
highest_cluster
<- paste0("mock_", highest_cluster)
c1_mock <- paste0("n_", highest_cluster)
c1_nasal <- paste0("m_", highest_cluster)
c1_musc <- "cluster2_sample"
category <- cluster_scd[["cluster2_sample"]] == c1_mock
mock_group sum(mock_group)
## [1] 1052
<- cluster_scd[["cluster2_sample"]] == c1_musc
musc_group sum(musc_group)
## [1] 2437
<- cluster_scd[["cluster2_sample"]] == c1_nasal
nasal_group sum(nasal_group)
## [1] 1556
<- FindMarkers(
nasal_clust1 group.by=category,
cluster_scd, ident.1=c1_nasal, ident.2=c1_mock)
<- as.data.frame(nasal_clust1)
nasal_clust1 rownames(nasal_clust1) <- toupper(rownames(nasal_clust1))
<- merge(nasal_clust1, brief, by="row.names", by.y="hgnc_symbol",
nasal_clust1 all.x=TRUE)
head(nasal_clust1)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 AY036118 3.390e-06 -0.2744 0.596 0.672 7.350e-02
## 2 BST2 2.210e-09 -0.3041 0.408 0.497 4.790e-05
## 3 DDX5 1.286e-14 -0.2630 0.843 0.885 2.788e-10
## 4 FAU 9.608e-32 0.3347 0.998 0.998 2.083e-27
## 5 FOS 8.581e-01 -0.2732 0.548 0.532 1.000e+00
## 6 GADD45B 3.533e-03 -0.3465 0.246 0.290 1.000e+00
## description
## 1 <NA>
## 2 bone marrow stromal cell antigen 2 [Source:HGNC Symbol;Acc:HGNC:1119]
## 3 DEAD-box helicase 5 [Source:HGNC Symbol;Acc:HGNC:2746]
## 4 FAU ubiquitin like and ribosomal protein S30 fusion [Source:HGNC Symbol;Acc:HGNC:3597]
## 5 Fos proto-oncogene, AP-1 transcription factor subunit [Source:HGNC Symbol;Acc:HGNC:3796]
## 6 growth arrest and DNA damage inducible beta [Source:HGNC Symbol;Acc:HGNC:4096]
## ensembl_gene_id
## 1 <NA>
## 2 ENSG00000130303
## 3 ENSG00000108654
## 4 ENSG00000149806
## 5 ENSG00000170345
## 6 ENSG00000099860
"highestvdj_nasal_vs_mock"]] <- nasal_clust1 de_lst[[
<- FindMarkers(
musc_clust1 group.by=category,
cluster_scd, ident.1=c1_musc, ident.2=c1_mock)
<- as.data.frame(musc_clust1)
musc_clust1 rownames(musc_clust1) <- toupper(rownames(musc_clust1))
<- merge(musc_clust1, brief, by="row.names", by.y="hgnc_symbol",
musc_clust1 all.x=TRUE)
head(musc_clust1)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 ATP5E 2.682e-19 0.2666 0.568 0.429 5.814e-15
## 2 BST2 1.366e-17 -0.3586 0.370 0.497 2.962e-13
## 3 CCL5 6.684e-05 -0.8266 0.068 0.106 1.000e+00
## 4 CD52 1.630e-28 0.4050 0.797 0.661 3.534e-24
## 5 COX8A 6.312e-16 0.2611 0.689 0.574 1.368e-11
## 6 FAU 5.237e-105 0.5753 1.000 0.998 1.135e-100
## description
## 1 <NA>
## 2 bone marrow stromal cell antigen 2 [Source:HGNC Symbol;Acc:HGNC:1119]
## 3 C-C motif chemokine ligand 5 [Source:HGNC Symbol;Acc:HGNC:10632]
## 4 CD52 molecule [Source:HGNC Symbol;Acc:HGNC:1804]
## 5 cytochrome c oxidase subunit 8A [Source:HGNC Symbol;Acc:HGNC:2294]
## 6 FAU ubiquitin like and ribosomal protein S30 fusion [Source:HGNC Symbol;Acc:HGNC:3597]
## ensembl_gene_id
## 1 <NA>
## 2 ENSG00000130303
## 3 ENSG00000271503
## 4 ENSG00000169442
## 5 ENSG00000176340
## 6 ENSG00000149806
"highestvdj_muscular_vs_mock"]] <- musc_clust1 de_lst[[
<- FindMarkers(
nm_clust1 group.by=category,
cluster_scd, ident.1=c1_musc, ident.2=c1_nasal)
<- as.data.frame(nm_clust1)
nm_clust1 rownames(nm_clust1) <- toupper(rownames(nm_clust1))
<- merge(nm_clust1, brief, by="row.names", by.y="hgnc_symbol",
nm_clust1 all.x=TRUE)
head(nm_clust1)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 AY036118 1.264e-13 0.3009 0.682 0.596 2.740e-09
## 2 CCL5 5.210e-18 -0.9787 0.068 0.151 1.129e-13
## 3 CCR7 7.185e-13 0.2799 0.728 0.629 1.558e-08
## 4 DDX5 7.724e-25 0.2788 0.897 0.843 1.675e-20
## 5 EEF2 9.395e-21 0.2690 0.988 0.976 2.037e-16
## 6 GM42418 8.020e-27 0.4369 0.998 0.995 1.739e-22
## description
## 1 <NA>
## 2 C-C motif chemokine ligand 5 [Source:HGNC Symbol;Acc:HGNC:10632]
## 3 C-C motif chemokine receptor 7 [Source:HGNC Symbol;Acc:HGNC:1608]
## 4 DEAD-box helicase 5 [Source:HGNC Symbol;Acc:HGNC:2746]
## 5 eukaryotic translation elongation factor 2 [Source:HGNC Symbol;Acc:HGNC:3214]
## 6 <NA>
## ensembl_gene_id
## 1 <NA>
## 2 ENSG00000271503
## 3 ENSG00000126353
## 4 ENSG00000108654
## 5 ENSG00000167658
## 6 <NA>
"highestvdj_muscular_vs_nasal"]] <- musc_clust1 de_lst[[
Conversely, we can just ask for differences among all vdj cells.
"vdj_sample"]] <- ""
cluster_scd[[<- paste0(cluster_scd@meta.data[["orig.ident"]], "_",
cluster_scd_vector @meta.data[["tcell"]])
cluster_scd"vdj_sample"]] <- cluster_scd_vector
cluster_scd[[
<- "mock_Yes"
v_mock <- "n_Yes"
v_nasal <- "m_Yes"
v_musc <- "vdj_sample"
category <- cluster_scd[[category]] == v_mock
mock_group sum(mock_group)
## [1] 3029
<- cluster_scd[[category]] == v_musc
musc_group sum(musc_group)
## [1] 4167
<- cluster_scd[[category]] == v_nasal
nasal_group sum(nasal_group)
## [1] 3020
<- cluster_scd
vdj_cluster Idents(vdj_cluster) <- cluster_scd[[category]]
<- NormalizeData(vdj_cluster) %>%
vdj_cluster FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors() %>%
FindClusters() %>%
RunTSNE() %>%
RunUMAP(reduction="pca", dims=1:10)
## Centering and scaling data matrix
## PC_ 1
## Positive: Sparc, Emp2, Cd63, Igfbp7, Mettl7a1, Ager, Gsn, Cldn18, Timp3, Limch1
## Cst3, Gpx3, Cystm1, Selenbp1, Cryab, Bgn, Npnt, Selenbp2, Bcam, Pmp22
## Clic3, Apoe, Hspb1, Krt7, Ces1d, Aqp5, Anxa3, Ifitm3, Alcam, Fhl1
## Negative: Rac2, Coro1a, Arhgdib, Laptm5, Cd3e, Ms4a4b, Selplg, Cd52, Cd3g, Ptprc
## Cd3d, Lat, Lck, Ms4a6b, Thy1, Ltb, Ptprcap, Sept1, Il7r, Itgb7
## Lsp1, Gramd3, Fyb, Cd2, Cd28, Icos, Cd69, Bin2, Il2rb, AW112010
## PC_ 2
## Positive: Bgn, Serping1, Sod3, Plpp3, Prelp, Rarres2, Col1a2, Sparcl1, Hsd11b1, Vim
## Olfml3, Col1a1, Loxl1, Tcf21, Clec3b, Ptgis, Mxra8, Spon1, Plxdc2, Fxyd1
## Gpx3, Fhl1, Ltbp4, Col3a1, Lrp1, Ppp1r14a, Cdo1, Pcolce2, Mfap4, Colec12
## Negative: Sftpd, Napsa, Wfdc2, Sftpb, Cxcl15, Slc34a2, Cbr2, Sftpa1, Cldn3, Atp1b1
## Chil1, Lamp3, Muc1, Ppp1r14c, Sftpc, Lgi3, Lyz2, Ctsh, Dram1, Scd1
## Hc, Car8, Egfl6, Lpcat1, S100g, Irx1, Pi4k2b, Abca3, Pla2g1b, Ank3
## PC_ 3
## Positive: Rtkn2, Scnn1g, Hopx, Cyp4b1, Igfbp2, Clic5, Spock2, Col4a3, Sema3e, Pdpn
## Cyp2b10, Scnn1b, Tspan8, Col4a4, Lama3, Flrt3, Cytl1, Krt7, Tacstd2, Sec14l3
## Cdkn2b, Lamb3, Itgb6, Krt19, Mab21l4, Pakap.1, Tmem37, Scnn1a, Lmo7, Akap5
## Negative: Mgst1, Col6a1, Nupr1, Loxl1, Col3a1, Col1a1, Sod3, Col1a2, C3, Dram1
## Ces1d, Col6a2, Apoe, Slc34a2, Mmp2, Plxdc2, Cxcl15, Pdgfra, Dpep1, Chil1
## Lamp3, Sftpb, Spon1, Prelp, Cd302, Sftpd, Olfml3, Lyz2, Atp1a2, Sftpa1
## PC_ 4
## Positive: Cldn5, Plvap, Hpgd, H2-Ab1, Cd93, Cd74, Tmem252, H2-Aa, H2-Eb1, Sema3c
## Pcdh1, Mcam, Gja4, Efnb2, Sox17, Dll4, Lyve1, Stmn2, Tm4sf1, Hspa1a
## Kdr, Ly6c1, Plk2, Hilpda, Tmcc2, Emcn, Hspa1b, Scn7a, Ly6a, Ifitm3
## Negative: S100a6, Lgals1, Birc5, Lgals3, Mki67, Anxa1, Pclaf, Prdx6, Tacstd2, Top2a
## Pdpn, Cyp2b10, Sema3e, Ccna2, Krt19, Scnn1g, Coro1a, Ube2c, Spock2, Rac2
## Scnn1b, S100a11, Rtkn2, Cdk1, Mab21l4, Arhgdig, Col4a4, Ccnb2, Rrm2, Col4a3
## PC_ 5
## Positive: Fam183b, Foxj1, Tmem212, 1700007K13Rik, 1110017D15Rik, Ccdc153, Gm19935, Cfap126, Spaca9, Tctex1d4
## BC051019, 1700094D03Rik, Odf3b, 1700001C02Rik, 1700016K19Rik, Tekt4, Dnali1, Rsph1, AU040972, Cyp2s1
## Ak7, Ccdc113, Sntn, Gm867, Nme5, Nme9, 2410004P03Rik, Dynlrb2, Pifo, Lrrc10b
## Negative: Col4a4, Aqp5, Scnn1g, Ager, Rtkn2, Ndnf, Col4a3, Slc39a8, Gprc5a, Spock2
## Tmem37, Pdpn, Scnn1b, Igfbp6, Igfbp2, Flrt3, Crlf1, Scd2, Lamb3, Fads3
## Cldn18, Rgcc, Itgb6, Mal2, Tmod1, Fam189a2, Rnase4, Hs2st1, Abca5, Lamc2
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 39084
## Number of edges: 1210541
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9102
## Number of communities: 25
## Elapsed time: 11 seconds
## 11:17:36 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:17:36 Read 39084 rows and found 10 numeric columns
## 11:17:36 Using Annoy for neighbor search, n_neighbors = 30
## 11:17:36 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:17:43 Writing NN index file to temp file /tmp/RtmpBNmUfO/filee94a2feb10c1
## 11:17:43 Searching Annoy index using 1 thread, search_k = 3000
## 11:18:05 Annoy recall = 100%
## 11:18:07 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:18:13 Initializing from normalized Laplacian + noise (using irlba)
## 11:18:21 Commencing optimization for 200 epochs, with 1600638 positive edges
## 11:18:44 Optimization finished
<- DimPlot(vdj_cluster, reduction="umap", group.by="ident", label=TRUE)
plotted plotted
DimPlot(vdj_cluster, reduction="tsne")
<- FindMarkers(
nasal_v group.by=category,
cluster_scd, ident.1=v_nasal, ident.2=v_mock)
<- as.data.frame(nasal_v)
nasal_v rownames(nasal_v) <- toupper(rownames(nasal_v))
<- merge(nasal_v, brief, by="row.names", by.y="hgnc_symbol",
nasal_v all.x=TRUE)
head(nasal_v)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 ANXA2 4.860e-02 0.2751 0.322 0.356 1.000e+00
## 2 ARL6IP1 1.968e-02 0.2702 0.689 0.728 1.000e+00
## 3 ATP5E 5.470e-27 0.3738 0.620 0.513 1.186e-22
## 4 ATP5G1 1.089e-15 0.3410 0.514 0.433 2.361e-11
## 5 ATP5G2 5.228e-18 0.3742 0.809 0.754 1.133e-13
## 6 ATP5G3 4.975e-03 0.3538 0.726 0.700 1.000e+00
## description
## 1 annexin A2 [Source:HGNC Symbol;Acc:HGNC:537]
## 2 ADP ribosylation factor like GTPase 6 interacting protein 1 [Source:HGNC Symbol;Acc:HGNC:697]
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## ensembl_gene_id
## 1 ENSG00000182718
## 2 ENSG00000170540
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
"highvdj_nasal_vs_mock"]] <- nasal_v
de_lst[[
<- FindMarkers(
musc_v group.by=category,
cluster_scd, ident.1=v_musc, ident.2=v_mock)
<- as.data.frame(musc_v)
musc_v rownames(musc_v) <- toupper(rownames(musc_v))
<- merge(musc_v, brief, by="row.names", by.y="hgnc_symbol",
musc_v all.x=TRUE)
head(musc_v)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 1810037I17RIK 1.110e-07 0.2520 0.464 0.422 2.407e-03
## 2 2410006H16RIK 3.189e-33 0.2746 0.536 0.405 6.914e-29
## 3 ACTB 2.581e-07 0.2547 0.998 0.999 5.595e-03
## 4 ACTG1 2.451e-03 0.2753 0.982 0.987 1.000e+00
## 5 AGER 1.620e-02 0.2952 0.137 0.118 1.000e+00
## 6 ANP32B 3.749e-10 0.3270 0.528 0.463 8.127e-06
## description
## 1 <NA>
## 2 <NA>
## 3 actin beta [Source:HGNC Symbol;Acc:HGNC:132]
## 4 actin gamma 1 [Source:HGNC Symbol;Acc:HGNC:144]
## 5 advanced glycosylation end-product specific receptor [Source:HGNC Symbol;Acc:HGNC:320]
## 6 acidic nuclear phosphoprotein 32 family member B [Source:HGNC Symbol;Acc:HGNC:16677]
## ensembl_gene_id
## 1 <NA>
## 2 <NA>
## 3 ENSG00000075624
## 4 ENSG00000184009
## 5 ENSG00000204305
## 6 ENSG00000136938
"highvdj_muscular_vs_mock"]] <- musc_v
de_lst[[
<- FindMarkers(
nm_v group.by=category,
cluster_scd, ident.1=v_musc, ident.2=v_nasal)
<- as.data.frame(nm_v)
nm_v rownames(nm_v) <- toupper(rownames(nm_v))
<- merge(nm_v, brief, by="row.names", by.y="hgnc_symbol",
nm_v all.x=TRUE)
head(nm_v)
## Row.names p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 ACE 2.707e-12 0.2564 0.142 0.087 5.869e-08
## 2 ATF3 1.435e-08 0.3619 0.207 0.157 3.110e-04
## 3 BST2 3.935e-11 -0.6438 0.546 0.608 8.530e-07
## 4 CALCRL 3.474e-12 0.2991 0.159 0.103 7.531e-08
## 5 CAV1 5.676e-12 0.2568 0.154 0.098 1.230e-07
## 6 CCL5 4.496e-50 -0.8236 0.162 0.309 9.747e-46
## description
## 1 angiotensin I converting enzyme [Source:HGNC Symbol;Acc:HGNC:2707]
## 2 activating transcription factor 3 [Source:HGNC Symbol;Acc:HGNC:785]
## 3 bone marrow stromal cell antigen 2 [Source:HGNC Symbol;Acc:HGNC:1119]
## 4 calcitonin receptor like receptor [Source:HGNC Symbol;Acc:HGNC:16709]
## 5 caveolin 1 [Source:HGNC Symbol;Acc:HGNC:1527]
## 6 C-C motif chemokine ligand 5 [Source:HGNC Symbol;Acc:HGNC:10632]
## ensembl_gene_id
## 1 ENSG00000159640
## 2 ENSG00000162772
## 3 ENSG00000130303
## 4 ENSG00000064989
## 5 ENSG00000105974
## 6 ENSG00000271503
"highvdj_muscular_vs_nasal"]] <- musc_v de_lst[[
<- write_xlsx(
written
de_lst,excel=glue::glue("excel/scd_de_results-v{ver}.xlsx"))
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## [1] "vdj_vs_all"
## [1] "vdj_vs_all"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [1] "vdj_vs_all" "samples_vs_all"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all" "samplecluster_vs_all"
## [1] "vdj_vs_all" "samples_vs_all" "samplecluster_vs_all"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all" "samplecluster_vs_all"
## [4] "highvdjcluster_vs_all"
## [1] "vdj_vs_all" "samples_vs_all" "samplecluster_vs_all"
## [4] "highvdjcluster_vs_all"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock"
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal"
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal" "highvdj_nasal_vs_mock"
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal" "highvdj_nasal_vs_mock"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal" "highvdj_nasal_vs_mock"
## [9] "highvdj_muscular_vs_mock"
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal" "highvdj_nasal_vs_mock"
## [9] "highvdj_muscular_vs_mock"
## Deleting the file excel/scd_de_results-v202301.xlsx before writing the tables.
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## Warning in .self$setColWidths(i): NAs introduced by coercion
## [1] "vdj_vs_all" "samples_vs_all"
## [3] "samplecluster_vs_all" "highvdjcluster_vs_all"
## [5] "highestvdj_nasal_vs_mock" "highestvdj_muscular_vs_mock"
## [7] "highestvdj_muscular_vs_nasal" "highvdj_nasal_vs_mock"
## [9] "highvdj_muscular_vs_mock" "highvdj_muscular_vs_nasal"
Ok, I finally figured out its logic, I weirdly assumed it is doing something different than the various ways I was concatenating clusters/samples and using FindMarkers() to compare them. Instead, this logic of this function is very similar to my own all_pairwise(); except it is performing the pairwise comparison of subsets of the primary cell identities.
DefaultAssay(filt_scd) <- "RNA"
<- FindConservedMarkers(
conserved_markers ident.1="control",
filt_scd, grouping.var="res0p2_clusters", only.pos=TRUE,
verbose=TRUE)
## Testing group 13: (control) vs (mock, m, n)
## Testing group 12: (control) vs (mock, m, n)
## Testing group 11: (control) vs (mock, m, n)
## Testing group 10: (control) vs (mock, m, n)
## Testing group 9: (control) vs (mock, m, n)
## Testing group 8: (control) vs (mock, m, n)
## Testing group 7: (control) vs (mock, m, n)
## Testing group 6: (control) vs (mock, m, n)
## Testing group 5: (control) vs (mock, m, n)
## Testing group 4: (control) vs (mock, m, n)
## Testing group 3: (control) vs (mock, m, n)
## Testing group 2: (control) vs (mock, m, n)
## Testing group 1: (control) vs (mock, m, n)
## Testing group 0: (control) vs (mock, m, n)
<- FindMarkers(cluster_scd, group.by="orig.ident",
mock_vs_control ident.1="mock",
ident.2="control")
head(mock_vs_control)
<- FindMarkers(cluster_scd, group.by="orig.ident",
muscular_vs_mock ident.1="m",
ident.2="mock")
summary(muscular_vs_mock)
<- FindMarkers(cluster_scd, group.by="orig.ident",
nasal_vs_mock min.pct=0.25, ident.1="n",
ident.2="mock")
summary(nasal_vs_mock)
FeaturePlot(cluster_scd, features=c("Rgcc"),
split.by="orig.ident", max.cutoff=3,
cols=c("darkgreen", "darkred"))
First load the mSigDB data for mouse cell types.
<- load_gmt_signatures(signatures="reference/m8.all.v2022.1.Mm.symbols.gmt")
broad_types <- list()
broad_list for (i in names(broad_types)) {
<- geneIds(broad_types[[i]])
broad_list[[i]]
}
<- grepl(x=names(broad_list), pattern="^DESC")
descartes_idx <- broad_list[descartes_idx]
descartes_list names(descartes_list) <- gsub(x=names(descartes_list),
pattern="^DESCARTES_ORGANOGENESIS_",
replacement="")
<- grepl(x=names(broad_list), pattern="^TAB")
senis_idx <- broad_list[senis_idx]
senis_list names(senis_list) <- gsub(x=names(senis_list),
pattern = "^(TABULA_MURIS_SENIS)(_)*",
replacement = "")
<- suppressWarnings(AddModuleScore(object=cluster_scd, features=descartes_list,
cluster_scd name="descartes"))
<- suppressWarnings(AddModuleScore(object=cluster_scd, features=senis_list,
cluster_scd name="senis"))
<- summarize_scd_clusters(cluster_scd, real_column_names=names(descartes_list),
cluster_heats cluster_column="res0p2_clusters", abbreviate=FALSE,
column_prefix="descartes")
::pheatmap(cluster_heats$z_df) pheatmap
<- summarize_scd_clusters(cluster_scd, real_column_names=names(senis_list),
cluster_heats cluster_column="res0p2_clusters", abbreviate=TRUE,
column_prefix="senis")
::pheatmap(cluster_heats$z_df) pheatmap
DimPlot(fnorm_scd, label=TRUE)
<- c(151, 152, 153)
chosen_numbers <- names(senis_list)[chosen_numbers]
chosen_names chosen_names
## [1] "SPLEEN_B_CELL_AGEING"
## [2] "SPLEEN_CD4_POSITIVE_ALPHA_BETA_T_CELL_AGEING"
## [3] "SPLEEN_CD8_POSITIVE_ALPHA_BETA_T_CELL_AGEING"
<- paste0("senis", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(154, 155, 156)
chosen_numbers <- names(senis_list)[chosen_numbers]
chosen_names chosen_names
## [1] "SPLEEN_NK_CELL_AGEING" "SPLEEN_T_CELL_AGEING"
## [3] "SPLEEN_GRANULOCYTE_AGEING"
<- paste0("senis", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(70, 80)
chosen_numbers <- names(senis_list)[chosen_numbers]
chosen_names chosen_names
## [1] "LIMB_MUSCLE_MACROPHAGE_AGEING"
## [2] "LUNG_CD4_POSITIVE_ALPHA_BETA_T_CELL_AGEING"
<- paste0("senis", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(81, 88)
chosen_numbers <- names(senis_list)[chosen_numbers]
chosen_names chosen_names
## [1] "LUNG_CD8_POSITIVE_ALPHA_BETA_T_CELL_AGEING"
## [2] "LUNG_ENDOTHELIAL_CELL_OF_LYMPHATIC_VESSEL_AGEING"
<- paste0("senis", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(93)
chosen_numbers <- names(senis_list)[chosen_numbers]
chosen_names chosen_names
## [1] "LUNG_NEUTROPHIL_AGEING"
<- paste0("senis", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(2,9,20)
chosen_numbers <- names(descartes_list)[chosen_numbers]
chosen_names chosen_names
## [1] "ENDOTHELIAL_CELLS" "WHITE_BLOOD_CELLS" "STROMAL_CELLS"
<- paste0("descartes", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- c(31,36)
chosen_numbers <- names(descartes_list)[chosen_numbers]
chosen_names chosen_names
## [1] "INHIBITORY_INTERNEURONS" "NUETROPHILS"
<- paste0("descartes", chosen_numbers)
chosen_columns VlnPlot(cluster_scd, features=chosen_columns, pt.size=0, group.by="res0p2_clusters")
<- cluster_scd
cluster2_sample_scd
Idents(cluster2_sample_scd) <- cluster_scd[["cluster2_sample"]]
<- suppressWarnings(AddModuleScore(object=cluster2_sample_scd, features=descartes_list,
cluster2_sample_scd name="descartes"))
<- suppressWarnings(AddModuleScore(object=cluster2_sample_scd, features=senis_list,
cluster2_sample_scd name="senis"))
<- summarize_scd_clusters(cluster2_sample_scd, real_column_names=names(descartes_list),
cluster2_heats cluster_column="cluster2_sample", abbreviate=FALSE,
column_prefix="descartes")
::pheatmap(cluster2_heats$z_df) pheatmap
::pander(sessionInfo()) pander
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C
attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: future(v.1.31.0), tibble(v.3.1.8), GSEABase(v.1.60.0), graph(v.1.76.0), annotate(v.1.76.0), XML(v.3.99-0.13), AnnotationDbi(v.1.60.0), skimr(v.2.1.5), purrr(v.1.0.1), ggplot2(v.3.4.0), SeuratObject(v.4.1.3), Seurat(v.4.3.0), hpgltools(v.1.0), testthat(v.3.1.6), reticulate(v.1.28), SummarizedExperiment(v.1.28.0), GenomicRanges(v.1.50.2), GenomeInfoDb(v.1.34.7), IRanges(v.2.32.0), S4Vectors(v.0.36.1), MatrixGenerics(v.1.10.0), matrixStats(v.0.63.0), Biobase(v.2.58.0) and BiocGenerics(v.0.44.0)
loaded via a namespace (and not attached): ica(v.1.0-3), ps(v.1.7.2), Rsamtools(v.2.14.0), foreach(v.1.5.2), lmtest(v.0.9-40), rprojroot(v.2.0.3), crayon(v.1.5.2), rbibutils(v.2.2.13), MASS(v.7.3-58.2), nlme(v.3.1-162), backports(v.1.4.1), sva(v.3.46.0), GOSemSim(v.2.24.0), rlang(v.1.0.6), XVector(v.0.38.0), HDO.db(v.0.99.1), ROCR(v.1.0-11), irlba(v.2.3.5.1), nloptr(v.2.0.3), callr(v.3.7.3), limma(v.3.54.1), filelock(v.1.0.2), BiocParallel(v.1.32.5), rjson(v.0.2.21), bit64(v.4.0.5), glue(v.1.6.2), pheatmap(v.1.0.12), sctransform(v.0.3.5), vipor(v.0.4.5), pbkrtest(v.0.5.2), parallel(v.4.2.0), processx(v.3.8.0), spatstat.sparse(v.3.0-0), DOSE(v.3.24.2), spatstat.geom(v.3.0-6), tidyselect(v.1.2.0), usethis(v.2.1.6), fitdistrplus(v.1.1-8), variancePartition(v.1.28.3), tidyr(v.1.3.0), zoo(v.1.8-11), qqconf(v.1.3.1), GenomicAlignments(v.1.34.0), xtable(v.1.8-4), magrittr(v.2.0.3), evaluate(v.0.20), Rdpack(v.2.4), cli(v.3.6.0), zlibbioc(v.1.44.0), sn(v.2.1.0), rstudioapi(v.0.14), miniUI(v.0.1.1.1), sp(v.1.6-0), bslib(v.0.4.2), mathjaxr(v.1.6-0), fastmatch(v.1.1-3), aod(v.1.3.2), treeio(v.1.22.0), shiny(v.1.7.4), xfun(v.0.37), multtest(v.2.54.0), pkgbuild(v.1.4.0), gson(v.0.0.9), cluster(v.2.1.4), caTools(v.1.18.2), tidygraph(v.1.2.2), KEGGREST(v.1.38.0), clusterGeneration(v.1.3.7), ggrepel(v.0.9.2), ape(v.5.6-2), listenv(v.0.9.0), Biostrings(v.2.66.0), png(v.0.1-8), withr(v.2.5.0), bitops(v.1.0-7), ggforce(v.0.4.1), plyr(v.1.8.8), pillar(v.1.8.1), gplots(v.3.1.3), cachem(v.1.0.6), GenomicFeatures(v.1.50.4), multcomp(v.1.4-20), fs(v.1.6.0), clusterProfiler(v.4.6.0), RUnit(v.0.4.32), vctrs(v.0.5.2), ellipsis(v.0.3.2), generics(v.0.1.3), devtools(v.2.4.5), metap(v.1.8), tools(v.4.2.0), remaCor(v.0.0.11), beeswarm(v.0.4.0), munsell(v.0.5.0), tweenr(v.2.0.2), fgsea(v.1.24.0), DelayedArray(v.0.24.0), fastmap(v.1.1.0), compiler(v.4.2.0), pkgload(v.1.3.2), abind(v.1.4-5), httpuv(v.1.6.8), rtracklayer(v.1.58.0), sessioninfo(v.1.2.2), plotly(v.4.10.1), GenomeInfoDbData(v.1.2.9), gridExtra(v.2.3), edgeR(v.3.40.2), lattice(v.0.20-45), deldir(v.1.0-6), mutoss(v.0.1-12), utf8(v.1.2.3), later(v.1.3.0), dplyr(v.1.1.0), BiocFileCache(v.2.6.0), jsonlite(v.1.8.4), scales(v.1.2.1), tidytree(v.0.4.2), pbapply(v.1.7-0), genefilter(v.1.80.3), lazyeval(v.0.2.2), promises(v.1.2.0.1), doParallel(v.1.0.17), R.utils(v.2.12.2), goftest(v.1.2-3), spatstat.utils(v.3.0-1), sandwich(v.3.0-2), rmarkdown(v.2.20), openxlsx(v.4.2.5.1), cowplot(v.1.1.1), Rtsne(v.0.16), pander(v.0.6.5), downloader(v.0.4), uwot(v.0.1.14), igraph(v.1.3.5), plotrix(v.3.8-2), numDeriv(v.2016.8-1.1), survival(v.3.5-0), yaml(v.2.3.7), htmltools(v.0.5.4), memoise(v.2.0.1), profvis(v.0.3.7), BiocIO(v.1.8.0), locfit(v.1.5-9.7), graphlayouts(v.0.8.4), viridisLite(v.0.4.1), digest(v.0.6.31), assertthat(v.0.2.1), RhpcBLASctl(v.0.21-247.1), mime(v.0.12), rappdirs(v.0.3.3), repr(v.1.1.6), RSQLite(v.2.2.20), yulab.utils(v.0.0.6), future.apply(v.1.10.0), remotes(v.2.4.2), data.table(v.1.14.6), urlchecker(v.1.0.1), blob(v.1.2.3), R.oo(v.1.25.0), labeling(v.0.4.2), splines(v.4.2.0), RCurl(v.1.98-1.10), broom(v.1.0.3), hms(v.1.1.2), colorspace(v.2.1-0), base64enc(v.0.1-3), mnormt(v.2.1.1), ggbeeswarm(v.0.7.1), aplot(v.0.1.9), ggrastr(v.1.0.1), sass(v.0.4.5), Rcpp(v.1.0.10), RANN(v.2.6.1), mvtnorm(v.1.1-3), enrichplot(v.1.18.3), fansi(v.1.0.4), tzdb(v.0.3.0), brio(v.1.1.3), parallelly(v.1.34.0), R6(v.2.5.1), grid(v.4.2.0), ggridges(v.0.5.4), lifecycle(v.1.0.3), TFisher(v.0.2.0), zip(v.2.2.2), curl(v.5.0.0), minqa(v.1.2.5), leiden(v.0.4.3), jquerylib(v.0.1.4), PROPER(v.1.30.0), Matrix(v.1.5-3), qvalue(v.2.30.0), TH.data(v.1.1-1), desc(v.1.4.2), RcppAnnoy(v.0.0.20), RColorBrewer(v.1.1-3), iterators(v.1.0.14), spatstat.explore(v.3.0-6), stringr(v.1.5.0), htmlwidgets(v.1.6.1), polyclip(v.1.10-4), biomaRt(v.2.54.0), shadowtext(v.0.1.2), gridGraphics(v.0.5-1), mgcv(v.1.8-41), globals(v.0.16.2), patchwork(v.1.1.2), spatstat.random(v.3.1-3), progressr(v.0.13.0), codetools(v.0.2-19), GO.db(v.3.16.0), gtools(v.3.9.4), prettyunits(v.1.1.1), dbplyr(v.2.3.0), R.methodsS3(v.1.8.2), gtable(v.0.3.1), DBI(v.1.1.3), highr(v.0.10), ggfun(v.0.0.9), tensor(v.1.5), httr(v.1.4.4), KernSmooth(v.2.23-20), vroom(v.1.6.1), stringi(v.1.7.12), progress(v.1.2.2), reshape2(v.1.4.4), farver(v.2.1.1), viridis(v.0.6.2), ggtree(v.3.6.2), xml2(v.1.3.3), boot(v.1.3-28.1), lme4(v.1.1-31), restfulr(v.0.0.15), readr(v.2.1.3), ggplotify(v.0.1.0), scattermore(v.0.8), bit(v.4.0.5), scatterpie(v.0.1.8), spatstat.data(v.3.0-0), ggraph(v.2.1.0), pkgconfig(v.2.0.3) and knitr(v.1.42)
message(paste0("This is hpgltools commit: ", get_git_commit()))
## If you wish to reproduce this exact build of hpgltools, invoke the following:
## > git clone http://github.com/abelew/hpgltools.git
## > git reset 140f18734ff1a7e49669ad91dd4f469f0560ca1a
## This is hpgltools commit: Sun Jan 29 15:36:20 2023 -0500: 140f18734ff1a7e49669ad91dd4f469f0560ca1a
<- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
this_save ##message(paste0("Saving to ", this_save))
##tmp <- sm(saveme(filename=this_save))