2013-08-24T08:38:50Z

Plessy: Rename output page.

2013-03-23T10:20:39Z

Plessy: Minor correction made by Al in December.

= Richness =

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

<pre>
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
</pre>

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

== 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/
<pre>
509092f4ecdfdf47a2bf337a929f43a6 tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.txt.gz
</pre>

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). ===

<pre>
# 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="tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt", osc)
</pre>

== 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.)''

<pre>
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')
</pre>

=== Richness ===

<pre>
# 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)),]

# Discard clusters of total expression lower than 10 tags (the scale of the richness sampling).
osc <- osc[rowSums(osc) >= 10,] # 10,168 peaks discarded

# Calculate richness, maximum and median expression for all the remaining peaks.
richness <- rarefy(osc, sample=10)
max_expr <- apply(osc, 1, max)
</pre>

=== Sent to Al for supplementary material ===

[[File:Richness_histogram.png‎]]

<pre>
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()
</pre>

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.

== 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/
<pre>
3d2357037e25987f37e1c7f28875ffaf robust_dominant_+-100.bed
b7284de4bf79dd50ca9aea4339ea15d3 robust_dominant_+-100.fa_0.7_TATA-Box.res
</pre>

From https://fantom5-collaboration.gsc.riken.jp/webdav/home/mdehoon/MotifActivity/hg19
<pre>
be5f523a574b928402b0d6123085d4a4 hg19.network.txt.bz2
</pre>

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/

<pre>
691ef45f34f96ed4d38c003df3467caf SYNC015_3_tpm_transposed_p80_m22.layout
4babb31b79799043f7068d37cb701db2 SYNC015_combined_clustertables.txt
</pre>

=== MARA network. ===

<pre>
# 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),]
</pre>

=== Mette's motifs ===

Data preparation in R

<pre>
# 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)]
</pre>

=== Cutoffs ===
<pre>
# 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.
</pre>

== Plotting ==

The following commands, inspired from [https://fantom5-collaboration.gsc.riken.jp/webdav/home/Al.Forrest/tissue_specificity_feb2012_human180Krobust.pdf Al's presentation] sent on March 14th (fantom5:01462) are needed to setup the parameters for the <code>smoothScatter</code> plots.

<pre>
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)
}
</pre>

Plot all with <code>for (TF in names(cutoffs)) plotTF(TF)</code>

== A bit of stats... ==
<pre>
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
</pre>

== Comparison of two replicates ==

Richness calculated for each replicate was saved as <code>richness.1</code> and <code>richness.2</code> respectively.

<pre>
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)
</pre>

Graphical comparison of the two replicates.

<pre>
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)))
)
</pre>

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

=== Sample groups ===

''DBscan'' clusters from Timo.

<pre>
ab34b28705e5c945decf5c15ab4ca753 F5_sample_clustering.xls
</pre>

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

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

=== Merge similar libraries ===

<pre>
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.
osc <- t(rrarefy((t(osc)), 1000000))

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

# Backup the expression file.
write.table(file="tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r1000000.1.txt", osc)

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

=== Analysis ===

<pre>
md5sum('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r2000000.1.txt.bz2') == '8d9e5cd5d48a629d98bfcea227c0b8d1'

sampleclusters <- read.table('tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.pooled.r2000000.1.txt.bz2')
sampleclusters <- sampleclusters[-grep('01STAT:MAPPED', rownames(sampleclusters)),]
sampleclusters <- sampleclusters[rowSums(sampleclusters) > 0,] # 17,203 peaks discarded.

sampleclusters <- sampleclusters[apply(sampleclusters,1,max) >= 49,] # 147,631 more peaks discarded ( 19,993 remaining).
#sampleclusters <- sampleclusters[rowSums(sampleclusters) >= 49, ] # 113,343 more peaks discarded ( 54,281 remaining).

richness.s <- rarefy(sampleclusters, sample=49)
max_expr.s <- apply(sampleclusters, 1, max) / 2 # We have 2,000,000 tags per library.

png('Richness_sample-coexpr-clusters_histogram.png')
hist(richness.s, br='fd')
dev.off()

png('Max-expr_Richness_sample-coexpr-clusters.png')
boxplot(log10(max_expr.s) ~ factor(ceiling(richness.s)), xlab='richness', ylab='log10(TPM)', main='B: Max expression level, binned by richness')
dev.off()
</pre>

[[File:Richness_sample-coexpr-clusters_histogram.png]]
[[File:Max-expr_Richness_sample-coexpr-clusters.png‎]]

In order to avoid artefacts, the subsampling must be done on a scale larger or equal to the number of libraries, and on promoters where at least one library has more tags than that scale. Here it makes a rough 25 TPM cutoff. Unfortunately, the conclusion is that clusters that have that high expression levels in at least one library are uniformly highly expressed regardless of their richness.

Clusters expressed only in exactly one co-expression group (that is: undetectable in the other groups) have a richness of exactly one and tend to have low expression scores.

<pre>
> summary(richness.s)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 11.86 22.56 19.72 28.13 30.75
</pre>

== Variant with clusters (not used in final analysis) ==

See <code> [fantom5:01529]</code>.

<pre>
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]
</pre>

<pre>
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)
</pre>

=== Cutoffs ===
<pre>
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.
</pre>

=== Plotting ===

The following commands, inspired from [https://fantom5-collaboration.gsc.riken.jp/webdav/home/Al.Forrest/tissue_specificity_feb2012_human180Krobust.pdf Al's presentation] sent on March 14th (fantom5:01462) are needed to setup the parameters for the <code>smoothScatter</code> plots.

<pre>
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)) }
</pre>

File:Richness DBscan-sample-coexpr-clusters histogram.png

2013-03-23T10:07:26Z

Plessy: uploaded a new version of "File:Richness DBscan-sample-coexpr-clusters histogram.png": Added abline.

2013-01-19T06:01:02Z

Plessy: Sort the chromosomes in alphanumeric order (for BEDTools).

= Charles Plessy =

= Works on =

* [[nanoCAGE]]
* [[CAGEscan]]
* [[Olfactory receptors]]
* [[OP-SOLEXA-5RACE-v1.0|Deep-RACE]]
* [[Extend rat gene models with CAGEscan]]

== [[Painting]] and other strange CAGE signals ==

Jessica pointed out that in Fantom5, some hemoglobin genes single out by the absence of locus painting.

https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=d1yxe9xqKJDltVi86B6CFD;loc=mm9::chr11:32168803..32199353

Question: when Timo aligns promoter sequences in which errors were introduced, some reads map on introns and exons. Do we find more of them in intron and exon-painted genes ?

How is it on polysome fractions ?

Need CAGE libraries without cap trapping.

A page is created to stimulate collaboration, [[Painting]].

Is the quantity of tags aligning to the [[rDNA]] correlated to the expression of the [[Ribosome|ribosomal proteins]] ?

= Shared files =

== Whole-genome BED files ==

These BED files recapitulate the expression of all libraries for each genome, and are available on the [[Dataset_introduction#How_to_share_my_analysis_results_with_the_consortium_member.3F|webdav]] for the [[File_release|file releases]]
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1.1/per-genome FREEZE_PHASE1.1] and
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE2/per-genome FREEZE_PHASE2], and previous updates
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_11/per-genome 011],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/per-genome 012],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_13/per-genome 013],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_014/per-genome 014]
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_015/per-genome 015], and
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_016/per-genome 016]. They are sorted by position, [http://samtools.sourceforge.net/tabix.shtml compressed with bgzip and indexed with tabix]. In CTSS files, column 4 contains the library ID, and the genome name in TSS files. For example, here is the result of the command <code>tabix rn4.bed.gz chr1:1923214-1963638</code>.

<pre>
chr1 1923214 1923215 CNhs11312 1 +
chr1 1923717 1923718 CNhs12070 1 +
chr1 1931875 1931876 CNhs10614 1 -
chr1 1963635 1963636 CNhs12070 1 -
chr1 1963636 1963637 CNhs11298 1 -
chr1 1963636 1963637 CNhs12070 3 -
chr1 1963637 1963638 CNhs11298 1 -
chr1 1963637 1963638 CNhs12070 2 -
</pre>

It indicates for instance, that in the library CNhs12070, there was 1 tag on chromosome 1 at position 1963635, 3 at position 1963636 and 2 at position 1963637. Conversely, it indicates that at position 1963636 on chromosome 1, there are 3 tags in library CNhs12070 and 2 tag in CNhs11298.

These files can take a particular advantage of the <code>groupby</code> command from the [https://github.com/arq5x/filo FILO] package:

<pre>
tabix rn4.bed.gz chr1:1923214-1963638 | sort -k4 | groupBy -g 4 -c 5 -ops sum
CNhs10614 1
CNhs11298 2
CNhs11312 1
CNhs12070 7
</pre>

=== Making of ===

==== Interrogating the list of all the HeliScopeCAGE libraries ====

The CNhs_lib function will be needed below.

This needs a FANTOM5 mirror or a flat directory with all [[Dataset_introduction#Where_can_I_get_all_the_details_of_the_libraries.3F|SDRF]] files. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*sdrf.txt</code> is useful for that purpose. The following [http://en.wikipedia.org/wiki/Bash_(Unix_shell) bash] function interrogates them, looking by default in <code>~/public_html/F5/sdrf</code>, but this can be overridden by setting the <code>F5</code> variable.

<pre>
# Return all libraries.
function SDRF_lib {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f14 | grep CNhs
done
}

# Returns only HeliScopeCAGE libraries.
function CNhs_lib {
for CNhs in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $CNhs | cut -f14 | grep CNhs
done
}

# Returns the description of a library.
function CNhs-desc {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f3
done
}
</pre>

==== Flat directory with modified BED CTSS files. ====

Original BED CTSS files are expected to be downloaded in a local directory called <code>longnames</code>. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*bed.gz</code> is useful for that purpose.

The following commands produces modified BED CTSS files, where the <code>name</code> field contain the library name instead of the coordinates. Beware that this format does not allow for the co-existence of more than one alignment file per library. In what follows, the galGal3 alignments are discarded and the galGal4 alignments are kept.

<pre>
# Make sure there are no duplicata.
ls longnames | cut -f2 -d '.' | sort | uniq -d

mkdir -p bgrezip

for LIB in $(ls longnames | cut -f2 -d '.')
do
echo -ne "${LIB}.ctss.bed.gz\t"
zcat longnames/*${LIB}*ctss.bed.gz |
awk -v LIB=$LIB '{OFS="\t"} {$4=LIB ; print}' |
sort -k1,1 -k2,2n |
bgzip |
tee bgrezip/${LIB}.ctss.bed.gz |
md5sum |
cut -f1 -d' '
done |
tee bgrezip.txt

# Index with tabix.
(cd bgrezip;
for BED in *ctss.bed.gz
do
tabix -p bed $BED
done)
</pre>

''Note that wikimedia inserts a strange non-breakable character in the <code>awk</code> command. Paste it from the source.''

Collect the files in the <code>bgrezip</code> directory to replace the original ones.

Note: perhaps <code>-k6,6</code> would make the sorting completely deterministic, but this may not be needed if the input already guarantees that for bases where both strands are transcribed, it is always either the <code>+</code> or <code>-</code> strand that comes first.

==== Construction of whole-genome BED CTSS files ====

Now that the <code>name</code> field contains the library name, the data can pooled in one file per genome.

<pre>
# Tip: run this with a low IO priority:
# exec ionice -c 3 bash

# Mouse
CHRLIST='chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrM chrX chrY'
ORGANISM='Mus musculus'
OUTPUT=mm9.ctss.bed.gz

rm -rf chrtmp
mkdir chrtmp

for CHR in $CHRLIST
do
for BED in $( CNhs_lib "$ORGANISM" )
do
if [ -e ${BED}.ctss.bed.gz ]
then
echo -e "${CHR}\t${BED}" 1>&2
tabix ${BED}.ctss.bed.gz ${CHR}
# zgrep -P "${CHR}\t" ${BED}.ctss.bed.gz # Delete after testing with the above, which is almost 10 times faster.
fi
done |
sort --temporary-directory=$(pwd) -k2,2n > chrtmp/${CHR} # Make sure there is enough space for sort's temporary files !
done
cd chrtmp
cat $CHRLIST |
bgzip > $OUTPUT
</pre>

<pre>
CHRLIST='chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 chr20 chr21 chr22 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrM chrX chrY'
ORGANISM='Homo sapiens'
OUTPUT=hg19.ctss.bed.gz
</pre>

<pre>
for LIB in $(CNhs_lib 'Canis lupus familiaris'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > canFam2.ctss.bed.gz
for LIB in $(CNhs_lib gallus); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > galGal3.ctss.bed.gz
for LIB in $(CNhs_lib 'Rattus'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > rn4.ctss.bed.gz
</pre>

==== Construction of whole-genome BED TSS files ====

Pooling all libraries reduces the size compressed files.

<pre>
for GENOME in hg19 mm9 rn4 canFam2 galGal3
do
zcat $GENOME.ctss.bed.gz |
groupBy -g 1,2,3,6 -c 5 -o sum |
awk -v GENOME=$GENOME '{OFS="\t"} {print $1, $2, $3, GENOME, $5, $4}' |
bgzip > $GENOME.tss.bed.gz
done
</pre>

==== In Zenbu ====

The whole-genome BED TSS files for human and mouse are loaded in [[ZENBU_Howto|Zenbu]]. Because the library information is collapsed, it is possible to display large intervals quiclkly, or even whole chromosome (~30 s for human chromosome 19).

== Indexed expression files ==

Source:

* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/110804-dpi-clusters-expression/00readme.html (UPDATE_012)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/00readme.html (FREEZE_PHASE1)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/description120126/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz

Four columns were added for sorting and indexing: ''chromosome'', ''start'', ''end'' and ''strand''.

The files for FREEZE_PHASE1 are available at https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/expression and was indexed with [http://samtools.sourceforge.net/tabix.shtml tabix] (see also [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/expression UPDATE_012]).

=== Making of ===
<pre>
TABLE=tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrM chrX chrY'
zgrep '^##' -A3 $TABLE > $(basename $TABLE .osc.txt.gz).header
# Edit the header to add four columns, named ‘chrom start end strand’.
for chrom in $CHRLIST
do
zgrep "^$chrom:" $TABLE |
perl -ne '/^(chr[\d\w]+):(\d+)..(\d+),([+-])/; print join "\t", $1, $2, $3, $4, $_' |
sort --field-separator "$(printf "\t")" -k2,2n -k3,3n -k4,4 > $chrom
done
cat $(basename $TABLE .osc.txt.gz).header $CHRLIST | bgzip > $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz
tabix $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz -s 1 -b 2 -e 3 -S 3 -0
</pre>

For mouse use the following chromosome list instead.
<pre>
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'
</pre>

= Draft notes / Brainstorm =

Assemble chrM and compare recorded ethnicity with mitochondrial haploypes ?

Ultraconserved elements: http://users.soe.ucsc.edu/~jill/ultra.html

* Liftover: http://genome.ucsc.edu/cgi-bin/hgLiftOver
* Or map the fasta files ?
* Or recalculate within hg19 ?

Do we see chromosomal aberrations in the CAGEscan data?

Search for CDR3 regions in the longest reads.

Can we basecall SNPs from the Helicos data ?

Sense-antisense pair at VIM ? HG19::chr10:17270611-17272991

Orphan promoter: HG19::chr11:27908185-27913031

Another orphan promoter, expressed only in nervous tissue: HG19::chr11:28891298-29036678

Long chain: HG19::chr11:31253864-31864905

Cancer cells seem to express more snoRNAs in the ''RPS3'' locus than other cells; it may be interesting to test if they express more ribosomal components in general, perhaps because they are fast-growing.

* hg19::chr15:66816186..66817714 : alternative RPL4 promoter ?

* Strong painting in ''IQSEC1''. https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:12915663..13276777
** And in ''NUP210''.

== Reads and CTSS distributions in arbitrary windows ==

[[File:plessy-Antisense-RefSeqTSS-CTSS-w2000.png|thumb|See the thread fantom5:00740 for more details]]

<pre>
PER_GENOME=/home/charles/public_html/F5/tabix_12
HALF_WINDOW=5000
WINDOW=$(( $HALF_WINDOW * 2 ))

slopBed -i refGeneTss.bed -g /usr/share/bedtools/genomes/human.hg19.genome -l $HALF_WINDOW -r $HALF_WINDOW |
while read chr start end name score strand
do
tabix $PER_GENOME/hg19.bed.gz ${chr}:$(( $start + 1))-${end} |
grep -v "${strand}" |
awk -v start=$start -v end=$end -v strand=$strand '{OFS="\t"} {if (strand == "+") print $2 - start, $5 ; else print end - $3, $5 }' |
groupBy -c 2 -g 1 -o sum -inheader |
perl -ne 'print unless /^0$/'
done > window-$WINDOW.txt

sort -nk1,1 window-$WINDOW.txt | groupBy -g 1 -c 2 -o sum > window-$WINDOW.profile

cut -f1 window-$WINDOW.txt | sort -n | uniq -c | sort -nk2 > window-$WINDOW.mapcount
</pre>

== [[CAGEscan]] ==

* ''MLT1A'' repeat bridging with ''EYPC'': https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=n3gh7Y8TOJGkOqG_q9EZB;loc=hg19::chr12:91344715..91439984

* Why did we lose this promoter in ''IQSEC1''? https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:13007279..13012828

== [[User:Plessy/Richness|Richness]] ==

[[File:Richness_FREEZE1_TBP.png|thumb|TBP]][[File:Richness_FREEZE1_SP1.png|thumb|SP1]]

In the <code>R</code> package called <code>vegan</code> there is a function, called <code>rarefy</code>, to reduce the data as if it had contained a fixed number of tags. Used on the libraries, it evaluates how complex they are (genome-bioinfo:02933). Another function, <code>rrarefy</code> makes a random sub-sample of a library with an arbitrary number of counts, respecting the original distribution. See [[User:Plessy/Richness]] for more details.

I calculated the ''richness'' of the robust clusters in FREEZE1 after normalising the libraries to one million tags by random rarefaction. The richness evaluates how ubiquitous the clusters are. For instance, with a sampling size of 100, the clusters with a TBP motif have a richness of 55.8 ± 29.2 (''n''=948) and the clusters with a SP1 motif have a richness of 76.4 ± 18.4 (''n''=16,134), a significant difference (''p'' < 2.2e-16).

The profiles for each motif are available at the following URL: [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/ https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/], as well as the summary panels [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness//T1.png T1.png] to [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/T8.png T8.png]. The sub-sampled data used for this analysis is available in the file [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2 tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2]

=== To Do ===

* Cross data with Mette's motifs and Sebastian's houskeeping list <code>fantom5:01568</code>
* Filter CAGE peak robust set with predicted TSS. <code>fantom5:01581</code>

= Misc =

* [[User:Plessy/Debian|Software used in FANTOM5 and available in the Debian operating system.]]

* Sorting a ''level1'' OscTable called OSCTABLE: <code>cat <(grep \# -A1 OSCTABLE) <(grep -v \# OSCTABLE | sed '1d' | sort --field-separator "$(printf "\t")" -k2.4,2n -k 2.4,2.4 -k3,3n -k4,4n -k5,5) | sponge OSCTABLE</code>

* Indexing a sorted and bgzipped ''level1'' OSCtable with tabix: <code>tabix -s2 -b3 -e4 OSCTABLE</code>

* Beware that <code>groupBy -o sum</code> prints <code>0</code> if there is no input.
<pre>
printf '' | groupBy -g 1 -c 2 -o sum
0
</pre>

* Regular expression to reduce a library name to its LS-A identifier: <code>s/\S+(CNhs\d{5})\S+/$1/g</code>

* Painting: association between ''non-promoter'' tags and sequence motifs ? (fantom5:00894)

* Very basic data de-duplication:

<pre>
for FILE in *
do
unset SUM
echo -ne "$FILE\t"
SUM=$(sha256sum $FILE | cut -f1 -d' ')
echo $SUM
[ -e ../sha256/$SUM ] || ln $FILE ../sha256/$SUM
[ -e ../sha256/$SUM ] && ln -f ../sha256/$SUM $FILE
done
</pre>

* Cancer-specific promoter chr10:129703392..129703416

*Normalisation factors

Download data:

<pre>
wget https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz
</pre>

Load it in R, edited the column names, and exported its two first lines.

<pre>
osc <- read.table("tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz", row.names=1, head=TRUE)

colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))

write.table(file='FREEZE_PHASE1.normalisation.tsv', t(osc[1:2,]), quote=F, sep="\t")
</pre>

User:Plessy

2013-01-15T07:24:34Z

Plessy: Added FREEZE_PHASE2 files.

= Charles Plessy =

= Works on =

* [[nanoCAGE]]
* [[CAGEscan]]
* [[Olfactory receptors]]
* [[OP-SOLEXA-5RACE-v1.0|Deep-RACE]]
* [[Extend rat gene models with CAGEscan]]

== [[Painting]] and other strange CAGE signals ==

Jessica pointed out that in Fantom5, some hemoglobin genes single out by the absence of locus painting.

https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=d1yxe9xqKJDltVi86B6CFD;loc=mm9::chr11:32168803..32199353

Question: when Timo aligns promoter sequences in which errors were introduced, some reads map on introns and exons. Do we find more of them in intron and exon-painted genes ?

How is it on polysome fractions ?

Need CAGE libraries without cap trapping.

A page is created to stimulate collaboration, [[Painting]].

Is the quantity of tags aligning to the [[rDNA]] correlated to the expression of the [[Ribosome|ribosomal proteins]] ?

= Shared files =

== Whole-genome BED files ==

These BED files recapitulate the expression of all libraries for each genome, and are available on the [[Dataset_introduction#How_to_share_my_analysis_results_with_the_consortium_member.3F|webdav]] for the [[File_release|file releases]]
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1.1/per-genome FREEZE_PHASE1.1] and
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE2/per-genome FREEZE_PHASE2], and previous updates
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_11/per-genome 011],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/per-genome 012],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_13/per-genome 013],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_014/per-genome 014]
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_015/per-genome 015], and
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_016/per-genome 016]. They are sorted by position, [http://samtools.sourceforge.net/tabix.shtml compressed with bgzip and indexed with tabix]. In CTSS files, column 4 contains the library ID, and the genome name in TSS files. For example, here is the result of the command <code>tabix rn4.bed.gz chr1:1923214-1963638</code>.

<pre>
chr1 1923214 1923215 CNhs11312 1 +
chr1 1923717 1923718 CNhs12070 1 +
chr1 1931875 1931876 CNhs10614 1 -
chr1 1963635 1963636 CNhs12070 1 -
chr1 1963636 1963637 CNhs11298 1 -
chr1 1963636 1963637 CNhs12070 3 -
chr1 1963637 1963638 CNhs11298 1 -
chr1 1963637 1963638 CNhs12070 2 -
</pre>

It indicates for instance, that in the library CNhs12070, there was 1 tag on chromosome 1 at position 1963635, 3 at position 1963636 and 2 at position 1963637. Conversely, it indicates that at position 1963636 on chromosome 1, there are 3 tags in library CNhs12070 and 2 tag in CNhs11298.

These files can take a particular advantage of the <code>groupby</code> command from the [https://github.com/arq5x/filo FILO] package:

<pre>
tabix rn4.bed.gz chr1:1923214-1963638 | sort -k4 | groupBy -g 4 -c 5 -ops sum
CNhs10614 1
CNhs11298 2
CNhs11312 1
CNhs12070 7
</pre>

=== Making of ===

==== Interrogating the list of all the HeliScopeCAGE libraries ====

The CNhs_lib function will be needed below.

This needs a FANTOM5 mirror or a flat directory with all [[Dataset_introduction#Where_can_I_get_all_the_details_of_the_libraries.3F|SDRF]] files. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*sdrf.txt</code> is useful for that purpose. The following [http://en.wikipedia.org/wiki/Bash_(Unix_shell) bash] function interrogates them, looking by default in <code>~/public_html/F5/sdrf</code>, but this can be overridden by setting the <code>F5</code> variable.

<pre>
# Return all libraries.
function SDRF_lib {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f14 | grep CNhs
done
}

# Returns only HeliScopeCAGE libraries.
function CNhs_lib {
for CNhs in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $CNhs | cut -f14 | grep CNhs
done
}

# Returns the description of a library.
function CNhs-desc {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f3
done
}
</pre>

==== Flat directory with modified BED CTSS files. ====

Original BED CTSS files are expected to be downloaded in a local directory called <code>longnames</code>. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*bed.gz</code> is useful for that purpose.

The following commands produces modified BED CTSS files, where the <code>name</code> field contain the library name instead of the coordinates. Beware that this format does not allow for the co-existence of more than one alignment file per library. In what follows, the galGal3 alignments are discarded and the galGal4 alignments are kept.

<pre>
# Make sure there are no duplicata.
ls longnames | cut -f2 -d '.' | sort | uniq -d

mkdir -p bgrezip

for LIB in $(ls longnames | cut -f2 -d '.')
do
echo -ne "${LIB}.ctss.bed.gz\t"
zcat longnames/*${LIB}*ctss.bed.gz |
awk -v LIB=$LIB '{OFS="\t"} {$4=LIB ; print}' |
sort -k1,1 -k2,2n |
bgzip |
tee bgrezip/${LIB}.ctss.bed.gz |
md5sum |
cut -f1 -d' '
done |
tee bgrezip.txt

# Index with tabix.
(cd bgrezip;
for BED in *ctss.bed.gz
do
tabix -p bed $BED
done)
</pre>

''Note that wikimedia inserts a strange non-breakable character in the <code>awk</code> command. Paste it from the source.''

Collect the files in the <code>bgrezip</code> directory to replace the original ones.

Note: perhaps <code>-k6,6</code> would make the sorting completely deterministic, but this may not be needed if the input already guarantees that for bases where both strands are transcribed, it is always either the <code>+</code> or <code>-</code> strand that comes first.

==== Construction of whole-genome BED CTSS files ====

Now that the <code>name</code> field contains the library name, the data can pooled in one file per genome.

<pre>
# Tip: run this with a low IO priority:
# exec ionice -c 3 bash

# Mouse
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'
ORGANISM='Mus musculus'
OUTPUT=mm9.ctss.bed.gz

rm -rf chrtmp
mkdir chrtmp

for CHR in $CHRLIST
do
for BED in $( CNhs_lib "$ORGANISM" )
do
if [ -e ${BED}.ctss.bed.gz ]
then
echo -e "${CHR}\t${BED}" 1>&2
tabix ${BED}.ctss.bed.gz ${CHR}
# zgrep -P "${CHR}\t" ${BED}.ctss.bed.gz # Delete after testing with the above, which is almost 10 times faster.
fi
done |
sort --temporary-directory=$(pwd) -k2,2n > chrtmp/${CHR} # Make sure there is enough space for sort's temporary files !
done
cd chrtmp
cat $CHRLIST |
bgzip > $OUTPUT
</pre>

<pre>
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrM chrX chrY'
ORGANISM='Homo sapiens'
OUTPUT=hg19.ctss.bed.gz
</pre>

<pre>
for LIB in $(CNhs_lib 'Canis lupus familiaris'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > canFam2.ctss.bed.gz
for LIB in $(CNhs_lib gallus); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > galGal3.ctss.bed.gz
for LIB in $(CNhs_lib 'Rattus'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > rn4.ctss.bed.gz
</pre>

==== Construction of whole-genome BED TSS files ====

Pooling all libraries reduces the size compressed files.

<pre>
for GENOME in hg19 mm9 rn4 canFam2 galGal3
do
zcat $GENOME.ctss.bed.gz |
groupBy -g 1,2,3,6 -c 5 -o sum |
awk -v GENOME=$GENOME '{OFS="\t"} {print $1, $2, $3, GENOME, $5, $4}' |
bgzip > $GENOME.tss.bed.gz
done
</pre>

==== In Zenbu ====

The whole-genome BED TSS files for human and mouse are loaded in [[ZENBU_Howto|Zenbu]]. Because the library information is collapsed, it is possible to display large intervals quiclkly, or even whole chromosome (~30 s for human chromosome 19).

== Indexed expression files ==

Source:

* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/110804-dpi-clusters-expression/00readme.html (UPDATE_012)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/00readme.html (FREEZE_PHASE1)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/description120126/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz

Four columns were added for sorting and indexing: ''chromosome'', ''start'', ''end'' and ''strand''.

The files for FREEZE_PHASE1 are available at https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/expression and was indexed with [http://samtools.sourceforge.net/tabix.shtml tabix] (see also [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/expression UPDATE_012]).

=== Making of ===
<pre>
TABLE=tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrM chrX chrY'
zgrep '^##' -A3 $TABLE > $(basename $TABLE .osc.txt.gz).header
# Edit the header to add four columns, named ‘chrom start end strand’.
for chrom in $CHRLIST
do
zgrep "^$chrom:" $TABLE |
perl -ne '/^(chr[\d\w]+):(\d+)..(\d+),([+-])/; print join "\t", $1, $2, $3, $4, $_' |
sort --field-separator "$(printf "\t")" -k2,2n -k3,3n -k4,4 > $chrom
done
cat $(basename $TABLE .osc.txt.gz).header $CHRLIST | bgzip > $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz
tabix $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz -s 1 -b 2 -e 3 -S 3 -0
</pre>

For mouse use the following chromosome list instead.
<pre>
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'
</pre>

= Draft notes / Brainstorm =

Assemble chrM and compare recorded ethnicity with mitochondrial haploypes ?

Ultraconserved elements: http://users.soe.ucsc.edu/~jill/ultra.html

* Liftover: http://genome.ucsc.edu/cgi-bin/hgLiftOver
* Or map the fasta files ?
* Or recalculate within hg19 ?

Do we see chromosomal aberrations in the CAGEscan data?

Search for CDR3 regions in the longest reads.

Can we basecall SNPs from the Helicos data ?

Sense-antisense pair at VIM ? HG19::chr10:17270611-17272991

Orphan promoter: HG19::chr11:27908185-27913031

Another orphan promoter, expressed only in nervous tissue: HG19::chr11:28891298-29036678

Long chain: HG19::chr11:31253864-31864905

Cancer cells seem to express more snoRNAs in the ''RPS3'' locus than other cells; it may be interesting to test if they express more ribosomal components in general, perhaps because they are fast-growing.

* hg19::chr15:66816186..66817714 : alternative RPL4 promoter ?

* Strong painting in ''IQSEC1''. https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:12915663..13276777
** And in ''NUP210''.

== Reads and CTSS distributions in arbitrary windows ==

[[File:plessy-Antisense-RefSeqTSS-CTSS-w2000.png|thumb|See the thread fantom5:00740 for more details]]

<pre>
PER_GENOME=/home/charles/public_html/F5/tabix_12
HALF_WINDOW=5000
WINDOW=$(( $HALF_WINDOW * 2 ))

slopBed -i refGeneTss.bed -g /usr/share/bedtools/genomes/human.hg19.genome -l $HALF_WINDOW -r $HALF_WINDOW |
while read chr start end name score strand
do
tabix $PER_GENOME/hg19.bed.gz ${chr}:$(( $start + 1))-${end} |
grep -v "${strand}" |
awk -v start=$start -v end=$end -v strand=$strand '{OFS="\t"} {if (strand == "+") print $2 - start, $5 ; else print end - $3, $5 }' |
groupBy -c 2 -g 1 -o sum -inheader |
perl -ne 'print unless /^0$/'
done > window-$WINDOW.txt

sort -nk1,1 window-$WINDOW.txt | groupBy -g 1 -c 2 -o sum > window-$WINDOW.profile

cut -f1 window-$WINDOW.txt | sort -n | uniq -c | sort -nk2 > window-$WINDOW.mapcount
</pre>

== [[CAGEscan]] ==

* ''MLT1A'' repeat bridging with ''EYPC'': https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=n3gh7Y8TOJGkOqG_q9EZB;loc=hg19::chr12:91344715..91439984

* Why did we lose this promoter in ''IQSEC1''? https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:13007279..13012828

== [[User:Plessy/Richness|Richness]] ==

[[File:Richness_FREEZE1_TBP.png|thumb|TBP]][[File:Richness_FREEZE1_SP1.png|thumb|SP1]]

In the <code>R</code> package called <code>vegan</code> there is a function, called <code>rarefy</code>, to reduce the data as if it had contained a fixed number of tags. Used on the libraries, it evaluates how complex they are (genome-bioinfo:02933). Another function, <code>rrarefy</code> makes a random sub-sample of a library with an arbitrary number of counts, respecting the original distribution. See [[User:Plessy/Richness]] for more details.

I calculated the ''richness'' of the robust clusters in FREEZE1 after normalising the libraries to one million tags by random rarefaction. The richness evaluates how ubiquitous the clusters are. For instance, with a sampling size of 100, the clusters with a TBP motif have a richness of 55.8 ± 29.2 (''n''=948) and the clusters with a SP1 motif have a richness of 76.4 ± 18.4 (''n''=16,134), a significant difference (''p'' < 2.2e-16).

The profiles for each motif are available at the following URL: [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/ https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/], as well as the summary panels [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness//T1.png T1.png] to [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/T8.png T8.png]. The sub-sampled data used for this analysis is available in the file [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2 tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2]

=== To Do ===

* Cross data with Mette's motifs and Sebastian's houskeeping list <code>fantom5:01568</code>
* Filter CAGE peak robust set with predicted TSS. <code>fantom5:01581</code>

= Misc =

* [[User:Plessy/Debian|Software used in FANTOM5 and available in the Debian operating system.]]

* Sorting a ''level1'' OscTable called OSCTABLE: <code>cat <(grep \# -A1 OSCTABLE) <(grep -v \# OSCTABLE | sed '1d' | sort --field-separator "$(printf "\t")" -k2.4,2n -k 2.4,2.4 -k3,3n -k4,4n -k5,5) | sponge OSCTABLE</code>

* Indexing a sorted and bgzipped ''level1'' OSCtable with tabix: <code>tabix -s2 -b3 -e4 OSCTABLE</code>

* Beware that <code>groupBy -o sum</code> prints <code>0</code> if there is no input.
<pre>
printf '' | groupBy -g 1 -c 2 -o sum
0
</pre>

* Regular expression to reduce a library name to its LS-A identifier: <code>s/\S+(CNhs\d{5})\S+/$1/g</code>

* Painting: association between ''non-promoter'' tags and sequence motifs ? (fantom5:00894)

* Very basic data de-duplication:

<pre>
for FILE in *
do
unset SUM
echo -ne "$FILE\t"
SUM=$(sha256sum $FILE | cut -f1 -d' ')
echo $SUM
[ -e ../sha256/$SUM ] || ln $FILE ../sha256/$SUM
[ -e ../sha256/$SUM ] && ln -f ../sha256/$SUM $FILE
done
</pre>

* Cancer-specific promoter chr10:129703392..129703416

*Normalisation factors

Download data:

<pre>
wget https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz
</pre>

Load it in R, edited the column names, and exported its two first lines.

<pre>
osc <- read.table("tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz", row.names=1, head=TRUE)

colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))

write.table(file='FREEZE_PHASE1.normalisation.tsv', t(osc[1:2,]), quote=F, sep="\t")
</pre>

User:Plessy

2013-01-11T10:05:42Z

Plessy: Updated to FREEZE_PHASE1.1

= Charles Plessy =

= Works on =

* [[nanoCAGE]]
* [[CAGEscan]]
* [[Olfactory receptors]]
* [[OP-SOLEXA-5RACE-v1.0|Deep-RACE]]
* [[Extend rat gene models with CAGEscan]]

== [[Painting]] and other strange CAGE signals ==

Jessica pointed out that in Fantom5, some hemoglobin genes single out by the absence of locus painting.

https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=d1yxe9xqKJDltVi86B6CFD;loc=mm9::chr11:32168803..32199353

Question: when Timo aligns promoter sequences in which errors were introduced, some reads map on introns and exons. Do we find more of them in intron and exon-painted genes ?

How is it on polysome fractions ?

Need CAGE libraries without cap trapping.

A page is created to stimulate collaboration, [[Painting]].

Is the quantity of tags aligning to the [[rDNA]] correlated to the expression of the [[Ribosome|ribosomal proteins]] ?

= Shared files =

== Whole-genome BED files ==

These BED files recapitulate the expression of all libraries for each genome, and are available on the [[Dataset_introduction#How_to_share_my_analysis_results_with_the_consortium_member.3F|webdav]] for the [[File_release|file releases]]
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1.1/per-genome FREEZE_PHASE1.1],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_016/per-genome UPDATE_016], and previous updates
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_11/per-genome 011],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/per-genome 012],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_13/per-genome 013],
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_014/per-genome 014], and
[https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_015/per-genome 015]. They are sorted by position, [http://samtools.sourceforge.net/tabix.shtml compressed with bgzip and indexed with tabix]. In CTSS files, column 4 contains the library ID, and the genome name in TSS files. For example, here is the result of the command <code>tabix rn4.bed.gz chr1:1923214-1963638</code>.

<pre>
chr1 1923214 1923215 CNhs11312 1 +
chr1 1923717 1923718 CNhs12070 1 +
chr1 1931875 1931876 CNhs10614 1 -
chr1 1963635 1963636 CNhs12070 1 -
chr1 1963636 1963637 CNhs11298 1 -
chr1 1963636 1963637 CNhs12070 3 -
chr1 1963637 1963638 CNhs11298 1 -
chr1 1963637 1963638 CNhs12070 2 -
</pre>

It indicates for instance, that in the library CNhs12070, there was 1 tag on chromosome 1 at position 1963635, 3 at position 1963636 and 2 at position 1963637. Conversely, it indicates that at position 1963636 on chromosome 1, there are 3 tags in library CNhs12070 and 2 tag in CNhs11298.

These files can take a particular advantage of the <code>groupby</code> command from the [https://github.com/arq5x/filo FILO] package:

<pre>
tabix rn4.bed.gz chr1:1923214-1963638 | sort -k4 | groupBy -g 4 -c 5 -ops sum
CNhs10614 1
CNhs11298 2
CNhs11312 1
CNhs12070 7
</pre>

=== Making of ===

==== Interrogating the list of all the HeliScopeCAGE libraries ====

The CNhs_lib function will be needed below.

This needs a FANTOM5 mirror or a flat directory with all [[Dataset_introduction#Where_can_I_get_all_the_details_of_the_libraries.3F|SDRF]] files. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*sdrf.txt</code> is useful for that purpose. The following [http://en.wikipedia.org/wiki/Bash_(Unix_shell) bash] function interrogates them, looking by default in <code>~/public_html/F5/sdrf</code>, but this can be overridden by setting the <code>F5</code> variable.

<pre>
# Return all libraries.
function SDRF_lib {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f14 | grep CNhs
done
}

# Returns only HeliScopeCAGE libraries.
function CNhs_lib {
for CNhs in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $CNhs | cut -f14 | grep CNhs
done
}

# Returns the description of a library.
function CNhs-desc {
for SDRF in $(find ${F5-$HOME/public_html/F5/sdrf} -name '*sdrf.txt')
do
grep "$1" $SDRF | cut -f3
done
}
</pre>

==== Flat directory with modified BED CTSS files. ====

Original BED CTSS files are expected to be downloaded in a local directory called <code>longnames</code>. The [[Dataset_introduction#How_to_download_a_series_of_files.3F|lftp]] command <code>mget ./f5pipeline/*/*bed.gz</code> is useful for that purpose.

The following commands produces modified BED CTSS files, where the <code>name</code> field contain the library name instead of the coordinates. Beware that this format does not allow for the co-existence of more than one alignment file per library. In what follows, the galGal3 alignments are discarded and the galGal4 alignments are kept.

<pre>
# Make sure there are no duplicata.
ls longnames | cut -f2 -d '.' | sort | uniq -d

mkdir -p bgrezip

for LIB in $(ls longnames | cut -f2 -d '.')
do
echo -ne "${LIB}.ctss.bed.gz\t"
zcat longnames/*${LIB}*ctss.bed.gz |
awk -v LIB=$LIB '{OFS="\t"} {$4=LIB ; print}' |
sort -k1,1 -k2,2n |
bgzip |
tee bgrezip/${LIB}.ctss.bed.gz |
md5sum |
cut -f1 -d' '
done |
tee bgrezip.txt

# Index with tabix.
(cd bgrezip;
for BED in *ctss.bed.gz
do
tabix -p bed $BED
done)
</pre>

''Note that wikimedia inserts a strange non-breakable character in the <code>awk</code> command. Paste it from the source.''

Collect the files in the <code>bgrezip</code> directory to replace the original ones.

Note: perhaps <code>-k6,6</code> would make the sorting completely deterministic, but this may not be needed if the input already guarantees that for bases where both strands are transcribed, it is always either the <code>+</code> or <code>-</code> strand that comes first.

==== Construction of whole-genome BED CTSS files ====

Now that the <code>name</code> field contains the library name, the data can pooled in one file per genome.

<pre>
# Tip: run this with a low IO priority:
# exec ionice -c 3 bash

# Mouse
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'
ORGANISM='Mus musculus'
OUTPUT=mm9.ctss.bed.gz

rm -rf chrtmp
mkdir chrtmp

for CHR in $CHRLIST
do
for BED in $( CNhs_lib "$ORGANISM" )
do
if [ -e ${BED}.ctss.bed.gz ]
then
echo -e "${CHR}\t${BED}" 1>&2
tabix ${BED}.ctss.bed.gz ${CHR}
# zgrep -P "${CHR}\t" ${BED}.ctss.bed.gz # Delete after testing with the above, which is almost 10 times faster.
fi
done |
sort --temporary-directory=$(pwd) -k2,2n > chrtmp/${CHR} # Make sure there is enough space for sort's temporary files !
done
cd chrtmp
cat $CHRLIST |
bgzip > $OUTPUT
</pre>

<pre>
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrM chrX chrY'
ORGANISM='Homo sapiens'
OUTPUT=hg19.ctss.bed.gz
</pre>

<pre>
for LIB in $(CNhs_lib 'Canis lupus familiaris'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > canFam2.ctss.bed.gz
for LIB in $(CNhs_lib gallus); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > galGal3.ctss.bed.gz
for LIB in $(CNhs_lib 'Rattus'); do ls $LIB.ctss.bed.gz; done | xargs zcat | sort -k1,1 -k2,2n | bgzip > rn4.ctss.bed.gz
</pre>

==== Construction of whole-genome BED TSS files ====

Pooling all libraries reduces the size compressed files.

<pre>
for GENOME in hg19 mm9 rn4 canFam2 galGal3
do
zcat $GENOME.ctss.bed.gz |
groupBy -g 1,2,3,6 -c 5 -o sum |
awk -v GENOME=$GENOME '{OFS="\t"} {print $1, $2, $3, GENOME, $5, $4}' |
bgzip > $GENOME.tss.bed.gz
done
</pre>

==== In Zenbu ====

The whole-genome BED TSS files for human and mouse are loaded in [[ZENBU_Howto|Zenbu]]. Because the library information is collapsed, it is possible to display large intervals quiclkly, or even whole chromosome (~30 s for human chromosome 19).

== Indexed expression files ==

Source:

* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/110804-dpi-clusters-expression/00readme.html (UPDATE_012)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/00readme.html (FREEZE_PHASE1)
* https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/description120126/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz

Four columns were added for sorting and indexing: ''chromosome'', ''start'', ''end'' and ''strand''.

The files for FREEZE_PHASE1 are available at https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/expression and was indexed with [http://samtools.sourceforge.net/tabix.shtml tabix] (see also [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/UPDATE_12/expression UPDATE_012]).

=== Making of ===
<pre>
TABLE=tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.desc.osc.txt.gz
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrM chrX chrY'
zgrep '^##' -A3 $TABLE > $(basename $TABLE .osc.txt.gz).header
# Edit the header to add four columns, named ‘chrom start end strand’.
for chrom in $CHRLIST
do
zgrep "^$chrom:" $TABLE |
perl -ne '/^(chr[\d\w]+):(\d+)..(\d+),([+-])/; print join "\t", $1, $2, $3, $4, $_' |
sort --field-separator "$(printf "\t")" -k2,2n -k3,3n -k4,4 > $chrom
done
cat $(basename $TABLE .osc.txt.gz).header $CHRLIST | bgzip > $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz
tabix $(basename $TABLE .osc.txt.gz).sorted.osc.txt.gz -s 1 -b 2 -e 3 -S 3 -0
</pre>

For mouse use the following chromosome list instead.
<pre>
CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'
</pre>

= Draft notes / Brainstorm =

Assemble chrM and compare recorded ethnicity with mitochondrial haploypes ?

Ultraconserved elements: http://users.soe.ucsc.edu/~jill/ultra.html

* Liftover: http://genome.ucsc.edu/cgi-bin/hgLiftOver
* Or map the fasta files ?
* Or recalculate within hg19 ?

Do we see chromosomal aberrations in the CAGEscan data?

Search for CDR3 regions in the longest reads.

Can we basecall SNPs from the Helicos data ?

Sense-antisense pair at VIM ? HG19::chr10:17270611-17272991

Orphan promoter: HG19::chr11:27908185-27913031

Another orphan promoter, expressed only in nervous tissue: HG19::chr11:28891298-29036678

Long chain: HG19::chr11:31253864-31864905

Cancer cells seem to express more snoRNAs in the ''RPS3'' locus than other cells; it may be interesting to test if they express more ribosomal components in general, perhaps because they are fast-growing.

* hg19::chr15:66816186..66817714 : alternative RPL4 promoter ?

* Strong painting in ''IQSEC1''. https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:12915663..13276777
** And in ''NUP210''.

== Reads and CTSS distributions in arbitrary windows ==

[[File:plessy-Antisense-RefSeqTSS-CTSS-w2000.png|thumb|See the thread fantom5:00740 for more details]]

<pre>
PER_GENOME=/home/charles/public_html/F5/tabix_12
HALF_WINDOW=5000
WINDOW=$(( $HALF_WINDOW * 2 ))

slopBed -i refGeneTss.bed -g /usr/share/bedtools/genomes/human.hg19.genome -l $HALF_WINDOW -r $HALF_WINDOW |
while read chr start end name score strand
do
tabix $PER_GENOME/hg19.bed.gz ${chr}:$(( $start + 1))-${end} |
grep -v "${strand}" |
awk -v start=$start -v end=$end -v strand=$strand '{OFS="\t"} {if (strand == "+") print $2 - start, $5 ; else print end - $3, $5 }' |
groupBy -c 2 -g 1 -o sum -inheader |
perl -ne 'print unless /^0$/'
done > window-$WINDOW.txt

sort -nk1,1 window-$WINDOW.txt | groupBy -g 1 -c 2 -o sum > window-$WINDOW.profile

cut -f1 window-$WINDOW.txt | sort -n | uniq -c | sort -nk2 > window-$WINDOW.mapcount
</pre>

== [[CAGEscan]] ==

* ''MLT1A'' repeat bridging with ''EYPC'': https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=n3gh7Y8TOJGkOqG_q9EZB;loc=hg19::chr12:91344715..91439984

* Why did we lose this promoter in ''IQSEC1''? https://fantom5-collaboration.gsc.riken.jp/zenbu/gLyphs/#config=zoZVMHT9D9ysb8G7Z84TED;loc=hg19::chr3:13007279..13012828

== [[User:Plessy/Richness|Richness]] ==

[[File:Richness_FREEZE1_TBP.png|thumb|TBP]][[File:Richness_FREEZE1_SP1.png|thumb|SP1]]

In the <code>R</code> package called <code>vegan</code> there is a function, called <code>rarefy</code>, to reduce the data as if it had contained a fixed number of tags. Used on the libraries, it evaluates how complex they are (genome-bioinfo:02933). Another function, <code>rrarefy</code> makes a random sub-sample of a library with an arbitrary number of counts, respecting the original distribution. See [[User:Plessy/Richness]] for more details.

I calculated the ''richness'' of the robust clusters in FREEZE1 after normalising the libraries to one million tags by random rarefaction. The richness evaluates how ubiquitous the clusters are. For instance, with a sampling size of 100, the clusters with a TBP motif have a richness of 55.8 ± 29.2 (''n''=948) and the clusters with a SP1 motif have a richness of 76.4 ± 18.4 (''n''=16,134), a significant difference (''p'' < 2.2e-16).

The profiles for each motif are available at the following URL: [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/ https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/per-TF/], as well as the summary panels [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness//T1.png T1.png] to [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/T8.png T8.png]. The sub-sampled data used for this analysis is available in the file [https://fantom5-collaboration.gsc.riken.jp/webdav/home/plessy/FREEZE_PHASE1/richness/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2 tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.counts.selected.sync015.clustername_update.osc.r1000000.1.txt.bz2]

=== To Do ===

* Cross data with Mette's motifs and Sebastian's houskeeping list <code>fantom5:01568</code>
* Filter CAGE peak robust set with predicted TSS. <code>fantom5:01581</code>

= Misc =

* [[User:Plessy/Debian|Software used in FANTOM5 and available in the Debian operating system.]]

* Sorting a ''level1'' OscTable called OSCTABLE: <code>cat <(grep \# -A1 OSCTABLE) <(grep -v \# OSCTABLE | sed '1d' | sort --field-separator "$(printf "\t")" -k2.4,2n -k 2.4,2.4 -k3,3n -k4,4n -k5,5) | sponge OSCTABLE</code>

* Indexing a sorted and bgzipped ''level1'' OSCtable with tabix: <code>tabix -s2 -b3 -e4 OSCTABLE</code>

* Beware that <code>groupBy -o sum</code> prints <code>0</code> if there is no input.
<pre>
printf '' | groupBy -g 1 -c 2 -o sum
0
</pre>

* Regular expression to reduce a library name to its LS-A identifier: <code>s/\S+(CNhs\d{5})\S+/$1/g</code>

* Painting: association between ''non-promoter'' tags and sequence motifs ? (fantom5:00894)

* Very basic data de-duplication:

<pre>
for FILE in *
do
unset SUM
echo -ne "$FILE\t"
SUM=$(sha256sum $FILE | cut -f1 -d' ')
echo $SUM
[ -e ../sha256/$SUM ] || ln $FILE ../sha256/$SUM
[ -e ../sha256/$SUM ] && ln -f ../sha256/$SUM $FILE
done
</pre>

* Cancer-specific promoter chr10:129703392..129703416

*Normalisation factors

Download data:

<pre>
wget https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/111220-DPI/hg19/tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz
</pre>

Load it in R, edited the column names, and exported its two first lines.

<pre>
osc <- read.table("tc.decompose_smoothing_merged.ctssMaxCounts11_ctssMaxTpm1.tpm.selected.clustername_update.osc.txt.gz", row.names=1, head=TRUE)

colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))

write.table(file='FREEZE_PHASE1.normalisation.tsv', t(osc[1:2,]), quote=F, sep="\t")
</pre>