-
Notifications
You must be signed in to change notification settings - Fork 4
Signal from telomeres
There is a dedicated function to collect signal from telomeres in hwglabr2
. It collects signal from all telomeres inward up to a specified distance.
Here is an example of how you can use it with wild type data mapped to "sacCer3". Start by importing the signal track data, here input-normalized Red1 ChIP-seq signal mapped to the sacCer3 reference genome assembly.
# Load signal track data
signal_file <- paste0('/Volumes/LabShare/HTGenomics/HiSeqOutputs/AveReps_SacCer3_MACS2_FE/',
'Red1-wildtype-71-34-199-29-Reps-SacCer3-B3W3-MACS2/',
'Red1-wildtype-71-34-199-29-Reps-SacCer3-2mis_B3W3_MACS2_FE.bdg.gz')
signal <- hwglabr2::import_bedGraph(signal_file)
You can now run the signal_from_telomeres2
function to collect the signal. The length to collect defaults to 100 kb.
# Collect signal
signal_at_tels <- hwglabr2::signal_from_telomeres2(signal,
length_to_collect = 100000,
genome = 'sacCer3')
This outputs a data frame with one row per chromosome arm and a variable number of columns:
-
chr
chromosome number -
arm
chromosome arm ("L" or "R") -
size_cat
chromosome size category ("small" or "large") -
t1:tn
signal columns for each position; number depends on the specified collection length (n = length_to_collect
)
You can also use parameter averaging_window
to get average signal at windows of the specified size, which will reduce the resolution and run faster. Here we are using the default, averaging_window=1
, to get signal at every position.
You now have signal at telomeres for all 32 telomeres in the genome. You can average different groups together and plot according to your needs. For example, let's say we want to plot average signal at all small and large chromosome telomeres normalized to the mean genome-wide signal. We would start by getting the signal columns and computing the average signal.
# Get signal columns for small and large chromosomes
small_chrs <- signal_at_tels[signal_at_tels$size_cat=='small',
4:ncol(signal_at_tels)]
large_chrs <- signal_at_tels[signal_at_tels$size_cat=='large',
4:ncol(signal_at_tels)]
# Calculate average signal
small_chrs <- colMeans(small_chrs, na.rm = T)
large_chrs <- colMeans(large_chrs, na.rm = T)
We can then normalize the signal to the genome-wide mean. The genome-wide mean can be obtained using function average_chr_signal
in the package. This function outputs both the individual chromosome mean signal values and the genome-wide mean (second element in the output list). Once the mean is computed we can simply divide the signal vectors computed above by this value.
# Get genome-wide mean
genome_wide_mean <- hwglabr2::average_chr_signal(signal)[[2]]
# Normalize signal
small_chrs <- small_chrs / genome_wide_mean
large_chrs <- large_chrs / genome_wide_mean
Finally, we can build a quick plot to visualize the signal.
# Plot
plot(x=seq(1, 100000), y=small_chrs, col='orange',
xlab='Distance to chr end', ylab='Average signal (genome mean = 1)',
type='l')
lines(x=seq(1, 100000), y=large_chrs, col='red')
abline(h=1, lty=3)