Richness

Rarefactions are available on the webdav: https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/

e168d33e79e0ebef4dc48d4341124315  tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r1000000.1.txt.bz2
b73412201512d47ec458ba1bd15a3047  tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2
7cf82f89c4ff969c495fc01e9556bb83  tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.2.txt.bz2
19d572c3105d030ee7de57facac9a286  tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.richnesses.txt.bz2

A Git repository is also available on request, or on OSC's intranet at http://anx191/gitweb/?p=F5/richness.git or osc-lm:/osc-fs_home/plessy/git/F5/richness.git.

Data preparation

Files

Checksum (MD5) of the files, downloaded from the webdav.

From https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/

509092f4ecdfdf47a2bf337a929f43a6  tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.txt.gz

This file contains 889 libraries. This is a subset of FREEZE_PHASE1 (988 libraries) where some libraries were discarded because of low quality or low count of tags.

CAGE peaks (DPI clusters).

Subsampling one million tags

In 2014, I have wrote a tutorial with similar commands on GitHub.

# Load human CAGE peaks (robust set, FREEZE_PHASE1).
library(tools)
md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.txt.gz') == '509092f4ecdfdf47a2bf337a929f43a6'
osc <- read.table("tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.txt.gz", row.names=1, head=TRUE)

# Substract the sum of the counts in the robust CAGE peaks, to the total tag count recorded in the '01STAT:MAPPED' row.
# The sums of the columns in the resulting matrix is the number of mapped tags.
osc['01STAT:MAPPED',] <- osc['01STAT:MAPPED',] - colSums(osc[-grep('01STAT:MAPPED', rownames(osc)),])

# Remove 57 libraries with less than 1,000,000 tags.
min.tags.per.lib <- 1000000
osc <- osc[,colSums(osc) >= min.tags.per.lib]

# Replace each library by a random sub-sample of 1,000,000 tags.
library(vegan)
osc <- t(rrarefy(t(osc),min.tags.per.lib))

# Backup the expression file.
write.table(file=paste('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r',min.tags.per.lib, '.1.txt', sep=''), osc)

Subsampling using the first quartile as treshold (currently being tested)

Same above, but using min.tags.per.lib <- quantile(colSums(osc), probs=0.25). This is a treshold of 2,514,503 tags that discards 222 libraries instead of 57.

Data analysis

Starting point using the saved rarefied data frames

Since the subsampling involves random numbers, each time the above commands are run the resulting data frame will be slightly different. Load the file from the webdav to work on the same numbers.

(Perhaps this may also be solved by setting the random seed.)

library(tools) # for md5sum
library(vegan) # for rarefy
md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2') == 'b73412201512d47ec458ba1bd15a3047'
osc <- read.table('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2')

Richness

# Rename libraries with their identifier only.
colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))

# Discard the '01STAT:MAPPED' row after the subsampling.
osc            <- osc[-grep('01STAT:MAPPED', rownames(osc)),]

# Calculate richness, maximum and sum expression for all the peaks.
richness   <- rarefy(osc, sample=10)
max_expr   <- apply(osc, 1, max)
sum_expr   <- rowSums(osc)

# Use a 10 tags treshold
treshold <- 10

hist(richness[sum_expr >= 10], br='fd')

Sent to Al for supplementary material

png('Richness_histogram.png')
hist(richness, br='fd', main="Expression breadth of robust CAGE peaks", xlab="Richness", ylab="number of peaks")
abline(v=median(richness), col='grey')
dev.off()

To quantify the breadth of expression of the robust CAGE peaks, we calculated a richness index (Hurlbert, 1971, The nonconcept of species diversity: a critique and alternative parameters. Ecology 52, 577-586.). This index represents the number of different libraries where the CAGE peaks would expected to be detected, if all the peaks would contain an arbitrary number of tags, here chosen to be 10. To avoid to under-estimate the expression breadth of peaks found in libraries that yielded less tags than average, and to obtain normalized expression values that are still round tag counts, we down-sampled without replacement each library to a total of one million mapped tags, using the rrarefy function of the R vegan package (Oksanen, J., F. et al., 2012. vegan: Community Ecology Package. R package version 2.0-3), discarding the 57 libraries where the total count was lower than one million. We also discarded 10,168 peaks that had less than 10 tags in total accross all the libraries after downsampling, as it is not possible to estimate their richness on scales smaller than 10. We then calculated the richness for the remaining 174,659 robust CAGE peak with a sample size of 10.

The distribution of richness indexes shows a peak for high values representing ubiquitously or very broadly expressed clusters. The median is 9.2 (vertical grey line), which is outside of the peak, showing that roughly half of the clusters are not ubiquitously expressed.

Potential shortcomings

The calculation of the richness score ignores the zero values. For instance, on a scale of 10, the following clusters would have the maximum score (10).

cluster1: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
cluster2: 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0
cluster3: 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0

The solution is to down-sample on a scale that is at least as high as the number of libraries: in that case, there is more than one tag to distribute per library, so zero values are echoed by increased tag counts in other libraries.

However, in order to properly evaluate the expression of the clusters, we need to consider only those where at least one library has as many tags as the sampling scale used in the richness calculations.

Combining both previous constraints, this places a cut-off at 832 TPM, which is way too high.

I tried to circumvent this by using the 49 library co-expression clusters ("Kenny's clusters"). The smallest cluster contains two libraries. Therefore, I can make pools of 2 million tags. This makes a cut-off of 49 / 2 = 24.5 TPM. The result is that clusters that have that high expression levels in at least one library are uniformly highly expressed regardless of their richness.

Other calculations not used in the final analysis

Files

Checksum (MD5) of the files, downloaded from the webdav.

From https://fantom5-collaboration.gsc.riken.jp/webdav/home/mette/

3d2357037e25987f37e1c7f28875ffaf  robust_dominant_+-100.bed
b7284de4bf79dd50ca9aea4339ea15d3  robust_dominant_+-100.fa_0.7_TATA-Box.res

From https://fantom5-collaboration.gsc.riken.jp/webdav/home/mdehoon/MotifActivity/hg19

be5f523a574b928402b0d6123085d4a4  hg19.network.txt.bz2

From https://fantom5-collaboration.gsc.riken.jp/webdav/home/baillie/Clustering_by_cell_type/ and https://fantom5-collaboration.gsc.riken.jp/webdav/home/baillie/V3_clustering_p75_m22/

691ef45f34f96ed4d38c003df3467caf  SYNC015_3_tpm_transposed_p80_m22.layout
4babb31b79799043f7068d37cb701db2  SYNC015_combined_clustertables.txt

MARA network.

# Load human MARA network.
md5sum('hg19.network.txt.bz2') == 'be5f523a574b928402b0d6123085d4a4'
net <- read.table("hg19.network.txt.bz2")
# Remove "FOX.C1.C2.", "GATA1..3", "KLF12" and "MSX1.2", for which all z-scores are null.
net <- net[,colSums(net) > 0]
# Fit the net and osc data frames.
net <- net[rownames(osc),]

Mette's motifs

Data preparation in R

# Load peak list.  This is needed as the motif files refer to the peaks by their position in the list.
md5sum( 'robust_dominant_+-100.bed') == '3d2357037e25987f37e1c7f28875ffaf'
mette.p           <- read.table('robust_dominant_+-100.bed', head=F, stringsAsFactors=F)[,4]

# Load position and scores for TATA motif.
md5sum( 'robust_dominant_+-100.fa_0.7_TATA-Box.res') == 'b7284de4bf79dd50ca9aea4339ea15d3'
mette.tata.df     <- read.table('robust_dominant_+-100.fa_0.7_TATA-Box.res', head=T, stringsAsFactors=F)

# Score cutoff, completely arbitrary for the moment.
mette.tata.cutoff <- 10

# Names of the peaks with a motif above the cutoff.
mette.tata        <- mette.p[with(mette.tata.df, unique(tcIndex[Score > mette.tata.cutoff]))]

# Removing names of the peaks discarded earlier. 
mette.tata        <- mette.tata[mette.tata %in% rownames(osc)]

Cutoffs

# Chose a z-score cutoff and a corresponding output directory.
# All non-null scores (already filtered in hg19.network.txt.bz2).
cutoffs <- apply(net,2,function (X) min(X[X >0]));  cutoff_dir <- "cutoff_0/"
#cutoffs <- apply(net,2,function (X) as.numeric(summary( X[X >0])[2])); cutoff_dir <- "cutoff_Q1/" # 1st Quartile.
#cutoffs <- apply(net,2,function (X) as.numeric(summary( X[X >0])[4])); cutoff_dir <- "cutoff_Q2/" # 2nd Quartile.
#cutoffs <- apply(net,2,function (X) as.numeric(summary( X[X >0])[5])); cutoff_dir <- "cutoff_Q3/" # 3rd Quartile.

Plotting

The following commands, inspired from Al's presentation sent on March 14th (fantom5:01462) are needed to setup the parameters for the smoothScatter plots.

library(RColorBrewer)
brewer.pal(11, "RdYlBu") 
buylrd <- c("#313695", "#4575B4", "#74ADD1", "#ABD9E9", "#E0F3F8", "#FFFFBF",  "#FEE090", "#FDAE61", "#F46D43", "#D73027", "#A50026")

smoothScatterRed <- function (clusters, richness, max_expr, samplesize, TF) {
  smoothScatter(
    richness[clusters],
    log10(max_expr[clusters]),
    nbin=600, colramp = colorRampPalette(c(buylrd)),
    xlab='Richness', ylab='log10(Max TPM)',
    main=paste(TF, ' (', length(clusters), ' clusters)', sep=''),
    xlim=c(0,samplesize), ylim=c(0,5.5))
}

plotTF <- function (TF) {
  clusters <- rownames(net)[net[,TF] > cutoffs[TF]]
  smoothScatterRed(clusters=clusters, richness=richness, max_expr=max_expr, samplesize=100, TF=TF)
}

pngTF <- function (TF) {
  png(paste(cutoff_dir,TF,".png", sep=''))
  plotTF(TF)
  dev.off() }

plotsTF <- function (TF) {
  clusters <- rownames(net)[net[,TF] > cutoffs[TF]]
  smoothScatterRed(clusters=clusters, richness=richness.s, max_expr=max_expr.s, samplesize=49, TF=TF)
}

Plot all with for (TF in names(cutoffs)) plotTF(TF)

A bit of stats...

t.test(richness[net$TBP > 0], richness[net$SP1 > 0])

        Welch Two Sample t-test

data:  richness[net$TBP > 0] and richness[net$SP1 > 0]
t = -21.4623, df = 991.734, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.46669 -18.70248
sample estimates:
mean of x mean of y
 55.77583  76.36041

Comparison of two replicates

Richness calculated for each replicate was saved as richness.1 and richness.2 respectively.

x <- rbind(
  data.frame(cluster=names(richness.1), id="richness.1", value=as.numeric(richness.1)),
  data.frame(cluster=names(richness.2), id="richness.2", value=as.numeric(richness.2)),
  data.frame(cluster=names(max_expr.1), id="max_expr.1", value=as.numeric(max_expr.1)),
  data.frame(cluster=names(max_expr.2), id="max_expr.2", value=as.numeric(max_expr.2)))
library(reshape)
richnesses <- cast(x, cluster ~ id)
rownames(richnesses) <- richnesses$cluster
richnesses <- richnesses[,-1]
write.table(file="tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.richnesses.txt.bz2", richnesses)

Graphical comparison of the two replicates.

with(richnesses,
  smoothScatter(
   richness.1, richness.2,
   main="Richness compared for two independant sub-samplings.",
   xlab="Sub-sampling 1", ylab="Sub-sampling 1",
  nbin=1000, colramp = colorRampPalette(c(buylrd)))
)

Variant where libraries in the same co-expression groups are collapsed (analysis ongoing)

Sample groups

DBscan clusters from Timo.

ab34b28705e5c945decf5c15ab4ca753  F5_sample_clustering.xls

Sheet 2 was exported as tabulation-separated values and the sample names were simplified to library IDs with the following command.

sed -i -e 's/\t.*CNhs/\tCNhs/' -e 's/\..*//' F5_sample_clustering.tsv

Merge similar libraries

library(tools) # for md5sum
library(vegan) # for rarefy

# Load the definition of the clusters.

md5sum('F5_sample_clustering.tsv') == '297c2e5edcafdebcf991937d6e0e8fe7'
sampleclusters <- read.table('F5_sample_clustering.tsv', head=F, col.names=c('cluster','sample'), colClasses=c('character', 'character'))
rownames(sampleclusters) <- sampleclusters$sample

# Load the tag counts of the rarefied libraries.

md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2') == 'b73412201512d47ec458ba1bd15a3047'
osc <- read.table('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2')
colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))

# Keep only the 832 common libraries
sampleclusters <- sampleclusters[intersect(colnames(osc),rownames(sampleclusters)),]

# 595/655 clusters have only one library (734/796 before matching libraries).
table(table(sampleclusters$cluster) == 1)
# FALSE  TRUE 
#    60   595 

# Pool the libraries that are in the same cluster.
osc <- t(rowsum(t(osc[,rownames(sampleclusters)]),sampleclusters$cluster))

# Rarefy as above.
min.tags.per.lib <- min(colSums(osc))
osc <- t(rrarefy(t(osc), min.tags.per.lib))

# Only 640 DPI clusters were lost in the process.
table(rowSums(osc) == 0)
#  FALSE   TRUE 
# 184188    640 

# Backup the expression file.
write.table(file=paste('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r',min.tags.per.lib, '.1.txt', sep=''), osc)

md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r1000000.1.txt') == 'dd7c8f19b4705f05dc2f1283dfb8057f'

md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r2514503.1.txt') == '1864303f6526854f9377a25112616a3c'

Analysis

md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r1000000.1.txt') == 'dd7c8f19b4705f05dc2f1283dfb8057f'
osc <- read.table('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r1000000.1.txt')

# Discard the '01STAT:MAPPED' row after the subsampling.
osc <- osc[-grep('01STAT:MAPPED', rownames(osc)),]

# Discard clusters of total expression lower than 10 tags (the scale of the richness sampling).
osc <- osc[apply(osc, 1, max) >= 10,] # 94,875 peaks discarded, 89,953 remain

# Calculate richness, maximum and median expression for all the remaining peaks.
richness <- rarefy(osc, sample=10)

png('Richness_DBscan-sample-coexpr-clusters_histogram.png')
hist(richness, br='fd', main="Expression breadth of robust CAGE peaks (DBscan sample clusters)", xlab="Richness", ylab="number of peaks")
abline(v=median(richness), col='grey')
dev.off()

In order to avoid low-richness artefacts, the subsampling must be done on a scale larger or equal to the number of libraries (not the case here...), and on promoters where at least one library has more tags than that scale. Here it makes a 10 TPM cutoff.

Variant with clusters (not used in final analysis)

See [fantom5:01529].

clustered.promoters <- read.table('SYNC015_combined_clustertables.txt', header=TRUE, row.names=1)
i <- intersect(intersect(names(richness),rownames(clustered.promoters)), rownames(net))
coexpr.clusters <- data.frame(mean=tapply(richness[i], clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt'], mean))
coexpr.clusters$sd <- tapply(richness[i], clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt'], sd)
coexpr.clusters$n <- tapply(richness[i], clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt'], length)
coexpr.clusters <-coexpr.clusters[order(coexpr.clusters$n, decreasing=TRUE),]
max_expr.c <- tapply(max_expr[i], clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt'], max)
net.c <- aggregate(net[i,],list(clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt']),mean)
rownames(net.c) <- net.c[,1]
net.c <-  net.c[,-1]

sampleclusters.c <- rowsum(sampleclusters[i,], clustered.promoters[i,'SYNC015_3_tpmALL_no_pooled_p75_m22.clustertable.txt'])
srichness.c <- rarefy(sampleclusters.c, 49)
smax_expr.c <- apply(sampleclusters.c, 1, max)

Cutoffs

cutoffs.c <- apply(net.c,2,function (X) min(X[X >0]));  cutoff_dir <- "cutoff_0/"
#cutoffs.c <- apply(net.c,2,function (X) as.numeric(summary( X[X >0])[2])); cutoff_dir <- "cutoff_Q1/" # 1st Quartile.
#cutoffs.c<- apply(net.c,2,function (X) as.numeric(summary( X[X >0])[4])); cutoff_dir <- "cutoff_Q2/" # 2nd Quartile.
#cutoffsc. <- apply(net.c,2,function (X) as.numeric(summary( X[X >0])[5])); cutoff_dir <- "cutoff_Q3/" # 3rd Quartile.

Plotting

The following commands, inspired from Al's presentation sent on March 14th (fantom5:01462) are needed to setup the parameters for the smoothScatter plots.

plotTFc <- function (TF) {
  clusters <- rownames(net.c)[net.c[,TF] > cutoffs.c[TF]]
  smoothScatter(
    coexpr.clusters[clusters,'mean'],
    log10(max_expr.c[clusters]),
    nbin=600, colramp = colorRampPalette(c(buylrd)),
    xlab='Richness', ylab='log10(Max TPM)',
    main=paste(TF, ' (', length(clusters), ' clusters)', sep=''),
    xlim=c(0,100), ylim=c(0,5.5)) }

plotsTFc <- function (TF) {
  clusters <- rownames(net.c)[net.c[,TF] > cutoffs.c[TF]]
  smoothScatter(
    srichness.c[clusters],
    log10(smax_expr.c[clusters]),
    nbin=600, colramp = colorRampPalette(c(buylrd)),
    xlab='Richness', ylab='log10(Max TPM)',
    main=paste(TF, ' (', length(clusters), ' clusters)', sep=''),
    xlim=c(0,49), ylim=c(0,5.5)) }

User:Plessy/Richness

Contents

Richness

Data preparation

Files

CAGE peaks (DPI clusters).

Subsampling one million tags

Subsampling using the first quartile as treshold (currently being tested)

Data analysis

Starting point using the saved rarefied data frames

Richness

Sent to Al for supplementary material

Potential shortcomings

Other calculations not used in the final analysis

Files

MARA network.

Mette's motifs

Cutoffs

Plotting

A bit of stats...

Comparison of two replicates

Variant where libraries in the same co-expression groups are collapsed (analysis ongoing)

Sample groups

Merge similar libraries

Analysis

Variant with clusters (not used in final analysis)

Cutoffs

Plotting

Navigation menu

Page actions

Page actions

Personal tools

Menu

Search

Special topics

Resources

ZENBU genome browser

UCSC Genome Browser RIKEN mirror

Navigation

Tools