User:Plessy

From Wiki
Jump to navigationJump to search

Charles Plessy

Works on

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 ribosomal proteins ?

Shared files

Whole-genome BED files

These BED files, announced in the message number 01884 on the fantom5 mailing list, recapitulate the expression of all libraries for each genome, and are available on the webdav for the file releases FREEZE_PHASE1.1 (public link) and FREEZE_PHASE2, and previous updates 011, 012, 013, 014 015, and 016. They are sorted by position, 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 tabix rn4.bed.gz chr1:1923214-1963638. 

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       -

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 groupby command from the FILO package: 

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

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 SDRF files. The lftp command mget ./f5pipeline/*/*sdrf.txt is useful for that purpose. The following bash function interrogates them, looking by default in ~/public_html/F5/sdrf, but this can be overridden by setting the F5 variable. 

# 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
}

Selection of libraries

Alternatively to the use of the CNhs_lib function (see above) one can make selections based on the md5sum.txt file distributed with the FANTOM5 data.

For FREEZE_PHASE2. 

HUMAN_hCAGE_LIBS='
  human.cell_line.LQhCAGE
  human.cell_line.hCAGE
  human.fractionation.hCAGE
  human.primary_cell.LQhCAGE
  human.primary_cell.hCAGE
  human.qualitycontrol.hCAGE
  human.timecourse.LQhCAGE
  human.timecourse.hCAGE
  human.tissue.hCAGE'

MOUSE_hCAGE_LIBS='
  mouse.cell_line.hCAGE
  mouse.primary_cell.LQhCAGE
  mouse.primary_cell.hCAGE
  mouse.qualitycontrol.LQhCAGE
  mouse.qualitycontrol.hCAGE
  mouse.timecourse.LQhCAGE
  mouse.timecourse.hCAGE
  mouse.tissue.LQhCAGE
  mouse.tissue.hCAGE'

Flat directory with modified BED CTSS files.

Original BED CTSS files are expected to be downloaded in a local directory called longnames. The lftp command mget ./f5pipeline/*/*bed.gz is useful for that purpose.

The following commands produces modified BED CTSS files, where the name 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. 

# 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 -k6,6 |
    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)

Note that wikimedia inserts a strange non-breakable character in the awk command. Paste it from the source.

Collect the files in the bgrezip directory to replace the original ones.

Construction of whole-genome BED CTSS files

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

# 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
LIBRARY_LIST=$( CNhs_lib "$ORGANISM" )
# or LIBRARY_LIST=$(for DIR in $MOUSE_hCAGE_LIBS ; do grep $DIR  ../md5sum.txt | grep 'ctss.bed.gz' | cut -f 5 -d '.' ; done)
# or LIBRARY_LIST=$(for DIR in $HUMAN_hCAGE_LIBS ; do grep $DIR  ../md5sum.txt | grep 'ctss.bed.gz' | cut -f 5 -d '.' ; done)

rm -rf chrtmp
mkdir chrtmp

for CHR in $CHRLIST
do
  for BED in $LIBRARY_LIST
  do
    if [ -e ${BED}.ctss.bed.gz ]
    then
      echo -e "${CHR}\t${BED}" 1>&2
      tabix  ${BED}.ctss.bed.gz ${CHR}
    fi
  done |
    sort --temporary-directory=$(pwd) -k2,2n -k6,6 > chrtmp/${CHR} # Make sure there is enough space for sort's temporary files !
done
cd chrtmp
cat $CHRLIST |
  bgzip > $OUTPUT

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

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

Construction of whole-genome BED TSS files

Pooling all libraries reduces the size compressed files. 

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

In Zenbu

The whole-genome BED TSS files for human and mouse are loaded in 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).

Quality control

  • Check that all the chromosomes are present in lexical order:
for BED in *tss.bed.gz
do
  zcat $BED | cut -f 1 | uniq -c
done
  • Check that all the expected libraries are present.

(this can generate some discrepancies with some unclassified libraries). 

for BED in *ctss.bed.gz
do
  grep $BED md5sum.txt | cut -f 5 -d '.' | wc -l
  zcat $BED | cut -f 4 | perl -E 'while (<>) {$h{$_}++} ; print keys(%h)' | wc -l
done
  • Check that only one name is found per TSS file
for BED in *.tss.bed.gz
do
  zcat $BED | cut -f 4 | perl -E 'while (<>) {$h{$_}++} ; print keys(%h)' | wc -l
done
  • The clusters in the same base are sorted by strand. Otherwise, the TSS file will have more than one entry per base and per strand. For instance, here is the interval chr1:564597-564598 when sorting was proper. In case of improper sorting, many more lines are produced.
tabix hg19.tss.bed.gz chr1:564597-564598

chr1	564596	564597	hg19	9	+
chr1	564597	564598	hg19	14	-
chr1	564597	564598	hg19	95	+

Indexed expression files

Source:

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 tabix (see also UPDATE_012).

Making of

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

For mouse use the following chromosome list instead. 

CHRLIST='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chrM chrX chrY'

Draft notes / Brainstorm

Assemble chrM and compare recorded ethnicity with mitochondrial haploypes ?

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

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 ?

Reads and CTSS distributions in arbitrary windows

Error creating thumbnail: Unable to run external programs, proc_open() is disabled. Error code: 1
See the thread fantom5:00740 for more details
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

CAGEscan

Richness

Error creating thumbnail: Unable to run external programs, proc_open() is disabled. Error code: 1
TBP
Error creating thumbnail: Unable to run external programs, proc_open() is disabled. Error code: 1
SP1

In the R package called vegan there is a function, called rarefy, 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, rrarefy makes a random sub-sample of a library with an arbitrary number of counts. 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/, as well as the summary panels T1.png to T8.png. The sub-sampled data used for this analysis is available in the file 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 fantom5:01568
  • Filter CAGE peak robust set with predicted TSS. fantom5:01581

Misc

  • Sorting a level1 OscTable called OSCTABLE: 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
  • Indexing a sorted and bgzipped level1 OSCtable with tabix: tabix -s2 -b3 -e4 OSCTABLE
  • Beware that groupBy -o sum prints 0 if there is no input.
printf '' | groupBy -g 1 -c 2 -o sum
0
  • Regular expression to reduce a library name to its LS-A identifier: s/\S+(CNhs\d{5})\S+/$1/g
  • Painting: association between non-promoter tags and sequence motifs ? (fantom5:00894)
  • Very basic data de-duplication:
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
  • Cancer-specific promoter chr10:129703392..129703416
  • Normalisation factors

Download data: 

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

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

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")
  • Normalisation factors for FREEZE_PHASE2
wget https://fantom5-collaboration.gsc.riken.jp/webdav/home/kawaji/121025-DPI-robust_phase1_pls_phase2/hg19/robust_phase1_pls_2.counts.osc.txt.gz
zgrep -m 1 -B1 01STAT:MAPPED robust_phase1_pls_2.counts.osc.txt.gz > robust_phase1_pls_2.counts.mapped.txt
x <- t(read.table('robust_phase1_pls_2.counts.mapped.txt', header=TRUE, row.names=1))
rownames(x) <- regmatches(rownames(x), regexpr('CNhs.....', rownames(x)))
Retrieved from "http://fantom5-collaboration.gsc.riken.jp/wiki/index.php?title=User:Plessy&oldid=6594"