Content Comparison

...

# Request a 180 minute interactiveidev node on the normal queue using w/our reservation
idev -m 120 -N 1 -A OTH21164 -r CoreNGS
idev -m 120 -N 1 -A TRA23004 -r CoreNGScore-ngs-class-0606

# Request a 120 minute idevinteractive node on the development queue 
idev -m 120 -N 1 -A OTH21164 -p development
idev -m 120 -N 1 -A TRA23004 -p development

The BED format

...

1. chrom (required) – string naming the chromosome or other contig
2. start (required) – the 0-based start position of the region
3. end (required) – the 1-based end position of the region
4. name (optional) – an arbitrary string describing the region
   • for BED files loaded as UCSC Genome Browser tracks, this text is displayed above the region
5. 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)
6. 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
  with
  • 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.

...

...

idev -m 120 -N 1 -A OTH21164 -r CoreNGScore-ngs-class-0606

    # or, -Awithout TRA23004a # orreservation
idev -m 90120 -N 1 -A OTH21164 -p development 

# or -A TRA23004

modulemodule load biocontainers
module load bedtools
bedtools --version   # should be bedtools v2.27.1

...

• Most BEDTools functions 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

...

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:

...

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.

...

idev -m 120 -N 1 -A OTH21164 -r CoreNGS   core-ngs-class-0606

# or, -Awithout TRA23004
# ora reservation:
idev -m 90120 -N 1 -A OTH21164 -p development  # or -A TRA23004

module load biocontainers
module load bedtools
bedtools --version   # should be bedtools v2.27.1

...

mkdir -p $SCRATCH/core_ngs/bedtools
cd $SCRATCH/core_ngs/bedtools 
cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* .  
cp $CORENGS/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

...

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.

...

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.

...

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 $1print $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:

...

...

• 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.

...

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

...

cut -f 5 merged.sc_genes.txt | sort | uniq -c | sort -k2,2n

   6374 1
    105 2
      4 3
      1 4
      1 7

cut -f 5 merged.sc_genes.txt | sort | uniq -c | sort -k2,2n | \
  | grep -v ' 1$' | awk 'BEGIN{ct=0}{ct=ct+$1}END{print ct}'

...

Make sure you're in an idev session, since we will be doing some significant computation, and make bedtools and samtools available.

idev -m 120 -N 1 -A OTH21164 -r 

core-ngs-class-0606

# or, 

without 

a 

reservation

idev -m 

120 -N 1 -A OTH21164 -p development 

module load biocontainers
module load bedtools
bedtools --version   # should 

be 

bedtools v2.27.1

Copy over the yeast RNA-seq files we'll need (also copy the GFF gene annotation file if you didn't make one).

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

...

module load samtools

cd $SCRATCH/core_ngs/bedtools_multicov
samtools flagstat yeast_mrna.sort.filt.bam | tee yeast_mrna.flagstat.txt

3347559 + 0 in total (QC-passed reads + QC-failed reads)
24317 + 0 secondary
0 + 0 supplementary
922114 + 0 duplicates
3347559 + 0 mapped (100.00% : N/A)
3323242 + 0 paired in sequencing
1661699 + 0 read1
1661543 + 0 read2
3323242 + 0 properly paired (100.00% : N/A)
3323242 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools view yeast_mrna.sort.filt.bam | cut -f 5 | sort | uniq -c 

    498 20
   6504 21
   1012 22
    355 23
   1054 24
   2800 25
    495 26
  14133 27
    282 28
    358 29
    954 30
   1244 31
    358 32
   6143 33
    256 34
    265 35
   1112 36
    905 37
    309 38
   4845 39
   5706 40
    427 41
   1946 42
   1552 43
   1771 44
   6140 45
   1771 46
   3049 47
   3881 48
   3264 49
   4475 50
  15692 51
  25378 52
  16659 53
  18305 54
   7108 55
   2705 56
  59867 57
   2884 58
   2392 59
3118705 60

samtools view yeast_mrna.sort.filt.bam | awk '
  BEGIN{FS="\t"; sum=0; tot=0}
  {sum = sum + $5; tot = tot + 1}
  END{printf("mapping quality average: %.1f for %d reads\n", sum/tot,tot) }'

Here's how to run bedtools multicov in stranded mode, directing the standard output to a file:

idev -m 120 -N 1 -A OTH21164 -r CoreNGScore-ngs-class-0606

# or, without a reservation:
idev -m 120 -N 1 -A TRA23004 OTH21164 -p development

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 .
cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* .

...

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?

wc -l yeast_mrna_gene_counts.bed

...

cut -f 7 yeast_mrna_gene_counts.bed | grep -v '^0' | wc -l
# or
cut -f 7 yeast_mrna_gene_counts.bed | grep -v -c '^0'
# or
cat yeast_mrna_gene_counts.bed | \
  awk '{if ($7 > 0) print $7}' | wc -l

...

• 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 Layered H3K27Ac
  • the 100 Vert. Cons 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 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"
  • 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.

Make sure you're in an idev session, then prepare a directory for this exercise.

idev -m 120 -N 1 -A OTH21164 -r 

core-ngs-class-0606

# or, 

without 

a 

reservation:

idev -m 

120 -N 1 -A OTH21164 -p development 

module load biocontainers
module load bedtools
bedtools --version   # should be bedtools v2.27.1

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

...