title: “infercnv_run” author: “Brian Haas” date: “9/5/2018” output: html_document: default pdf_document: default —
Create the InferCNV Object
infercnv_obj = CreateInfercnvObject(
raw_counts_matrix="sim.data",
annotations_file="sim.sample.annots.txt",
delim="\t",
gene_order_file="gencode_v19_gene_pos.txt",
ref_group_names=c("normal") )Filtering genes
Removing those genes that are very lowly expressed or present in very few cells
# filter out low expressed genes
cutoff=2
infercnv_obj <- require_above_min_mean_expr_cutoff(infercnv_obj, cutoff)## INFO [2018-09-11 12:01:20] ::above_min_mean_expr_cutoff:Start
## INFO [2018-09-11 12:01:20] ::process_data:Averages (counts).
## INFO [2018-09-11 12:01:20] Removing 16138 genes from matrix as below mean expr threshold: 2
## INFO [2018-09-11 12:01:20] validating infercnv_obj# 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-09-11 12:01:20] 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
infercnv_obj <- infercnv:::normalize_counts_by_seq_depth(infercnv_obj)perform Anscombe normalization
Suggested by Matan for removing noisy variation at low counts
infercnv_obj <- infercnv:::anscombe_transform(infercnv_obj)
save('infercnv_obj', file='infercnv_obj.anscombe')
plot_mean_chr_expr_lineplot(infercnv_obj)
log transform the normalized counts:
infercnv_obj <- log2xplus1(infercnv_obj)
save('infercnv_obj', file='infercnv_obj.log_transformed')Apply maximum bounds to the expression data to reduce outlier effects
threshold = mean(abs(get_average_bounds(infercnv_obj)))
infercnv_obj <- apply_max_threshold_bounds(infercnv_obj, threshold=threshold)## INFO [2018-09-11 12:01:22] ::process_data:setting max centered expr, threshold set to: +/-: 4.13549794979821perform smoothing across chromosomes
infercnv_obj = smooth_by_chromosome(infercnv_obj, window_length=101, smooth_ends=TRUE)## INFO [2018-09-11 12:01:22] ::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-09-11 12:01:25] ::center_smooth across chromosomes per cellsave('infercnv_obj', file='infercnv_obj.cells_recentered')
#plot_mean_chr_expr_lineplot(infercnv_obj)
plot_cnv(infercnv_obj, output_filename='infercnv.chr_smoothed', x.range="auto", title = "chr smoothed")## INFO [2018-09-11 12:01:26] ::plot_cnv:Start
## INFO [2018-09-11 12:01:26] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=-1261.60403741342 Min=-1.09939888679842 Max=2.26375113985155.
## INFO [2018-09-11 12:01:26] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-09-11 12:01:26] plot_cnv_observation:Start
## INFO [2018-09-11 12:01:26] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:01:26] clustering observations via method: average
## INFO [2018-09-11 12:01:26] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:01:42] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:01:43] plot_cnv_references:Start
## INFO [2018-09-11 12:01:43] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:01:43] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:01:43] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:01:47] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::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)## INFO [2018-09-11 12:01:49] ::subtract_ref_expr_from_obs:Startsave('infercnv_obj', file='infercnv_obj.ref_subtracted')
#plot_mean_chr_expr_lineplot(infercnv_obj, sep_obs_types = TRUE)plot_cnv(infercnv_obj, output_filename='infercnv.ref_subtracted', x.range="auto", title="ref subtracted")## INFO [2018-09-11 12:02:11] ::plot_cnv:Start
## INFO [2018-09-11 12:02:11] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=-4373.91939505554 Min=-0.991876041263262 Max=2.37127398538671.
## INFO [2018-09-11 12:02:11] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-09-11 12:02:11] plot_cnv_observation:Start
## INFO [2018-09-11 12:02:11] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:02:11] clustering observations via method: average
## INFO [2018-09-11 12:02:11] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:02:27] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:02:28] plot_cnv_references:Start
## INFO [2018-09-11 12:02:28] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:02:28] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:02:28] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:02:32] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::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)
save('infercnv_obj', file='infercnv_obj.inverted_log')
#plot_mean_chr_expr_lineplot(infercnv_obj, sep_obs_types = TRUE)
plot_cnv(infercnv_obj, output_filename='infercnv.inverted', color_safe_pal = FALSE, x.range="auto", x.center=1, title = "inverted log FC to FC")## INFO [2018-09-11 12:02:35] ::plot_cnv:Start
## INFO [2018-09-11 12:02:35] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=1316088.95692313 Min=0.502823491766739 Max=5.17397823583264.
## INFO [2018-09-11 12:02:35] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-09-11 12:02:35] plot_cnv_observation:Start
## INFO [2018-09-11 12:02:35] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:02:35] clustering observations via method: average
## INFO [2018-09-11 12:02:35] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:02:50] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:02:51] plot_cnv_references:Start
## INFO [2018-09-11 12:02:51] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:02:51] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:02:51] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:02:55] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::include_graphics("infercnv.inverted.png")
Removing noise
infercnv_obj <- clear_noise_via_ref_mean_sd(infercnv_obj, sd_amplifier = 1.5)## INFO [2018-09-11 12:02:57] :: **** clear_noise_via_ref_quantiles **** : removing noise between bounds: 0.901609084265536 - 1.10270432073886save('infercnv_obj', file='infercnv_obj.denoised')
#plot_mean_chr_expr_lineplot(infercnv_obj, sep_obs_types = TRUE)
plot_cnv(infercnv_obj, output_filename='infercnv.denoised', x.range="auto", x.center=1, title="denoised")## INFO [2018-09-11 12:02:57] ::plot_cnv:Start
## INFO [2018-09-11 12:02:57] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=1317838.1551309 Min=0.502823491766739 Max=5.17397823583264.
## INFO [2018-09-11 12:02:57] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-09-11 12:02:57] plot_cnv_observation:Start
## INFO [2018-09-11 12:02:57] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:02:57] clustering observations via method: average
## INFO [2018-09-11 12:02:57] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:03:10] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:03:11] plot_cnv_references:Start
## INFO [2018-09-11 12:03:11] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:03:11] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:03:11] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:03:15] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::include_graphics("infercnv.denoised.png")
Remove outlier data points
This generally improves on the visualization
infercnv_obj = remove_outliers_norm(infercnv_obj)## INFO [2018-09-11 12:03:16] ::remove_outlier_norm:Start out_method: average_bound lower_bound: NA upper_bound: NA
## INFO [2018-09-11 12:03:16] ::remove_outlier_norm using method: average_bound for defining outliers.save('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")## INFO [2018-09-11 12:03:17] ::plot_cnv:Start
## INFO [2018-09-11 12:03:17] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=1317750.43088188 Min=0.635456670448359 Max=2.05752707171781.
## INFO [2018-09-11 12:03:17] ::plot_cnv:Depending on the size of the matrix this may take a moment.
## INFO [2018-09-11 12:03:17] plot_cnv_observation:Start
## INFO [2018-09-11 12:03:17] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:03:17] clustering observations via method: average
## INFO [2018-09-11 12:03:17] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:03:31] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:03:32] plot_cnv_references:Start
## INFO [2018-09-11 12:03:32] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:03:32] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:03:32] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:03:35] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::include_graphics("infercnv.outliers_removed.png")
Tweaking the range to bring out the colors
plot_cnv(infercnv_obj, output_filename='infercnv.outliers_removed.adj', color_safe_pal = FALSE, x.range=c(0.7, 1.3), x.center=1, title = "outliers removed (color range adjustment)")## INFO [2018-09-11 12:03:37] ::plot_cnv:Start
## INFO [2018-09-11 12:03:37] ::plot_cnv:Current data dimensions (r,c)=5847,225 Total=1317750.43088188 Min=0.635456670448359 Max=2.05752707171781.
## INFO [2018-09-11 12:03:37] ::plot_cnv:Depending on the size of the matrix this may take a moment.## Warning in if (length(x.range) == 1 & x.range == "auto") {: the condition
## has length > 1 and only the first element will be used## INFO [2018-09-11 12:03:37] plot_cnv_observation:Start
## INFO [2018-09-11 12:03:37] Observation data size: Cells= 180 Genes= 5847
## INFO [2018-09-11 12:03:37] clustering observations via method: average
## INFO [2018-09-11 12:03:37] plot_cnv_observation:Writing observation groupings/color.## INFO [2018-09-11 12:03:50] plot_cnv_references:Writing observation data to ./observations.txt
## INFO [2018-09-11 12:03:51] plot_cnv_references:Start
## INFO [2018-09-11 12:03:51] Reference data size: Cells= 45 Genes= 5847
## INFO [2018-09-11 12:03:51] plot_cnv_references:Number reference groups= 1
## INFO [2018-09-11 12:03:51] plot_cnv_references:Plotting heatmap.## INFO [2018-09-11 12:03:55] plot_cnv_references:Writing reference data to ./references.txt## quartz_off_screen
## 2knitr::include_graphics("infercnv.outliers_removed.adj.png")