| Tip | ||
|---|---|---|
| ||
Use our today's reservation ( core-ngs-class-0606) when submitting batch jobs to get higher priority on the ls6 normal queue
|
| Table of Contents |
|---|
The BED format
...
.
|
| Table of Contents |
|---|
The BED format
BED (Browser Extensible Data) format is a simple text format for location-oriented data (genomic regions) developed to support UCSC Genome Browser (GenBrowse) tracks. Standard BED files have 3 to 6 Tab-separated columns, although up to 12 columns are defined. (Read more about the UCSC Genome Browser's official BED format.)
...
- chrom (required) – string naming the chromosome or other contig
- start (required) – the 0-based start position of the region
- end (required) – the 1-based end position of the region
- name (optional) – an arbitrary string describing the region
- for BED files loaded as UCSC Genome Browser tracks, this text is displayed above the region
- score (optional) – an integer score for the region
- for BED files to be loaded as UCSC Genome Browser tracks, this should be a number between 0 and 1000, higher = "better"
- for non-GenBrowse BED files, this can be any integer value (e.g. the length of the region)
- strand (optional) - a single character describing the region's strand
- + – plus strand (Watson strand) region
- - – minus strand (Crick strand) region
- . – no strand – the region is not associated with a strand (e.g. a transcription factor binding region)
...
- The number of fields per line must be consistent throughout any single BED file
- e.g. they must all have 3 fields or all have 6 fields
- The first base on a contig is numbered 0
- versus 1 for BAM file positions
- so the a BED start of 99 is actually the 100th base on the contig
- but end positions are 1-based
- so a BED end of 200 is the 200th base on the contig
- the length of a BED region is end - start
- not end - start + 1, as it would be if both coordinates with 0-based or both 1-based
- this difference is one of the single greatest source of errors dealing with BED files!
...
- A BED3+ file contains the 3 required BED fields, followed by some number of user-defined columns (
- all records
- having the same number
- number of columns
- A BED6+ file contains the 3 required BED fields, 3 additional standard BED fields (name, score, strand), followed by some number of user-defined columns
- all records having the same number number columns
As we will see, BEDTools functions require BED3+ input files, or BED6+ if strand-specific operations are requested.
...
The BEDTools suite is a set of utilities for manipulating BED and BAM files. We call it the "Swiss army knife" for genomic region analyses because its sub-commands are so numerous and versatile. Some of the most common bedtools operations perform set-theory functions on regions: intersection (intersect), union (merge), set difference (subtract) – but there are many others. The table below lists some of the most useful sub-commands along with applicable use cases.
| Sub-command | Description | Use case(s) |
|---|---|---|
| bamtobed | Convert BAM files to BED format. | You want to have the contig, start, end, and strand information for each mapped alignment record in separate fields. Recall that the strand is encoded in a BAM flag (0x10) and the exact end coordinate requires parsing the CIGAR string. |
| bamtofastq | Extract FASTQ sequences from BAM alignment records. | You have downloaded a BAM file from a public database, but it was not aligned against the reference version you want to use (e.g. it is hg19 and you want an hg38 alignment). To re-process, you need to start with the original FASTQ sequences. |
| getfasta | Get FASTA entries corresponding to regions. | You want to run motif analysis, which requires |
| FASTA sequences, on a set of regions of interest. |
| In addition to |
| a BED or BAM file, you must provide FASTA file(s) for the genome/reference used for alignment (e.g. the FASTA file used to build the aligner index). | ||
| genomecov | Generate per-base genome-wide signal trace | Produce a per-base genome-wide signal (in bedGraph format), for example for a ChIP-seq or ATAC-seq experiment. After conversion to binary bigWig format, such tracks can be visualized in the Broad's IGV (Integrative Genome Browser) application, or configured in the UCSC Genome Browser as custom tracks. |
| coverage | Compute coverage of your regions |
|
|
In either case, regions (e.g. chromosomes or transcripts) are provided as a BED or GFF/GTF file. | ||
| multicov | Count overlaps between one or more BAM files and a set of regions of interest. | Count RNA-seq alignments that overlap a set of genes of interest. While this task is usually done with a specialized RNA-seq quantification tool (e.g. featureCounts or HTSeq), bedtools multicov can provide a quick estimate, e.g. for QC purposes. |
| intersect | Determine the overlap between two sets of regions. | Similar to multicov, but can also report the overlapping regions, not just count them. |
| merge | Combine a set of possibly-overlapping regions into a single set of non-overlapping regions. | Collapse overlapping gene annotations into per-strand non-overlapping regions. For example, to create non-overlapping transcipt regions before counting RNA-seq reads (e.g with featureCounts or HTSeq). If this is not done, the source regions will potentially be counted multiple times, once for each (overlapping) target region it intersects. |
| subtract | Remove unwanted regions. |
|
|
| ||
| closest | Find the genomic features nearest to a set of regions. | For a set of significant ChIP-seq |
| Transcription Factor (TF) binding regions ("peaks") that have been identified, determine nearby genes that may be targets of TF regulation. |
We will explore a few of these functions in our exercises.
...
Login to ls6, start and idev session, then load the BioContainers bedtools module, and then check its version.
| Code Block | ||||
|---|---|---|---|---|
| ||||
idev -m 120 -N 1 -A OTH21164 -r CoreNGS-Fricore-ngs-class-0606 # or, without a reservation idev -m 90120 -N 1 -A OTH21164 -p development module load biocontainers module load bedtools bedtools --version # should be bedtools v2.27.1 |
Input format considerations
- Most BEDTools functions now accept either BAM or BED files as input.
- BED format files must be BED3+, or BED6+ if strand-specific operations are requested.
- When comparing against a set of regions, those regions are usually supplied in either BED or GTF/GFF.
- All text-format input files (BED, GTF/GFFGFF3, VCF) should use Unix line endings (linefeed only).
The most important thing to remember about comparing regions using BEDTools, is that all input files must share the same set of contig names and be based on the same reference! For example, if an alignment was performed against a human GRCh38 reference genome from Gencode, use annotations from the corresponding GFFGTF/GTFGFF3 annotations.
About strandedness
By default many bedtools utilities that perform overlapping, consider reads overlapping the feature on either strand, but can be made strand-specific with the -s or -S option. This strandedness options for bedtools utilities refers the orientation of the R1 read with respect to the feature's (gene's) strand.
- -s says the R1 read is sense stranded (on the same strand as the gene).
- -S says the R1 read is antisense stranded (the opposite strand as the gene).
...
Which type of RNA-seq library you have depends on the library preparation method – so ask your sequencing center! Our yeast RNA-seq library is sense stranded (; however note that most RNA-seq libraries these days, including ones prepared by GSAF, are more often antisense stranded).
If you have a stranded RNA-seq library, you should use either -s or -S to avoid false counting against a gene on the wrong strand.
...
Annotation files that you retrieve from public databases are often in GTF (Gene Transfer Format) or one of the in GFF (General Feature Format) formats (usually GFF3 these days).
Unfortunately, both formats are obscure and hard to work with directly. While bedtools does accept annotation files in GTF/GFF/GTF format, you will not like the results. This is because the most useful information in a GTF/GFF/GTF file is in a loosely-structured attributes field.
Also unfortunately, there are a number of variations of both annotation formats However ; however both GTF and GFF share the first 8 Tab-separated fields:
- seqname - The name of the chromosome or scaffoldcontig.
- source - Name of the program that generated this feature, or other data source (e.g. public database)
- feature_type - Type of the feature. Examples of common feature types include:Some examples of common feature types are:, for example:
- chromosome
- CDS (coding sequence), exon
- gene, transcript
- start_codon, stop_codon
- start - Start position of the feature, with sequence numbering starting at 1.
- end - End position of the feature, with sequence numbering starting at 1.
- score - A numeric value. Often but not always an integer. strand - Meaning differs and not usually important.
- strand - Defined as + (forward), - (reverse), or . (no relevant strand)
- frame - For a CDS, one of 0, 1 or 2, specifying the reading frame of the first base; otherwise '.'
The Tab-separated columns will care about are (1) seqname, (3) feature_type and , (4,5) start, end and (7) strand. The reason we care is that when working with annotations, we usually only want to look at annotations of a particular type, most commonly gene, but also transcript or exon.
...
| Expand | |||||||
|---|---|---|---|---|---|---|---|
| |||||||
|
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
mkdir -p $SCRATCH/core_ngs/bedtools cd $SCRATCH/core_ngs/bedtools cp $CORENGS/yeast_rnaseq/sacCer3.R64-1-1_20110208.gff .yeast_mrna.sort.filt.bam* . cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* ./catchup/references/gff/sacCer3.R64-1-1_20110208.gff . # Use the less pager to look at multiple lines less sacCer3.R64-1-1_20110208.gff # Look at just the most-important Tab-separated columns cat sacCer3.R64-1-1_20110208.gff | grep -v '#' | cut -f 1,3-5,7 | head -20 # Include the ugly 9th column where attributes are stored cat sacCer3.R64-1-1_20110208.gff | grep -v '#' | cut -f 1,3,9 | head |
...
One of the first things you want to know about your annotation file is what gene features it contains. Here's how to find that: (Read more about what's going on here at Piping a histogram.)
| Expand | |||||
|---|---|---|---|---|---|
| |||||
| |||||
| Code Block | |||||
|
Read more about what's going on here at piping a histogram.
| Code Block | ||||
|---|---|---|---|---|
| ||||
cd $SCRATCH/core_ngs/bedtools cat sacCer3.R64-1-1_20110208.gff | grep -v '^#' | cut -f 3 | \ sort | uniq -c | sort -k1,1nr | more |
...
| Code Block | ||
|---|---|---|
| ||
cat sacCer3.R64-1-1_20110208.gff | grep -v '#^#' | \ awk 'BEGIN{FS=OFS="\t"}{ if($3=="gene"){print} }' \ > sc_genes.gff wc -l sc_genes.gff |
...
What most folks to is find some way to convert their GTF/GFF/GTF file to a BED file, parsing out some (or all) of the name/value attribute pairs into BED file columns after the standard 6. You can find such conversion programs on the web – or write one yourself. Or you could use the BioITeam conversion script, /work/projects/BioITeam/common/script/gtf_to_bed.pl. While it will not work 100% of the time, it manages to do a decent job on most GFFGTF/GTF filesGFF3 files. And it's pretty easy to run.
...
| Expand | |||||
|---|---|---|---|---|---|
| |||||
|
| Code Block | ||||
|---|---|---|---|---|
| ||||
/work/projects/BioITeam/common/script/gtf_to_bed.pl sc_genes.gff 1 \ > sc_genes.converted.bed |
...
The final transformation is to do a bit of re-ordering, dropping some fields. We'll do this with awk, because cut can't re-order fields. While this is not strictly required, it can be helpful to have the critical fields (including the gene ID) in the 1st 6 columns. We do this separately for the header line and the rest of the file so that the BED file we give bedtools does not have a header (– but we know what those fields are). We would normally preserve valuable annotation information such as Ontology_term, dbxref and Note, but drop them here for simplicity.
| Expand | |||||
|---|---|---|---|---|---|
| |||||
|
| Code Block | ||||
|---|---|---|---|---|
| ||||
head -1 sc_genes.converted.bed | sed 's/\r//' | awk '
BEGIN{FS=OFS="\t"}{print $1,$2,$3,$10,$5,$6,$15,$16}
' > sc_genes.bed.hdr
tail -n +2 sc_genes.converted.bed | sed 's/\r//' | awk '
BEGIN{FS=OFS="\t"}
{ # make sure gene name is populated
if($15 == "") {$15 = $10}
# make sure gene name is populated print $1,$2,$3,$10,$5,$6,$15,$16}
' > sc_genes.bed |
One final detail. Annotation files you download may have non-Unix (linefeed-only) line endings. Specifically, they may use Windows line endings (carriage return + linefeed). (Read about Line ending nightmares.) The expression sed 's/\r//' uses the sed (substitution editor) tool to replace carriage return characters ( \r ) with nothing, removing them from the output if they are present.
Finally, the 8 re-ordered attributes are:
...
| Code Block |
|---|
chrI 334 649 YAL069W 315 + YAL069W Dubious chrI 537 792 YAL068W-A 255 + YAL068WYAL068W-A Dubious chrI 1806 2169 YAL068C 363 - PAU8 Verified chrI 2479 2707 YAL067W-A 228 + YAL067W-A Uncharacterized chrI 7234 9016 YAL067C 1782 - SEO1 Verified chrI 10090 10399 YAL066W 309 + YAL066W Dubious chrI 11564 11951 YAL065C 387 - YAL065C Uncharacterized chrI 12045 12426 YAL064W-B 381 + YAL064W-B Uncharacterized chrI 13362 13743 YAL064C-A 381 - YAL064C-A Uncharacterized chrI 21565 21850 YAL064W 285 + YAL064W Verified chrI 22394 22685 YAL063C-A 291 - YAL063C-A Uncharacterized chrI 23999 27968 YAL063C 3969 - FLO9 Verified chrI 31566 32940 YAL062W 1374 + GDH3 Verified chrI 33447 34701 YAL061W 1254 + BDH2 Uncharacterized chrI 35154 36303 YAL060W 1149 + BDH1 Verified chrI 36495 36918 YAL059C-A 423 - YAL059C-A Dubious chrI 36508 37147 YAL059W 639 + ECM1 Verified chrI 37463 38972 YAL058W 1509 + CNE1 Verified chrI 38695 39046 YAL056C-A 351 - YAL056C-A Dubious chrI 39258 41901 YAL056W 2643 + GPB2 GPB2 Verified |
Note that value in the 8th column. In the yeast annotations from SGD there are 3 gene classifications: Verified, Uncharacterized and Dubious. The Dubious ones have no experimental evidence so are generally excluded.
| Expand | |||||
|---|---|---|---|---|---|
| |||||
|
Exercise: How many genes in our sc_genes.bed file are in each category?
| Expand | ||
|---|---|---|
| ||
Use cut to isolate that field (8), sort to sort the resulting values into blocks, then uniq -c to count the members of each block. |
...
Using the -c (column) and -o (operation) options, you can have add information added in subsequent fields. Each comma-separated column number following -c specifies a column to operate on, and the corresponding comma-separated function name following the -o specifies the operation to perform on that column in order to produce an additional output field.
...
We can do this with -c 6,4,4 -o distinct,count,collapse, which says that three custom output columns should be added:
- the 1st custom column (output column 4) should result from collapsing distinct (unique) values of gene file column 6 (the strand column, + or -)
- since we will ask for stranded merging, the merged regions will always be on the same strand, so this value will always be + or -
- the 2nd custom output column should result from counting the gene names in column 4 for all genes that were merged, and
- the 3rd custom output should be a comma-separated collapsed list of those same column 4 gene names.
bedtools merge also requires that the input BED file be sorted by locus (chrom + start), so we do that first, then we request a strand-specific merge (-s):
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
cd $SCRATCH/core_ngs/bedtools
sort -k1,1 -k2,2n sc_genes.bed > sc_genes.sorted.bed
bedtools merge -i sc_genes.sorted.bed -s \
-c 6,4,4 -o distinct,count,collapse > merged.sc_genes.txt |
The first few lines of the merged.sc_genes.txt file look like this (I've tidied it up a bit):
| Code Block |
|---|
2-micron 251 1523 + 1 R0010W 2-micron 1886 3008 - 1 R0020C 2-micron 3270 3816 + 1 R0030W 2-micron 5307 6198 - 1 R0040C chrI 334 792 + 2 YAL069W,YAL068W-A chrI 1806 2169 - 1 YAL068C chrI 2479 2707 + 1 YAL067W-A chrI 7234 9016 - 1 YAL067C chrI 10090 10399 10090 10399 + 1 YAL066W chrI 11564 11951 - 1 YAL065C |
As we specified:
- Output column 4 has the region's strand.
- Column 5 is the count of merged regions
...
- Column 6 is a collapsed, comma-separated list of the merged gene names
...
Exercise: Compare the number of regions in the merged and before-merge gene files.
...
Exercise: How many regions represent only 1 gene, 2 genes, or more?
| Expand | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| ||||||||||
| Expand | ||||||||||
| ||||||||||
Output column 5 has the gene count.
Produces this histogram:
There are 111 regions (105 + 4 + 1 + 1) where more than one gene contributed. |
Exercise: Repeat the steps above, but first create a good.sc_genes.bed file that does not include Dubious ORFs.
Or being fancy:
|
Exercise: Repeat the steps above, but first create a good.sc_genes.bed
...
file that does not include Dubious ORFs.
| Expand | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| |||||||||||||||
There were 5797 "good" (non-Dubious) genes before merging and 5770 after.
Produces this histogram:
Now there are only 20 regions where more than one gene was collapsed. Clearly eliminating the Dubious ORFs helped. |
...
To make a valid BED6 file, we'll include the first 3 output columns of merged.good.sc_genes.txt (chrom, start, end), but if strand is to be included, it should be in column 6. Column 4 should be name (we'll put the collapsed gene name list there), and column 5 a score (we'll put the region count there).
...
Make sure you're in an idev session, since we will be doing some significant computation, and make bedtools and samtools available.
| Expand | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| |||||||||
|
...
|
Copy over the yeast RNA-seq files we'll need (also copy the GFF gene annotation file if you didn't make one).
| Code Block | ||||
|---|---|---|---|---|
| ||||
# Get the merged yeast genes bed file if you didn't create one mkdir -p $SCRATCH/core_ngs/bedtools_multicov cd $SCRATCH/core_ngs/bedtools_multicov cp $CORENGS/catchup/bedtools_merge/merged*bed . # Copy the BAM file cd $SCRATCH/core_ngs/bedtools_multicov cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* . |
Exercises:
- How many reads are represented in the yeast_mrna.sort.filt.bam file?
- How many mapped? How many proper pairs? How many duplicates?
- What is the distribution of mapping qualities? What is the average mapping quality?
| Expand | ||
|---|---|---|
| ||
samtools flagstat for the different read counts. samtools view + cut + sort + uniq -c for mapping quality distribution samtools view + awk for average mapping quality |
...
| Expand | |||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| |||||||||||||||||||||||||
There are 3323242 total reads, all mapped and all properly paired. So this must be a quality-filtered BAM. There are 922114 duplicates, or about 28%. To get the distribution of mapping qualities :(BAM field 5)
To compute average mapping quality:
Mapping qualities range from 20 to 60 – excellent quality! Because the majority reads have mapping quality 60, the average is 59.2. So again, there must have been quality filtering performed on upstream alignment records. |
Here's how to run bedtools multicov in stranded mode, directing the standard output to a file:
| 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 |
ExerciseExercises:
- How may records of output were written?
- Where is the count of overlaps per output record?
| Expand | |||||
|---|---|---|---|---|---|
| |||||
6485 records were written, one for each feature in the merged.sc_genes.bed file. The overlap count was added as the last field in each output record (. So here it is field 7 , since the input annotation file had 6 columns). |
Exercise: How many features have non-zero overlap counts?
| Expand | |||||
|---|---|---|---|---|---|
| |||||
Most of the genes (6141/6485) have non-zero read overlap counts. |
...
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/)./download) or in the UCSC Genome Browser (https://genome.ucsc.edu/).
- Go to the UCSC Genome Browser: https://genome.ucsc.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 100 Vert. Cons track is a signal track
- the x-axis is the genome position
- the y-axis represents the base-wise conservation among vertebrates
- customize the 100 Vert. Cons track
- right-click on "100 Vert. Cons" text in the left margin,
- select "Configure 100 Vert. Cons" from the menu
- in the 100 Vert. Cons Track Settings dialog:
- change "Track height" to 100
- change "Data view scaling" to "auto-scale to data view"
- click "OK"
- right-click on "100 Vert. Cons" text in the left margin,
- 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 100 Vert. Cons track is a signal track
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.
Make sure you're in an idev session, then prepare a directory for this exercise.
| Expand | ||||||
|---|---|---|---|---|---|---|
| ||||||
|
...
| Code Block | |||||
|---|---|---|---|---|---|
| |||||
|
| Code Block | ||||
|---|---|---|---|---|
| ||||
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 bedtools genomecov -bg -ibam yeast_mrna.sort.filt.bam > yeast_mrna.genomecov.bedGraph wc -l yeast_mrna.genomecov.bedGraph # 1519274 lines |
...
Because this bedGraph file is for the small-ish (12Mb) yeast genome, and for reads that cover only part of that genome, it is not too big – only ~34M. But depending on the species and read depth, bedGraph files can get very large, so there is a coresponding corresponding binary format called bigWig (see https://genome.ucsc.edu/goldenpath/help/bigWig.html). The program to covert a bedGraph file to bigWig format is part of the UCSC Tools suite of programs. Look for it with module spider, and note that you can get information about all the tools in it using module spider with a specific container version:
...
Looking at the help for bedGraphToBigWig, we'll need a file of chromosome sizes. We can create one from our BAM header, using a Perl substitution script, which I prefer to sed(see Tips and tricks#perlpatternsubstitution):
| Code Block | ||
|---|---|---|
| ||
module load ucsc_tools cd $SCRATCH/core_ngs/bedtools_genomecov bedGraphToBigWig # look at its usage # create the needed chromosome sizes file from our BAM header module load samtools samtools view -H yeast_mrna.sort.filt.bam | grep -P 'SN[:]' | \ perl -pe 's/.*SN[:]//' | perl -pe 's/LN[:]//' > sc_chrom_sizes.txt cat sc_chrom_sizes.txt # displays: chrI 230218 chrII 813184 chrIII 316620 chrIV 1531933 chrV 576874 chrVI 270161 chrVII 1090940 chrVIII 562643 chrIX 439888 chrX 745751 chrXI 666816 chrXII 1078177 chrXIII 924431 chrXIV 784333 chrXV 1091291 chrXVI 948066 chrM 85779 |
...
| Code Block | ||
|---|---|---|
| ||
cd $SCRATCH/core_ngs/bedtools_genomecov export LC_COLLATE=C; # may or may not need this... 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.bedGraph > yeast_mrna.genomecov.sorted.bedGraph bedGraphToBigWig.bw |
See 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.genomecov.bedGraph sc_chrom_sizes.txt yeast_mrna.genomecov.bw |
...
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. Unfortunately, the ucsc-bigwigtobedgraph BioContainer seems to be brokenThat's a different BioContainers module (ucsc-bigwigtobedgraph) and the default container version doesn't work, so we'll use a version in the BioITeam area insteadspecifically 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 |
...