Skip to content

Commit f9a61f9

Browse files
darlanminussidarlanminussi
committed
adding new files
1 parent 9e41509 commit f9a61f9

12 files changed

+1753
-0
lines changed

‎05-knnsmoothing.Rmd

+14
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
# kNN smoothing
2+
3+
scDNA-Seq experiments are often sparsely sequenced and present inherent noise from low-coverage datasets. To mitigate this, CopyKit is capable of smoothing profiles based on the k-nearest neighbors of each single cell.
4+
5+
This way, CopyKit aggregates the genomic bin counts based on a cell k-nearest neighbors, followed by re-segmentation of copy number profiles.
6+
7+
In order to ensure optimal performane it is essential to carefully consider the number of neighbors (k) used in the smoothing process. We recommend to use conservative k values, which are below the number of cells that compose the smallest subclone, and to visually inspect and compare smoothed single cells to the original profiles
8+
9+
```{r knn_smoothing}
10+
tumor <- knnSmooth(tumor)
11+
```
12+
13+
This step is recommended after the filtering process. Though we have noticed that smoothing can also rescue cells with low-depth quality. If done so, additional inspection is recommended due to the possibility of increased noise.
14+

‎06-analysis_visualization.Rmd

+365
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,365 @@
1+
# Analysis and Visualization module
2+
3+
The analysis and visualization module from CopyKit work in synergy to help you analyze your single cell copy number results.
4+
5+
## plotMetrics()
6+
7+
The `plotMetrics()` function can plot any information store in `colData()`, which is passedto the `metric` argument. If a `label` argument is supplied, the points will be colored based on that information.
8+
```{r}
9+
plotMetrics(tumor, metric = c("overdispersion",
10+
"breakpoint_count",
11+
"reads_total",
12+
"reads_duplicates",
13+
"reads_assigned_bins",
14+
"percentage_duplicates"),
15+
label = "percentage_duplicates")
16+
```
17+
18+
## plotRatio()
19+
20+
Visualizing the seegmeentation to eensure it follows the ratios as expected is crucial to a copy number analysis.
21+
This helps to verify the accuracy of the segmentation an assess the quality of the data visually.
22+
23+
The `plotRatio()` function is a useful tool for this, offering two different modes. When the input is the CopyKit object, an interactive app will open, allowing you to choose which cell to visualize from the drop-down option menu.
24+
```{r plot_ratio_app, eval = FALSE}
25+
plotRatio(tumor)
26+
```
27+
<center>
28+
![Ratio Plot App](images/plot_ratio_shiny.png)
29+
</center>
30+
31+
Providing a sample name to the `plotRatio()` function with the argument `sample_name`, will display the plot only for the selected cell.
32+
33+
```{r plot_ratio_sample, warning = FALSE}
34+
plotRatio(tumor, "PMTC6LiverC117AL4L5S1_S885_L003_R1_001")
35+
```
36+
37+
## runUmap()
38+
39+
The `runUmap()` function generates a [UMAP](https://umap-learn.readthedocs.io/en/latest/) embedding, which is stored in the `reducedDim` slot.
40+
41+
`runUmap()` is an important pre-processing step for the `findClusters()` feature.
42+
43+
```{r run_umap}
44+
tumor <- runUmap(tumor)
45+
```
46+
47+
You can pass additional arguments to control UMAP parameters using the `'...'` argument. The full list of additional arguments that can be passed on to `uwot::umap` with the '...' argument can be found in the [uwot manual](https://cran.r-project.org/web/packages/uwot/uwot.pdf) and information on their influence on clustering can be seen in the [UMAP webpage](https://umap-learn.readthedocs.io/en/latest/clustering.html)
48+
49+
## Clustering
50+
51+
The `findClusters()` function performs clustering using the UMAP embedding generated by `runUmap()`. For detecting subclones, CopyKit uses
52+
different clustering algorithms such as:
53+
54+
1) [hdbscan](https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html) (recommended)
55+
56+
2) [Leiden](https://www.nature.com/articles/s41598-019-41695-z)
57+
58+
3) [Louvain](https://arxiv.org/abs/0803.0476)
59+
60+
The hdbscan method is recommended and has been previously successfully applied in the work from Laks *et al.* and Minussi *et al.*.
61+
62+
### findSugestedK()
63+
64+
The `findSuggestedK()` is a helper function that provides guidance to choose the `k` parameter for clustering algorithms. The function
65+
`findSuggestedK` bootstraps clustering over a range of k values, and returns the value that maximizes the jaccard similarity, with *median* as the default metric (argument `metric`).
66+
67+
While `findSuggestedK` does not guarantee optimal clustering. It provides a guide to select k values.
68+
69+
```{r find_suggested_k}
70+
tumor <- findSuggestedK(tumor)
71+
```
72+
73+
The `plotSuggestedK()` function can be used to inspect the results of `findSuggestedK()`. It can plot a boxplot, heatmap (tile), or scatterplot
74+
75+
```{r plot_suggested_k_boxplot}
76+
plotSuggestedK(tumor)
77+
```
78+
79+
If the argument `geom` is set to *tile*, `plotSuggestedK()` plots a heatmap where each row is a detected subclone, each column is a k assessed during the grid search and the color represents the jaccard similarity for a given clone.
80+
81+
The suggested value can be accessed from the metadata:
82+
```{r find_sug_meta}
83+
S4Vectors::metadata(tumor)$suggestedK
84+
```
85+
86+
### findClusters()
87+
88+
To run `findClusters()`, pass the CopyKit object as its input. If `findSuggestedK()` was run beforehand, `findClusters()` will automatically use the stored value from `findSuggestedK()` as the argument for k_subclones, unless otherwise specified.
89+
90+
91+
```{r findCL_findSuggested}
92+
tumor <- findClusters(tumor)
93+
```
94+
95+
If hdbscan is used for clustering, singleetons [outliers](https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html) may bee idenetifieed and added to subgroup `c0`. Copykit will notify you of any cells classified as c0 after running findClusters(). To remove outliers, subset the CopyKit object with:
96+
97+
```{r filter_outliers}
98+
tumor <- tumor[,colData(tumor)$subclones != 'c0']
99+
```
100+
101+
### plotUmap()
102+
103+
`plotUmap()` can be used to visualize the reduced dimensional embedding.
104+
`plotUmap` can be colored by any element of the colData with the argument 'label'.
105+
106+
```{r}
107+
plotUmap(tumor)
108+
```
109+
110+
Cells can be colored with the argument 'label' by any element of `colData()`:
111+
112+
```{r plot_umap_subclones}
113+
plotUmap(tumor, label = 'subclones')
114+
```
115+
116+
## runPhylo()
117+
118+
To run a phylogenetic analysis of cells' copy number profiles, use the function `runPhylo()`.
119+
Available methods are [Neighbor Joining](https://pubmed.ncbi.nlm.nih.gov/3447015/) and [Balanced Minimum evolution](https://pubmed.ncbi.nlm.nih.gov/8412650/).
120+
121+
The resulting tree is stored within the CopyKit object in the phylo slot:
122+
```{r run_phylo}
123+
tumor <- runPhylo(tumor, metric = 'manhattan')
124+
```
125+
126+
### plotPhylo()
127+
128+
To visualize the resulting phylogenetic trees from `runPhylo()` we can use the `plotPhylo()` function.
129+
130+
```{r plot_phylo, fig.height=8}
131+
plotPhylo(tumor)
132+
```
133+
134+
`plotPhylo()` can use any element of the colData to color the leaves of the tree.
135+
136+
```{r plot_phylo_label, fig.height=8}
137+
plotPhylo(tumor, label = 'subclones')
138+
```
139+
140+
## calcInteger()
141+
142+
The true underlying copy number states for a given region of the genome are given in integer states.
143+
We can use the segment ratio means to make an inference of the copy number state for a given genomic region.
144+
Keep in mind that those are inferences and, as such, are subjected to errors of different sources, like erroneous inference or rounding errors.
145+
146+
Segment ratio means can be scaled to integer values with `calcInteger()`.
147+
148+
CopyKit supports different methods of calculating integer copy number profiles.
149+
150+
To calculate computational ploidies, a method is in development and should be available soon within copykit.
151+
152+
By setting the argument `method` to *fixed*, a fixed value of ploidy (generally determined using Flow Cytometry) will scale all cells.
153+
154+
```{r calc_integer_fixed}
155+
tumor <- calcInteger(tumor, method = 'fixed', ploidy_value = 4.3)
156+
```
157+
158+
The integer values are stored in the slot `integer` that can be accessed with the function `assay()`.
159+
160+
If the assay _integer_ exists, `plotRatio()` will use it to plot the integer copy number values as a secondary axis.
161+
```{r plot_ratio_integer, warning=FALSE}
162+
plotRatio(tumor, "PMTC6LiverC117AL4L5S1_S885_L003_R1_001")
163+
```
164+
165+
## plotHeatmap()
166+
167+
To visualize copy number profiles with a heatmap we can use `plotHeatmap()`.
168+
169+
The heatmap can be annotated with elements of colData.
170+
171+
To order subclones, one option is to calculate a consensus phylogeny, explained in later sections:
172+
```{r}
173+
tumor <- calcConsensus(tumor)
174+
tumor <- runConsensusPhylo(tumor)
175+
```
176+
177+
To plot the heatmap.
178+
```{r plot_heatmap_subclones, fig.height=8}
179+
plotHeatmap(tumor, label = 'subclones', order_cells = 'consensus_tree')
180+
```
181+
182+
Genes can be annotated to a heatmap in `plotHeatmap()` with the _genes_ argument.
183+
The genes argument is a vector of gene HUGO symbols:
184+
```{r plot_heatmap_subclones_genes, fig.height=8}
185+
plotHeatmap(tumor, label = 'subclones', genes = c("TP53", "BRAF", "MYC"), order_cells = 'consensus_tree')
186+
```
187+
188+
To plot integer copy number heatmaps pass the argument `assay = 'integer'`, importantly the integer matrix must be in the `'integer'`slot.
189+
```{r plot_ht_int, fig.height=8}
190+
plotHeatmap(tumor, assay = 'integer', order_cells = 'consensus_tree')
191+
```
192+
193+
New information can be added to `colData()` and used in conjunction with the plotting functions.
194+
195+
The example dataset has macro-spatial information.
196+
The information is encoded in the sample name by the letter L followed by a number.
197+
198+
We can extract that information and add it as an extra column to the metadata:
199+
200+
```{r add_spatial_info}
201+
colData(tumor)$spatial_info <- stringr::str_extract(colData(tumor)$sample, "L[0-9]")
202+
```
203+
204+
Once the information has been added, we can use it to color the umap by their spatial information:
205+
206+
```{r umap_spatial}
207+
plotUmap(tumor, label = 'spatial_info')
208+
```
209+
210+
It is also possible to annotate the heatmap with that information:
211+
The 'label' argument for `plotHeatmap()` can add as many annotations as specified by the user as long as they are elements in `colData()` of the CopyKit object.
212+
```{r heat_spatial, fig.height=8}
213+
plotHeatmap(tumor, label = c("spatial_info", "subclones"), order_cells = 'consensus_tree')
214+
```
215+
216+
## plotFreq()
217+
218+
Computing the frequencies of genomic gain or losses across the genome can be useful to visualize differences between groups. This can be done with `plotFreq()`.
219+
For every region of the genome, `plotFreq()` will calculate the frequency of gain or losses according to a threshold across all samples.
220+
The thresholds are controlled with the arguments `low_threshold` (values below will be counted as genomic losses) and `high_threshold` (values above will be counted as genomic gains).
221+
222+
Ideally thresholds will be set according to the ploidy of the sample.
223+
224+
The argument `assay` can be provided to pass on the `integer` assay instead of the `segment_ratios` assay (adjust thresholds accordingly).
225+
226+
Two geoms are available `area` or `line`.
227+
```{r plot_freq_sample}
228+
plotFreq(tumor)
229+
```
230+
231+
Group information can be set to any categorical element of the `colData()` and is provided with the argument `group`.
232+
233+
```{r plot_freq_label, fig.height=6}
234+
plotFreq(tumor,
235+
group = 'subclones')
236+
```
237+
238+
## calcConsensus()
239+
240+
Consensus sequences can help visualize the different segments across subclones. To calculate consensus matrices we can use `calcConsensus()`.
241+
```{r calc_consensus}
242+
tumor <- calcConsensus(tumor)
243+
```
244+
245+
`plotHeatmap()` can plot a consensus heatmap:
246+
```{r plot_heat_consensus}
247+
plotHeatmap(tumor, consensus = TRUE, label = 'subclones')
248+
```
249+
250+
`plotHeatmap()` can annotate the consensus heatmap with information from the metadata as long as label is the same as the information used to build the consensus matrix:
251+
```{r plot_heat_freq_bar}
252+
plotHeatmap(tumor, consensus = TRUE, label = 'subclones', group = 'spatial_info')
253+
```
254+
255+
By default `calcConsensus()` uses the subclones information to calculate a consensus for each subclone.
256+
Any element of the `colData()` can be used to calculate the consensus.
257+
258+
**Note:** Consensus matrices can be calculated from the integer assay. Importantly, the integer matrix must be in the `assay(tumor, 'integer')` slot. Check `calcInteger()` for more info.
259+
```{r eval = FALSE}
260+
tumor <- calcConsensus(tumor, consensus_by = 'subclones', assay = 'integer')
261+
```
262+
263+
## plotConsensusLine()
264+
265+
To compare the differences among subclones, `plotConsensusLine()` opens an interactive app where the consensus sequences are plotted as lines.
266+
267+
```{r calc_consensus_again, echo = FALSE}
268+
tumor <- calcConsensus(tumor)
269+
```
270+
271+
```{r plot_cons_line, eval=FALSE}
272+
plotConsensusLine(tumor)
273+
```
274+
275+
<center>
276+
![Plot Consensus Line](images/plot_consensus_shiny.png)
277+
</center>
278+
279+
## plotGeneCopy()
280+
281+
To check copy number states across of genes we can use `plotGeneCopy()`. Two different geoms: “swarm” (default) or “violin” can be applied.
282+
283+
As with other plotting functions, points can be colored with the argument 'label'.
284+
285+
```{r}
286+
plotGeneCopy(tumor, genes = c("CDKN2A",
287+
"FGFR1",
288+
"TP53",
289+
"PTEN",
290+
"MYC",
291+
"CDKN1A",
292+
"MDM2",
293+
"AURKA",
294+
"PIK3CA",
295+
"CCND1",
296+
"KRAS"),
297+
label = 'spatial_info')
298+
```
299+
300+
A positional dodge can be added to facilitate the visualization across groups:
301+
```{r}
302+
plotGeneCopy(tumor,
303+
genes = c("CDKN2A",
304+
"FGFR1",
305+
"TP53",
306+
"PTEN",
307+
"MYC",
308+
"CDKN1A",
309+
"MDM2",
310+
"AURKA",
311+
"PIK3CA",
312+
"CCND1",
313+
"KRAS"),
314+
label = 'subclones',
315+
dodge.width = 0.8)
316+
```
317+
318+
319+
A barplot geom is also provided to visualize the integer data as a frequency barplot for each gene:
320+
```{r plot_gene_copy_integer}
321+
plotGeneCopy(tumor, genes = c("CDKN2A",
322+
"FGFR1",
323+
"TP53",
324+
"PTEN",
325+
"MYC",
326+
"CDKN1A",
327+
"MDM2",
328+
"AURKA",
329+
"PIK3CA",
330+
"CCND1",
331+
"KRAS"),
332+
geom = 'barplot',
333+
assay = 'integer')
334+
```
335+
336+
## plotAlluvial()
337+
338+
To visualize frequencies across elements of the metadata we can use `plotAlluvial()`
339+
```{r plot_alluvial}
340+
plotAlluvial(tumor, label = c("subclones", "spatial_info"))
341+
```
342+
343+
## runPca()
344+
345+
Principal Component Analysis can be run with `runPca()`
346+
```{r run_pca}
347+
tumor <- runPca(tumor)
348+
```
349+
350+
The results are stored in the `reducedDim(tumor, 'PCA')` slot.
351+
352+
Principal Components can be used as an alternative pre-processing step to the `findClusters()` function by changing `findClusters()` _embedding_ argument to 'PCA'. The number of principal components to be used for clustering can also be set within `findClusters() with the argument ncomponents. At this time better results are achieved by using the default UMAP embedding and it is the recommend pre-processing step. However PCA remains as a viable alternative
353+
354+
### plotScree
355+
To determine the number of useful principal components a scree plot can be shown with:
356+
```{r plot_scree}
357+
plotScree(tumor)
358+
```
359+
360+
### plotPCA
361+
To visualize the first 2 principal components you can use `plotPca()`.
362+
```{r}
363+
plotPca(tumor, label = 'spatial_info')
364+
```
365+

0 commit comments

Comments
 (0)