...
Expand | |||||
---|---|---|---|---|---|
| |||||
|
Code Block | ||||
---|---|---|---|---|
| ||||
cd $SCRATCH/core_ngs/bedtools_multicov bedtools multicov -s -bams yeast_mrna.sort.filt.bam \ -bed merged.good.sc_genes.bed > yeast_mrna_gene_counts.bed |
...
A signal track is a bedGraph (BED3+) file with an entry for every base in a defined set of regions that shows the count of overlapping bases for the regions (see https://genome.ucsc.edu/goldenpath/help/bedgraph.html). bedGraph files can be visualized in the Broad's IGV (Integrative Genomics Viewer) application (https://software.broadinstitute.org/software/igv/download) or in the UCSC Genome Browser (https://genome.ucsc.edu/).The bedtools genomecov function (
- Go to the UCSC Genome Browser: https://
...
...
...
- edu/
- Select Genomes from the top menu bar
- Select Human from POPULAR SPECIES
- under Human Assembly select Feb 2009 (GrCh37/hg19)
- select GO
- In the hg19 browser page, the Layered H3K27Ac track is a signal track
- the x-axis is the genome position
- the y-axis represents the count of ChIP-seq reads that overlap each position
- where the ChIP'd protein is H3K27AC (histone H3, acetylated on the Lysine at amino acid position 27)
The bedtools genomecov function (https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html), with the -bg (bedgraph) option produces output in bedGraph format. Here we'll analyze the per-base coverage of yeast RNAseq reads in our merged yeast gene regions.
...
Code Block | ||||
---|---|---|---|---|
| ||||
idev -m 120 -N 1 -A OTH21164 -r CoreNGSday5CoreNGS-day5 module# loador biocontainersidev module-m load90 bedtools -N 1 -A OTH21164 -p development module load biocontainers module load bedtools mkdir -p $SCRATCH/core_ngs/bedtools_genomecov cd $SCRATCH/core_ngs/bedtools_genomecov cp $CORENGS/catchup/bedtools_merge/merged*bed $SCRATCH/core_ngs/bedtools_genomecov/ . cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* $SCRATCH/core_ngs/bedtools_genomecov/. |
Then calling bedtools genomecov is easy. The -bg option says to report the depth in bedGraph format.
...
Code Block | ||
---|---|---|
| ||
cd $SCRATCH/core_ngs/bedtools_genomecov export LC_COLLATE=C; sort -k1,1 -k2,2n -k3,3n yeast_mrna.genomecov.bedGraph > yeast_mrna.genomecov.sorted.bedGraph bedGraphToBigWig yeast_mrna.genomecov.sorted.bedGraph sc_chrom_sizes.txt yeast_mrna.genomecov.bw |
Just like zcat converts gzip'd files to text, and samtools view convets binary BAM files to text, the bigWigToBedGraph program can convert binary bigWig format to text. Unfortunately, the ucsc-bigwigtobedgraph BioContainer seems to be broken, so we'll use a version in the BioITeam area insteadSee the size difference between the bedGraph and the bigWig files. The bigWig (9.7M) is less that 1/3 the size of the bedGraph (34M).
Code Block | ||
---|---|---|
| ||
cd $SCRATCH/core_ngs/bedtools_genomecov
ls -lh yeast_mrna.genome*
|
Since the bigWig file is binary, not text, you can't use commands like cat, head, tail on them directly and get meaningful output. Instead, just as zcat converts gzip'd files to text, and samtools view convets binary BAM files to text, the bigWigToBedGraph program can convert binary bigWig format to text. That's a different BioContainers module (ucsc-bigwigtobedgraph) and the default container version doesn't work, so we'll specifically load one that does:
Code Block | ||
---|---|---|
| ||
cd $SCRATCH/core_ngs/bedtools_genomecov# The default version of is broken, so load this specific biocontainers version module load ucsc-bigwigtobedgraph/ctr-357--1 # see usage for bigWigToBedGraph: bigWigToBedGraph cd $SCRATCH/work/projects/BioITeam/common/opt/UCSC_utils.2019_08/bigWigToBedGraph core_ngs/bedtools_genomecov # use the program to view a few lines of the binary bigWig file /work/projects/BioITeam/common/opt/UCSC_utils.2019_08/bigWigToBedGraph \ yeast_mrna.genomecov.bw stdout | head |
...