infercnv_run
Create the InferCNV Object
Reading in the raw counts matrix and meta data, populating the infercnv object
infercnv_obj = CreateInfercnvObject(
raw_counts_matrix="oligodendroglioma_expression_downsampled.counts.matrix",
annotations_file="oligodendroglioma_annotations_downsampled.txt",
delim="\t",
gene_order_file="gencode_downsampled.txt",
ref_group_names=c("Microglia/Macrophage","Oligodendrocytes (non-malignant)"))## INFO [2018-11-07 15:21:50] ::order_reduce:Start.
## INFO [2018-11-07 15:21:50] .order_reduce(): expr and order match.
## INFO [2018-11-07 15:21:50] ::process_data:order_reduce:Reduction from positional data, new dimensions (r,c) = 10338,184 Total=18322440.6799817 Min=0 Max=34215.
## INFO [2018-11-07 15:21:50] validating infercnv_objFiltering genes
Removing those genes that are very lowly expressed or present in very few cells
# filter out low expressed genes
cutoff=1
infercnv_obj <- require_above_min_mean_expr_cutoff(infercnv_obj, cutoff)## INFO [2018-11-07 15:21:50] ::above_min_mean_expr_cutoff:Start
## INFO [2018-11-07 15:21:50] Removing 1510 genes from matrix as below mean expr threshold: 1
## INFO [2018-11-07 15:21:50] validating infercnv_obj
## INFO [2018-11-07 15:21:50] There are 8828 genes and 184 cells remaining in the expr matrix.# filter out bad cells
min_cells_per_gene=3
infercnv_obj <- require_above_min_cells_ref(infercnv_obj, min_cells_per_gene=min_cells_per_gene)## INFO [2018-11-07 15:21:50] no genes removed due to min cells/gene filter## for safe keeping
infercnv_orig_filtered = infercnv_obj
#plot_mean_chr_expr_lineplot(infercnv_obj)
save('infercnv_obj', file = 'infercnv_obj.orig_filtered')Normalize each cell’s counts for sequencing depth
Perform a total sum normalization. Generates counts-per-million or counts-per-100k, depending on the overall sequencing depth.
infercnv_obj <- infercnv:::normalize_counts_by_seq_depth(infercnv_obj)## INFO [2018-11-07 15:21:51] Computed total sum normalization factor as: 100000.000000Spike in artificial variation for tracking purposes
Add ~0x and 2x variation to an artificial spike-in data set based on the normal cells so we can track and later scale residual expression data to this level of variation.
infercnv_obj <- spike_in_variation_chrs(infercnv_obj)## INFO [2018-11-07 15:21:51] Selecting longest chrs for adding spike: chr1,chr2
## INFO [2018-11-07 15:21:51] processing group: malignant_MGH36
## INFO [2018-11-07 15:21:52] processing group: malignant_MGH53
## INFO [2018-11-07 15:21:52] processing group: malignant_93
## INFO [2018-11-07 15:21:52] processing group: malignant_97
## INFO [2018-11-07 15:21:53] processing group: Microglia/Macrophage
## INFO [2018-11-07 15:21:53] processing group: Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:21:53] processing group: malignant_MGH36
## INFO [2018-11-07 15:21:53] processing group: malignant_MGH53
## INFO [2018-11-07 15:21:53] processing group: malignant_93
## INFO [2018-11-07 15:21:53] processing group: malignant_97
## INFO [2018-11-07 15:21:54] processing group: Microglia/Macrophage
## INFO [2018-11-07 15:21:54] processing group: Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:22:01] validating infercnv_objperform Anscombe normalization
<<<<<<< HEAD Suggested for removing noisy variation at low counts ======= Useful noise reduction method.
See: https://en.wikipedia.org/wiki/Anscombe_transform >>>>>>> 29a0b973d2701fe5ea2834efcd6a82dd542e0308
infercnv_obj <- infercnv:::anscombe_transform(infercnv_obj)
save('infercnv_obj', file='infercnv_obj.anscombe')log transform the normalized counts:
infercnv_obj <- log2xplus1(infercnv_obj)## INFO [2018-11-07 15:22:02] transforming log2xplus1()save('infercnv_obj', file='infercnv_obj.log_transformed')Apply maximum bounds to the expression data to reduce outlier effects
Here we define a threshold by taking the mean of the bounds of expression data across all cells. This is then use to define a cap for the bounds of all data.
threshold = mean(abs(get_average_bounds(infercnv_obj)))
infercnv_obj <- apply_max_threshold_bounds(infercnv_obj, threshold=threshold)## INFO [2018-11-07 15:22:04] ::process_data:setting max centered expr, threshold set to: +/-: 4.23476201265372Initial view, before inferCNV operations:
plot_cnv(infercnv_obj,
output_filename='infercnv.logtransf',
x.range="auto",
title = "Before InferCNV (filtered & log2 transformed)",
color_safe_pal = FALSE,
x.center = mean(infercnv_obj@expr.data))knitr::include_graphics("infercnv.logtransf.png")
perform smoothing across chromosomes
The expression values are
infercnv_obj = smooth_by_chromosome(infercnv_obj, window_length=101, smooth_ends=TRUE)## INFO [2018-11-07 15:22:35] ::smooth_window:Start.
## INFO [2018-11-07 15:22:35] ::smooth_window:Start.
## INFO [2018-11-07 15:22:36] ::smooth_window:Start.
## INFO [2018-11-07 15:22:36] ::smooth_window:Start.
## INFO [2018-11-07 15:22:36] ::smooth_window:Start.
## INFO [2018-11-07 15:22:36] ::smooth_window:Start.
## INFO [2018-11-07 15:22:37] ::smooth_window:Start.
## INFO [2018-11-07 15:22:37] ::smooth_window:Start.
## INFO [2018-11-07 15:22:37] ::smooth_window:Start.
## INFO [2018-11-07 15:22:38] ::smooth_window:Start.
## INFO [2018-11-07 15:22:38] ::smooth_window:Start.
## INFO [2018-11-07 15:22:38] ::smooth_window:Start.
## INFO [2018-11-07 15:22:38] ::smooth_window:Start.
## INFO [2018-11-07 15:22:39] ::smooth_window:Start.
## INFO [2018-11-07 15:22:39] ::smooth_window:Start.
## INFO [2018-11-07 15:22:39] ::smooth_window:Start.
## INFO [2018-11-07 15:22:40] ::smooth_window:Start.
## INFO [2018-11-07 15:22:40] ::smooth_window:Start.
## INFO [2018-11-07 15:22:40] ::smooth_window:Start.
## INFO [2018-11-07 15:22:41] ::smooth_window:Start.
## INFO [2018-11-07 15:22:41] ::smooth_window:Start.
## INFO [2018-11-07 15:22:41] ::smooth_window:Start.
## INFO [2018-11-07 15:22:41] ::smooth_window:Start.
## INFO [2018-11-07 15:22:42] ::smooth_window:Start.save('infercnv_obj', file='infercnv_obj.smooth_by_chr')
# re-center each cell
infercnv_obj <- center_cell_expr_across_chromosome(infercnv_obj, method = "median")## INFO [2018-11-07 15:22:43] ::center_smooth across chromosomes per cellsave('infercnv_obj', file='infercnv_obj.cells_recentered')plot_cnv(infercnv_obj,
output_filename='infercnv.chr_smoothed',
x.range="auto",
title = "chr smoothed and cells re-centered",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.chr_smoothed.png")
subtract the reference values from observations, now have log(fold change) values
infercnv_obj <- subtract_ref_expr_from_obs(infercnv_obj, inv_log=TRUE)## INFO [2018-11-07 15:23:17] ::subtract_ref_expr_from_obs:Start
## INFO [2018-11-07 15:23:17] subtracting mean(normal) per gene per cell across all datasave('infercnv_obj', file='infercnv_obj.ref_subtracted')plot_cnv(infercnv_obj,
output_filename='infercnv.ref_subtracted',
x.range="auto",
title="ref subtracted",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.ref_subtracted.png")
invert log values
Converting the log(FC) values to regular fold change values, centered at 1 (no fold change)
This is important because we want (1/2)x to be symmetrical to 1.5x, representing loss/gain of one chromosome region.
infercnv_obj <- invert_log2(infercnv_obj)## INFO [2018-11-07 15:23:50] invert_log2(), computing 2^xsave('infercnv_obj', file='infercnv_obj.inverted_log')plot_cnv(infercnv_obj,
output_filename='infercnv.inverted',
color_safe_pal = FALSE,
x.range="auto",
x.center=1,
title = "inverted log FC to FC")knitr::include_graphics("infercnv.inverted.png")
Removing noise
infercnv_obj <- clear_noise_via_ref_mean_sd(infercnv_obj, sd_amplifier = 1.0)## INFO [2018-11-07 15:24:26] :: **** clear_noise_via_ref_quantiles **** : removing noise between bounds: 0.95005640663898 - 1.0511317026573save('infercnv_obj', file='infercnv_obj.denoised')plot_cnv(infercnv_obj,
output_filename='infercnv.denoised',
x.range="auto",
x.center=1,
title="denoised",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.denoised.png")
Remove outlier data points
This generally improves on the visualization
infercnv_obj = remove_outliers_norm(infercnv_obj)## INFO [2018-11-07 15:24:59] ::remove_outlier_norm:Start out_method: average_bound lower_bound: NA upper_bound: NA
## INFO [2018-11-07 15:24:59] ::remove_outlier_norm:Start out_method: average_bound lower_bound: NA upper_bound: NA
## INFO [2018-11-07 15:24:59] ::remove_outlier_norm using method: average_bound for defining outliers.
## INFO [2018-11-07 15:25:00] outlier bounds defined between: 0.655402 - 1.36342save('infercnv_obj', file="infercnv_obj.outliers_removed")plot_cnv(infercnv_obj,
output_filename='infercnv.outliers_removed',
color_safe_pal = FALSE,
x.range="auto",
x.center=1,
title = "outliers removed")knitr::include_graphics("infercnv.outliers_removed.png")
Scale residual expression values according to the Spike-in
Perform rescaling of the data according to the spike-in w/ preset variation levels. Then, remove the spike-in data.
# rescale
infercnv_obj <- scale_cnv_by_spike(infercnv_obj)
# remove the spike-in
infercnv_obj <- remove_spike(infercnv_obj)## INFO [2018-11-07 15:25:38] Removing spikeMask out those genes that are not signficantly different from the normal cells
Runs a Wilcoxon rank test comparing tumor/normal for each patient and normal sample, and masks out those genes that are not significantly DE.
infercnv_obj <- infercnv:::mask_non_DE_genes_basic(infercnv_obj, test.use = ‘t’, center_val=1)
save(‘infercnv_obj’, file=“infercnv_obj.non_DE_masked”)
## INFO [2018-11-07 15:25:38] Finding DE genes between malignant_MGH36 and Microglia/Macrophage
## INFO [2018-11-07 15:25:42] Found 5010 genes / 8828 total as DE
## INFO [2018-11-07 15:25:42] Finding DE genes between malignant_MGH36 and Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:25:47] Found 4272 genes / 8828 total as DE
## INFO [2018-11-07 15:25:47] Finding DE genes between malignant_MGH53 and Microglia/Macrophage
## INFO [2018-11-07 15:25:51] Found 3943 genes / 8828 total as DE
## INFO [2018-11-07 15:25:51] Finding DE genes between malignant_MGH53 and Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:25:55] Found 3185 genes / 8828 total as DE
## INFO [2018-11-07 15:25:55] Finding DE genes between malignant_93 and Microglia/Macrophage
## INFO [2018-11-07 15:25:59] Found 4657 genes / 8828 total as DE
## INFO [2018-11-07 15:25:59] Finding DE genes between malignant_93 and Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:26:03] Found 3888 genes / 8828 total as DE
## INFO [2018-11-07 15:26:03] Finding DE genes between malignant_97 and Microglia/Macrophage
## INFO [2018-11-07 15:26:07] Found 4824 genes / 8828 total as DE
## INFO [2018-11-07 15:26:07] Finding DE genes between malignant_97 and Oligodendrocytes (non-malignant)
## INFO [2018-11-07 15:26:11] Found 4471 genes / 8828 total as DEplot_cnv(infercnv_obj,
output_filename='infercnv.non-DE-genes-masked',
color_safe_pal = FALSE,
x.range=c(0,2), # want 0-2 post scaling by the spike-in
x.center=1,
title = "non-DE-genes-masked")knitr::include_graphics("infercnv.non-DE-genes-masked.png")
Brighten it up by changing the scale threshold to our liking:
plot_cnv(infercnv_obj,
output_filename='infercnv.finalized_view',
color_safe_pal = FALSE,
x.range=c(0.7, 1.3),
x.center=1,
title = "InferCNV")## INFO [2018-11-07 15:26:31] ::plot_cnv:Start
## INFO [2018-11-07 15:26:31] ::plot_cnv:Current data dimensions (r,c)=8828,184 Total=1617829.28244849 Min=0 Max=2.05658892668798.
## INFO [2018-11-07 15:26:31] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-11-07 15:26:31] plot_cnv_observation:Start
## INFO [2018-11-07 15:26:31] Observation data size: Cells= 142 Genes= 8828
## INFO [2018-11-07 15:26:31] clustering observations via method: ward.D
## INFO [2018-11-07 15:26:31] Number of genes in group(1) is 33
## INFO [2018-11-07 15:26:31] group size being clustered: 33,8828
## INFO [2018-11-07 15:26:31] Number of genes in group(2) is 34
## INFO [2018-11-07 15:26:31] group size being clustered: 34,8828
## INFO [2018-11-07 15:26:31] Number of genes in group(3) is 40
## INFO [2018-11-07 15:26:31] group size being clustered: 40,8828
## INFO [2018-11-07 15:26:31] Number of genes in group(4) is 35
## INFO [2018-11-07 15:26:31] group size being clustered: 35,8828
## INFO [2018-11-07 15:26:31] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-11-07 15:26:44] Colors for breaks: #00008B,#24249B,#4848AB,#6D6DBC,#9191CC,#B6B6DD,#DADAEE,#FFFFFF,#EEDADA,#DDB6B6,#CC9191,#BC6D6D,#AB4848,#9B2424,#8B0000
## INFO [2018-11-07 15:26:44] Quantiles of plotted data range: 0.7,1,1,1.0017271440522,1.3## INFO [2018-11-07 15:26:45] plot_cnv_references:Writing observation data to ./infercnv.finalized_view.observations.txt
## INFO [2018-11-07 15:26:46] plot_cnv_references:Start
## INFO [2018-11-07 15:26:46] Reference data size: Cells= 42 Genes= 8828
## INFO [2018-11-07 15:26:46] plot_cnv_references:Number reference groups= 2
## INFO [2018-11-07 15:26:46] plot_cnv_references:Plotting heatmap.## INFO [2018-11-07 15:26:50] Colors for breaks: #00008B,#24249B,#4848AB,#6D6DBC,#9191CC,#B6B6DD,#DADAEE,#FFFFFF,#EEDADA,#DDB6B6,#CC9191,#BC6D6D,#AB4848,#9B2424,#8B0000
## INFO [2018-11-07 15:26:50] Quantiles of plotted data range: 0.7,1.0017271440522,1.0017271440522,1.0017271440522,1.3## INFO [2018-11-07 15:26:50] plot_cnv_references:Writing reference data to ./infercnv.finalized_view.references.txt## quartz_off_screen
## 2knitr::include_graphics("infercnv.finalized_view.png")
And that’s it. You can experiment with each step to fine-tune your data exploration. See the documentation for uploading the resulting data matrix into the Next Generation Clustered Heatmap Viewer for more interactive exploration of the infercnv-processed data: https://github.com/broadinstitute/inferCNV/wiki/Next-Generation-Clustered-Heat-Map