How can I count read that falls into exon-intron and exon-exon junction?

Question

[biplab]

I have gtf file as follows

VI  ensembl gene    53260   54696   .   -   .   gene_id "YFL039C"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding";
VI  ensembl transcript  53260   54696   .   -   .   gene_id "YFL039C"; transcript_id "YFL039C"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding";
VI  ensembl exon    54687   54696   .   -   .   gene_id "YFL039C"; transcript_id "YFL039C"; exon_number "1"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "YFL039C.1";
VI  ensembl CDS 54687   54696   .   -   0   gene_id "YFL039C"; transcript_id "YFL039C"; exon_number "1"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "YFL039C";
VI  ensembl start_codon 54694   54696   .   -   0   gene_id "YFL039C"; transcript_id "YFL039C"; exon_number "1"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding";
VI  ensembl exon    53260   54377   .   -   .   gene_id "YFL039C"; transcript_id "YFL039C"; exon_number "2"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "YFL039C.2";
VI  ensembl CDS 53260   54377   .   -   1   gene_id "YFL039C"; transcript_id "YFL039C"; exon_number "2"; gene_name "ACT1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_name "ACT1"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "YFL039C";

How can I use this gtf file to make a bed file which I will be able to use for counting reads that falls into exon-exon and exon-intron junction?

[lieven.sterck]

Since it looks like your mapping to genomic sequences, how (and why?) would you count reads to fall into exon-exon? you will need some pretty long reads and/or very small introns to get into such a situation. Moreover any read that spans a exon-exon junction will be definition also span (2) exon-intron junctions.

I can understand the exon-intron part and then Alex Reynolds will work nicely

Answer 1

[Alex Reynolds]

Steps via awk and BEDOPS.

Build a BED file of exons:

$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gtf.gz \
    | gunzip -c - \
    | awk '($3=="exon")' - \
    | gtf2bed - \
    | cut -f1-6 - \
    > gencode.v28.annotation.exons.c1t6.bed

Convert transcripts to merged exons:

$ awk -f transcripts2mergedExons.awk gencode.v28.annotation.exons.c1t6.bed > gencode.v28.annotation.mergedExons.c1t6.bed

Make an exon-intron list:

$ awk -f mergedExons2exonIntronList.awk gencode.v28.annotation.mergedExons.c1t6.bed > gencode.v28.annotation.exonsAndIntrons.c1t6.bed

Convert these to junctions:

$ awk -f exonIntronList2JunctionList.awk gencode.v28.annotation.exonsAndIntrons.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.c1t6.bed

Pad them, e.g. by 25 nt around the junction (adjust as needed):

$ bedops --everything --range 25 gencode.v28.annotation.exonsIntronJunctions.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed

Map reads to the padded junctions:

$ bedmap --echo --count --delim '\t' gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed <(bam2bed < reads.bam) > answer.bed

The file answer.bed will contain the junction and the number of reads that map to — overlap with — the junction, by one or more bases.

Some links to Github Gists with listed awk scripts:

• 	BEGIN {
  	FS="\t";
  	old_chr = "chrN";
  	old_start = -1;
  	old_stop = -1;
  	old_id = "*";
  	old_score = "*";
  	old_strand = "*";
  	new_element_flag = 0;
  	}
  	{
  	new_chr = $1;
  	new_start = $2;
  	new_stop = $3;
  	new_id = $4;
  	new_score = $5;
  	new_strand = $6;
  	if (old_id != new_id) {
  	if (old_id != "*") {
  	old_score = 1; # exon
  	print old_chr"\t"old_start"\t"old_stop"\t"old_id"\t"old_score"\t"old_strand;
  	}
  	new_score = 1; # exon
  	print new_chr"\t"new_start"\t"new_stop"\t"new_id"\t"new_score"\t"new_strand;
  	old_id = new_id;
  	new_element_flag = 1;
  	}
  	else {
  	# find genomic difference between new and old elements
  	diff_start = new_start - old_stop;
  	if (diff_start > 0) {
  	# we have an intron between new and old elements
  	if (new_element_flag == 0) {
  	# print the current old element
  	old_score = 1; # we know it is an exon
  	print old_chr"\t"old_start"\t"old_stop"\t"old_id"\t"old_score"\t"old_strand;
  	}
  	else {
  	# we are no longer at the start of a new element
  	new_element_flag = 0;
  	}
  	}
  	else {
  	# alternate splice? some other transcript of the same exon?
  	# preserve the old start coordinate
  	if (old_start != -1) {
  	new_start = old_start;
  	}
  	if (new_stop < old_stop) {
  	new_stop = old_stop;
  	}
  	}
  	}
  	old_chr = new_chr;
  	old_start = new_start;
  	old_stop = new_stop;
  	old_id = new_id;
  	old_score = new_score;
  	old_strand = new_strand;
  	}
  	END {
  	}

  view raw transcripts2mergedExons.awk hosted with ❤ by GitHub
  (transcripts2mergedExons.awk)

• 	BEGIN {
  	FS="\t";
  	old_chr = "chrN";
  	old_start = -1;
  	old_stop = -1;
  	old_id = "*";
  	old_score = "*";
  	old_strand = "*";
  	}
  	{
  	new_chr = $1;
  	new_start = $2;
  	new_stop = $3;
  	new_id = $4;
  	new_score = $5;
  	new_strand = $6;
  	if (old_id != new_id) {
  	new_score = 1; # exon
  	print new_chr"\t"new_start"\t"new_stop"\t"new_id"\t"new_score"\t"new_strand;
  	old_id = new_id;
  	}
  	else {
  	# we construct and print an intron, then print the exon
  	old_score = 0;
  	print old_chr"\t"old_stop"\t"new_start"\t"old_id"\t"old_score"\t"old_strand;
  	new_score = 1; # exon
  	print new_chr"\t"new_start"\t"new_stop"\t"new_id"\t"new_score"\t"new_strand;
  	}
  	old_chr = new_chr;
  	old_start = new_start;
  	old_stop = new_stop;
  	old_id = new_id;
  	old_score = new_score;
  	old_strand = new_strand;
  	}
  	END {
  	}

  view raw mergedExons2exonIntronList.awk hosted with ❤ by GitHub
  (mergedExons2exonIntronList.awk)

• 	BEGIN {
  	FS="\t";
  	old_chr = "chrN";
  	old_start = -1;
  	old_stop = -1;
  	old_id = "*";
  	old_score = "*";
  	old_strand = "*";
  	}
  	{
  	new_chr = $1;
  	new_start = $2;
  	new_stop = $3;
  	new_id = $4;
  	new_score = $5;
  	new_strand = $6;
  	if (old_id != new_id) {
  	# new element
  	old_id = new_id;
  	}
  	else {
  	# construct and print junction
  	old_score = 2;
  	print old_chr"\t"old_stop"\t"(old_stop+1)"\t"old_id"\t"old_score"\t"old_strand;
  	}
  	old_chr = new_chr;
  	old_start = new_start;
  	old_stop = new_stop;
  	old_id = new_id;
  	old_score = new_score;
  	old_strand = new_strand;
  	}
  	END {
  	}

  view raw exonIntronList2JunctionList.awk hosted with ❤ by GitHub
  (exonIntronList2JunctionList.awk)

[biplab]

Thank you so much. I will try this.