4
$\begingroup$

I making a bed file for RSeQC, so it can do things like compute the number of reads from exons, introns, 5"UTRs, etc.

I want to use a bed file that corresponds to my GTF file, so I use gtf2bed to make a bed file, like this:

awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' Homo_sapiens.GRCh38.gtf | gtf2bed - > $Homo_sapiens.GRCh38.bed

This file throws some errors, among other when using it with junction_annotation.py:

Reading reference bed file:  /rst1/2017-0205_illuminaseq/data/references/Reference_Genomes/GRCh38.95/Homo_sapiens.GRCh38.87.bed  ... Traceback (most recent call last):
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 125, in <module>
    main()
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 109, in main
    obj.annotate_junction(outfile=options.output_prefix,refgene=options.ref_gene_model,min_intron=options.min_intron, q_cut = options.map_qual)
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/lib/python2.7/site-packages/qcmodule/SAM.py", line 3762, in annotate_junction
    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
ValueError: invalid literal for int() with base 10: 'gene_version'

and when I compare this bed file to the RSeQC provide (dated) one these are the differences:

My bed file I created:

1       11868   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       11868   14409   ENSG00000223972 .       +       havana  gene    .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_id "";
1       11868   14409   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
1       12009   12057   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
1       12009   13670   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
1       12178   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12697   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12721   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       12974   13052   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       13220   13374   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001746346"; exon_version "2"; tag "basic"; transcript_support_level "NA";

The RSeQC provided bed file:

1       11868   14409   ENST00000456328 0       +       11868   14409   0       3       359,109,1189,   0,744,1352,     DDX11L1
1       12009   13670   ENST00000450305 0       +       12009   13670   0       6       48,49,85,78,154,218,    0,169,603,965,1211,1443,        DDX11L1
1       17368   17436   ENST00000619216 0       -       17368   17436   0       1       68,     0,      MIR6859-1
1       14403   29570   ENST00000488147 0       -       14403   29570   0       11      98,34,152,159,198,136,137,147,99,154,37,        0,601,1392,2203,2454,2829,3202,3511,3864,10334,15130,   WASH7P
1       29553   31097   ENST00000473358 0       +       29553   31097   0       3       486,104,122,    0,1010,1422,    MIR1302-2
1       30266   31109   ENST00000469289 0       +       30266   31109   0       2       401,134,        0,709,  MIR1302-2
1       30365   30503   ENST00000607096 0       +       30365   30503   0       1       138,    0,      MIR1302-2
1       34553   36081   ENST00000417324 0       -       34553   36081   0       3       621,205,361,    0,723,1167,     FAM138A
1       35244   36073   ENST00000461467 0       -       35244   36073   0       2       237,353,        0,476,  FAM138A
1       52472   53312   ENST00000606857 0       +       52472   53312   0       1       840,    0,      OR4G4P

Why are the files so different? How would I make a simpler bed file form my GTF file?

Improve this question
$\endgroup$
1
  • $\begingroup$ Please show us an excerpt of your GTF file (a few lines) and i) the output you get from it and ii) the output you want to get from it. As it stands, I don't understand what you mean by "simpler bed file" nor why the files are different (can't know that if we don't know what data you are starting from). $\endgroup$
    – terdon ♦
    Commented Feb 21, 2019 at 16:40

1 Answer 1

Reset to default
3
$\begingroup$

UCSC Utilities

Such output results when using the UCSC utilities gtfToGenePred and genePredToBed in series. Builds for macOS (x86_64) and Linux are available directly from UCSC. Alternatively, Conda builds are available through the Bioconda channel, under the package names ucsc-gtftogenepred and ucsc-genepredtobed.

Example: Homo sapiens (Ensembl Release 106)

## (optional install through Conda)
#mamba install -c conda-forge -c bioconda ucsc-gtftogenepred ucsc-genepredtobed

## download GTF
wget ftp://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz

gzip -cd Homo_sapiens.GRCh38.106.gtf.gz |\
  gtfToGenePred /dev/stdin /dev/stdout |\
  genePredToBed /dev/stdin /dev/stdout |\
  head

Output BED

1   1211339 1214153 ENST00000379236 0   -   1211554 1214127 0   7   286,129,197,67,102,123,171, 0,364,602,1298,1652,2323,2643,
1   1211339 1214138 ENST00000497869 0   -   1214138 1214138 0   5   493,197,67,794,156, 0,602,1298,1652,2643,
1   1212018 1213498 ENST00000453580 0   -   1213498 1213498 0   4   120,67,102,104, 0,619,973,1376,
1   1203507 1206571 ENST00000328596 0   -   1203590 1206571 0   4   453,88,123,187, 0,891,1862,2877,
1   1203507 1206592 ENST00000379268 0   -   1203843 1206571 0   5   461,203,88,123,208, 0,526,891,1862,2877,
1   1203843 1205680 ENST00000486728 0   -   1203843 1205463 0   4   125,203,88,311, 0,190,555,1526,
1   1203843 1206571 ENST00000379265 0   -   1203843 1206571 0   5   125,182,88,123,187, 0,211,555,1526,2541,
1   1471764 1497848 ENST00000673477 0   +   1471884 1495817 0   16  325,77,102,60,70,166,70,156,57,126,125,52,71,168,109,2364,  0,5509,6879,7284,9102,10373,10780,13251,14017,14345,14779,16098,17439,18492,18798,23720,
1   1478025 1497848 ENST00000472194 0   +   1497848 1497848 0   14  720,60,70,166,70,156,57,126,125,52,71,168,109,2364, 0,1023,2841,4112,4519,6990,7756,8084,8518,9837,11178,12231,12537,17459,
1   1479048 1482662 ENST00000378736 0   +   1482662 1482662 0   4   60,70,166,118,  0,1818,3089,3496,
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanx for this answer, I'm not really in a position to check this quickly anymore, but I'll label this as solved. $\endgroup$
    – Freek
    Commented Apr 22, 2022 at 13:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.