-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy path3_RNAFucciPseudotime.py
214 lines (192 loc) · 9.4 KB
/
3_RNAFucciPseudotime.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
# -*- coding: utf-8 -*-
"""
Analysis of transcript abundance in individual cells over cell division time.
- Cell division time is measured with FUCCI markers and modeled in log-log space using polar coordinates.
- RNA abundance was measured with single-cell RNA sequencing.
- FUCCI marker intensities are measured for each individual cell with fluorescence assisted cell sorting (FACS)
@author: Anthony J. Cesnik, [email protected]
"""
from SingleCellProteogenomics import (FucciPseudotime,
RNACellCycleDependence,
RNADataPreparation,
RNAVelocity,
utils)
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
import shutil
import os
import argparse
if __name__ == "__main__":
description = "Single cell proteogenomic analysis script -- protein cell cycle dependence"
parser = argparse.ArgumentParser(description=description)
parser.add_argument(
"--quicker", action="store_true", help="Skip some plotting and analyses, namely for HPA releases."
)
args = parser.parse_args()
doAllPlotsAndAnalyses = not args.quicker
# Make PDF text readable
plt.rcParams["pdf.fonttype"] = 42
plt.rcParams["ps.fonttype"] = 42
plt.rcParams["savefig.dpi"] = 300
bioccd = np.genfromtxt(
"input/ProteinData/BiologicallyDefinedCCD.txt", dtype="str"
) # from mitotic structures
wp_ensg = np.load("output/pickles/wp_ensg.npy", allow_pickle=True)
ccd_comp = np.load("output/pickles/ccd_comp.npy", allow_pickle=True)
nonccd_comp = np.load("output/pickles/nonccd_comp.npy", allow_pickle=True)
u_plates = ["355","356","357"]
#%% Convert FACS intensities for FUCCI markers to pseudotime using the same polar coordinate methods as for protein
# Idea: Use the polar coordinate pseudotime calculations to calculate the pseudotime for each cell
# Execution: Adapt Devin's code for the cells sorted for RNA-Seq
# Output: Make log-log fucci intensity plots for the cells analyzed by RNA-Seq; Plot of all fucci pseudotimes; table of pseudotimes for each cell
adata, phases_filt = RNADataPreparation.read_counts_and_phases(
"Counts", False, "protein_coding", u_plates
)
adata = RNADataPreparation.zero_center_fucci(adata)
FucciPseudotime.pseudotime_rna(adata)
#%% Single cell RNA-Seq data preparation and general analysis
if doAllPlotsAndAnalyses:
RNADataPreparation.general_plots(u_plates)
RNADataPreparation.analyze_noncycling_cells(u_plates)
adata, phasesfilt = RNADataPreparation.qc_filtering(
adata, do_log_normalize=True, do_remove_blob=False
)
adata = RNADataPreparation.zero_center_fucci(adata)
RNADataPreparation.plot_pca_for_batch_effect_analysis(adata, "BeforeRemovingNoncycling")
#%% Idea: Similar to mock-bulk analysis for proteins, we can evaluate each gene bundled by phase across cells
# Execution: Make boxplots of RNA expression by phase
# Output: boxplots for each gene
valuetype, use_spikeins, biotype_to_use = "Tpms", False, "protein_coding"
adata, phases = RNADataPreparation.read_counts_and_phases(
valuetype, use_spikeins, biotype_to_use, u_plates
)
adata, phasesfilt = RNADataPreparation.qc_filtering(
adata, do_log_normalize=True, do_remove_blob=True
)
adata = RNADataPreparation.zero_center_fucci(adata)
if doAllPlotsAndAnalyses:
RNADataPreparation.plot_markers_vs_reads(adata)
RNADataPreparation.plot_pca_for_batch_effect_analysis(adata, "AfterRemovingNoncycling")
g1 = adata.obs["phase"] == "G1"
s = adata.obs["phase"] == "S-ph"
g2 = adata.obs["phase"] == "G2M"
do_make_boxplots = False
if do_make_boxplots:
for iii, ensg in enumerate(adata.var_names):
maxtpm = np.max(
np.concatenate((adata.X[g1, iii], adata.X[s, iii], adata.X[g2, iii]))
)
RNACellCycleDependence.boxplot_result(
adata.X[g1, iii] / maxtpm,
adata.X[s, iii] / maxtpm,
adata.X[g2, iii] / maxtpm,
"figures/RNABoxplotByPhase",
ensg,
)
#%% Idea: Display general RNA expression patterns in single cells using UMAP dimensionality reduction, and display with FUCCI pseudotime overlayed
if doAllPlotsAndAnalyses:
FucciPseudotime.pseudotime_umap(adata) # Generate a UMAP with the pseudotime overlayed
# We can also show that the cycle pattern remains when the curated CCD genes or CCD proteins are removed,
# demonstrating that there's still valuable information about cell cycling beyond what was called CCD
RNADataPreparation.demonstrate_umap_cycle_without_ccd(adata)
# Read in the curated CCD genes / CCD proteins from the present work / Non-CCD genes from the present work; filter for genes that weren't filtered in QC of RNA-Seq
ccd_regev_filtered, ccd_filtered, nonccd_filtered = utils.ccd_gene_lists(adata)
adata_ccdprotein, adata_nonccdprotein, adata_regevccdgenes = RNADataPreparation.is_ccd(
adata, wp_ensg, ccd_comp, nonccd_comp, bioccd, ccd_regev_filtered
)
# Generate plots with expression of genes overlayed
expression_data = adata.X
normalized_exp_data = (expression_data.T / np.max(expression_data, axis=0)[:, None]).T
# Log-log FUCCI plot with RNA expression overlayed
if doAllPlotsAndAnalyses:
RNACellCycleDependence.plot_expression_facs(
adata,
wp_ensg[np.isin(wp_ensg, adata.var_names)],
normalized_exp_data,
"figures/GeneExpressionFucci",
)
# Cluster the expression into phases and analyze it that way
bulk_phase_tests = RNACellCycleDependence.analyze_ccd_variation_by_phase_rna(
adata, normalized_exp_data, biotype_to_use
)
# RNACellCycleDependence.plot_expression_boxplots(adata, wp_ensg[np.isin(wp_ensg, adata.var_names)], bulk_phase_tests, "figures/GeneExpressionBoxplots")
# RNACellCycleDependence.plot_expression_umap(adata, wp_ensg[np.isin(wp_ensg, adata.var_names)], "figures/GeneExpressionUmap")
#%% Moving average calculations and randomization analysis for RNA
rna_ccd_analysis_results = RNACellCycleDependence.analyze_ccd_variation_by_mvavg_rna(
adata,
wp_ensg,
ccd_comp,
bioccd,
adata_nonccdprotein,
adata_regevccdgenes,
biotype_to_use,
make_mvavg_plots_isoforms=doAllPlotsAndAnalyses,
)
percent_ccd_variance, total_gini, mean_diff_from_rng, pass_meandiff, eq_percvar_adj, fucci_time_inds, norm_exp_sort, moving_averages, mvavg_xvals, perms, ccdtranscript, ccdprotein, mvpercs = (
rna_ccd_analysis_results
)
RNACellCycleDependence.make_plotting_dataframe(
adata,
ccdtranscript,
wp_ensg,
bioccd,
norm_exp_sort[np.argsort(fucci_time_inds), :],
mvavg_xvals,
moving_averages,
mvpercs,
)
if doAllPlotsAndAnalyses:
RNACellCycleDependence.plot_umap_ccd_cutoffs(adata, mean_diff_from_rng)
RNACellCycleDependence.figures_ccd_analysis_rna(
adata,
percent_ccd_variance,
mean_diff_from_rng,
pass_meandiff,
eq_percvar_adj,
wp_ensg,
ccd_comp,
ccd_regev_filtered,
)
RNACellCycleDependence.plot_overall_and_ccd_variances(
adata,
biotype_to_use,
total_gini,
percent_ccd_variance,
pass_meandiff,
adata_ccdprotein,
adata_nonccdprotein,
adata_regevccdgenes,
)
RNACellCycleDependence.compare_to_lasso_analysis(adata, ccdtranscript)
RNACellCycleDependence.analyze_cnv_calls(adata, ccdtranscript)
#%% Moving average calculations and randomization analysis for the spike-in internal controls
if doAllPlotsAndAnalyses:
adata_spikeins, phases_spikeins = RNADataPreparation.read_counts_and_phases(
valuetype, use_spike_ins=True, biotype_to_use="", u_plates=u_plates
)
sc.pp.filter_genes(adata_spikeins, min_cells=100)
print(f"data shape after filtering: {adata_spikeins.X.shape}")
RNACellCycleDependence.ccd_analysis_of_spikeins(adata_spikeins, perms)
#%% Analyze RNA velocity
valuetype, use_spikeins, biotype_to_use = "Tpms", False, "protein_coding"
adata, phases = RNADataPreparation.read_counts_and_phases(
valuetype, use_spikeins, biotype_to_use, u_plates, load_velocities=True
)
adata_blobless, phasesfilt = RNADataPreparation.qc_filtering(
adata, do_log_normalize=True, do_remove_blob=False
)
RNAVelocity.analyze_rna_velocity(adata_blobless, mean_diff_from_rng, doAllPlotsAndAnalyses)
#%% Analyze isoforms
if doAllPlotsAndAnalyses:
adata_isoform, ccdtranscript_isoform = RNACellCycleDependence.analyze_isoforms(
adata, ccdtranscript, wp_ensg, ccd_comp, nonccd_comp, u_plates, make_mvavg_plots_isoforms=doAllPlotsAndAnalyses
)
RNACellCycleDependence.compare_genes_to_isoforms(
adata,
ccdprotein,
ccdtranscript,
adata_nonccdprotein,
adata_isoform,
ccdtranscript_isoform,
)