I have a mapped CAGE file (hg19), which gives CTSS positions along with the digital count of reads supporting each CTSS. There are 429428 reads.
I need to retrieve the rRNA with UCSC, so as to quantify the CTSS in the rRNA genomic regions.
I have created a table of 'rmsk', but am not sure what kind of filter I need to apply and what other parameters to set? How do I know how many rRNA regions to expect?
Thank you
1 answer
You can download the rRNA data into a sorted BED-like file:
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/rmsk.txt.gz \
| gunzip -c - \
| grep rRNA - \
| cut -f6- - \
| awk '{print $1"\t"$2"\t"$3"\t"$6"\t"$7"\t"$5}' - \
| sort-bed - \
> rRNA.bed
Let's say you have an unsorted BED file of CTSS positions, which you have prepared as a sorted BED file:
$ sort-bed unsorted.CTSS.bed > CTSS.bed
Once you have these two files, you can use BEDOPS bedmap --count to count the number of CTSS positions overlapping the rRNA elements:
$ bedmap --echo --count rRNA.bed CTSS.bed > answer.bed
Each line of the file answer.bed will contain the rRNA element and the number of CTSS positions that overlap that element.
If you want more than just a count of CTSS, the bedmap tool offers a variety of quantitative and qualitative options: the documentation offers much more detail.
•
written
by
Alex Reynolds
Perhaps you could load your CTSS positions and rRNA annotations as custom tracks and then intersect them within the browser.
•
written
by
Alex Reynolds
Which are you trying? The command-line approach or the UCSC browser approach?
•
written
by
Alex Reynolds
I often do as follows:
- using bedtools intersect, create a file where each line indicates which CTSS overlaps with which RNA region. Include dummy entries for CTSS overlapping with nothing. Ensure the file is sorted the same way as the CTSS expression table.
- In R, load the expression table and the file created above (let's call it an overlap table). Transform the CTSS expression table in a RNA region expression table using the
rowsumfunction (do not confuse it withrowSums). Use the overlap table to provide the factor needed byrowsum.
You can see a variant of this strategy in my tutorial under construction on GitHub. Sorry if it does not build from scratch, it may already have bitrot a bit.
•
written
by
Charles Plessy
My data looks as such:
V1 V2 V3 V4 V5 V6
1 chr1 108113121 108113122 chr1:108113121-108113122,- 3 -
2 chr1 108113470 108113471 chr1:108113470-108113471,- 1 -
3 chr1 237766677 237766678 chr1:237766677-237766678,+ 1 +
I tried to add the colnames, but it came up for the first 2 columns and I am not entirely sure which is for CTSS and which the rRNA.
inter_result <- read.table("Intersect_on_CTSS_and_rRNA.bed", stringsAsFactors = FALSE)
colnames(inter_result) <- c("TSS", "Cluster")
head(inter_result)
inter2 <- rowsum(GBGR[c(HeLa, THP1)], linter_result$Cluster)
I am not really sure about what all your files and variables are doing in the tutorial you gave.
Sorry, the tutorial covers a much wider scope; I understand that it is confusing in the context of our discussion. Here is an example using your test with 3 TSSs. I stored them in a file called toto.bed.
cat toto.bed
chr1 108113121 108113122 chr1:108113121-108113122,- 3 -
chr1 108113470 108113471 chr1:108113470-108113471,- 1 -
chr1 237766677 237766678 chr1:237766677-237766678,+ 1 +
Then I downloaded GENCODE23 and extracted position of features such as promoters, exons, genes, etc. (Do not run this command without checking by yourself that there is no malicious code !). It produces a file called "gencode.v23.annotation.bed".
curl --silent https://gist.githubusercontent.com/charles-plessy/9dbc8bc98fb773bf71b6/raw/getAndParseGencode.bash | tee getAndParseGencode.bash | bash
Then, let's intersect the TSSes and the annotation, re-transform the output in a simple association table, and collapse the entries that have the same coordinates.
bedtools intersect -a toto.bed -b gencode.v23.annotation.bed -s -loj |
awk '{OFS="\t"}{print $1":"$2"-"$3$6,$10}' |
bedtools groupby -g 1 -c 2 -o distinct > annot.bed
Here are the results. "." stand for "no annotation".
chr1:108113121-108113122- .
chr1:108113470-108113471- .
chr1:237766677-237766678+ gene,protein_coding_RYR2
Then, in R:
> tss <- read.table('toto.bed', row.names=4)
> annot <- read.table('annot.bed')
> rowsum(tss[,c("V5")], annot$V2)
[,1]
. 4
gene,protein_coding_RYR2 1
The example is a bit contrived, but I hope it gives you a good start.
•
written
by
Charles Plessy
I would just use the command line and BEDOPS bedmap --count:
$ bedmap --echo --count rRNA.bed reads.bed > answer.bed
I guess you could use third-party libraries in R like GenomicRanges, but doing genomic set operations is likely slower this way, since they probably don't use the enhancements put in bedops and bedmap tools to do things efficiently. I wouldn't recommend it.
•
written
by
Alex Reynolds
The big advantage of doing part of the work in R is that it can process a whole expression table at a time, instead of one sample per BED file. I found the rowsum comment to be surprisingly efficient: it works even with millions of entries.
•
written
by
Charles Plessy
Log in to answer this question.
-
Extracting promotes from UCSC
written by peter •
-
twoBitToFa retrieving longer sequence than expected in hg19
written by biosol
-
Get positions from FASTA alignment in bed format
written by kalahari •
-
quality control rRNA mapping
written by ferozef •
-
[bwa sai2sam_pe] fail to locate the index
written by genomics Newbie •
-
Quantification of Maps to rRNA and mtDNA
written by dec986
-
How to query sequences from reference genome with R?
written by amir.jebelli •
-
What galaxy tool should I use to compare these regions?
written by espop23
-
Why Does Cufflinks With Mask (-M) Option Have Lower Fpkm For Mrna Genes?
written by Xianjun
-
Rrna Removal In Rna-Seq Data
written by Nicolas Rosewick
I do not have the answer for the UCSC browser, but it is possible to intersect CAGE data with annotation data in the Zebu browser. There is an extensive example in its documentation. Once such a track has been created, it is possible to download it in various formats, where the expression score will reflect the number of CTSS that intersect with a given feature.