Versions Compared

Key

  • This line was added.
  • This line was removed.
  • Formatting was changed.

...

Expand
titleSetup (if needed)


Code Block
languagebash
idev -m 120 -N 1 -A OTH21164 -r CoreNGSday5
module load biocontainers
module load samtools
module load bedtools

mkdir -p $SCRATCH/core_ngs/bedtools_multicov
cd $SCRATCH/core_ngs/bedtools_multicov
cp $CORENGS/catchup/bedtools_merge/merged*bed      $SCRATCH/core_ngs/bedtools_multicov/.
cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* $SCRATCH/core_ngs/bedtools_multicov/.



Code Block
languagebash
titleRun bedtools multicov to count BAM alignments overlapping a set of genes
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 (

...

...

...

  • 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
languagebash
titlePrepare for bedtools coverage
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
languagebash
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
languagebash
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
languagebash
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

...