search

dieterich-lab / rp-bp

Rp-Bp is a Bayesian approach to predict, at base-pair resolution, ribosome occupancy and translation.
MIT License
7 stars 5 forks source link

Error using de novo assembly #104

Closed HeyLifeHD closed 4 years ago

HeyLifeHD commented 4 years ago

Dear authors,

We are able to run rp-bp now with "default" options. However when only adding the de_novo_gtf in the config file, we are facing some problems during prepare_rpbp_genome. The error is due to the gtf-to-bed12 script:

gtf-to-bed12 /data/LTR/191204_Riboseq/191217_rbbp/gffCompare.annotated.sorted.gff /data/LTR/191204_Riboseq/191217_rbbp/index/hg19ensembl75_deNovo/hg19ensembl75.de-novo.bed.gz --chr-name-file /data/LTR/191204_Riboseq/191217_rbbp/index/hg19ensembl75star_deNovo/chrName.txt --num-cpus 10 --logging-level WARNING --file-logging-level NOTSET --stdout-logging-level NOTSET --stderr-logging-level NOTSET
100%|#################################################################################################################################################################################| 500/500 [00:59<00:00,  8.44it/s]
Traceback (most recent call last):
  File "/usr/local/bin/gtf-to-bed12", line 8, in <module>
    sys.exit(main())
  File "/usr/local/lib/python3.6/site-packages/pbio/utils/pgrms/gtf_to_bed12.py", line 206, in main
    bed12_df = bed12_df.merge(cds_start_df, on='id', how='left')
  File "/usr/local/lib/python3.6/site-packages/pandas/core/frame.py", line 6877, in merge
    copy=copy, indicator=indicator, validate=validate)
  File "/usr/local/lib/python3.6/site-packages/pandas/core/reshape/merge.py", line 47, in merge
    validate=validate)
  File "/usr/local/lib/python3.6/site-packages/pandas/core/reshape/merge.py", line 529, in __init__
    self.join_names) = self._get_merge_keys()
  File "/usr/local/lib/python3.6/site-packages/pandas/core/reshape/merge.py", line 846, in _get_merge_keys
    left_keys.append(left._get_label_or_level_values(lk))
  File "/usr/local/lib/python3.6/site-packages/pandas/core/generic.py", line 1697, in _get_label_or_level_values
    raise KeyError(key)

It seems as it is expecting the specification of CDS. However, as explained in your documentation, this is a de novo assembly based on RNAseq. CDS are supposed to be specified via rp-bp. As an example here are the first 6 lines of each gtf: Proper reference:

chr1    pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene";
chr1    processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana";
chr1    processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";
chr1    processed_transcript    exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_number "2"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00003582793";
chr1    processed_transcript    exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_number "3"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002312635";

De novo reference:

GL000192.1      HAVANA  transcript      493155  495565  .       -       .       transcript_id "ENST00000614199.1_1"; gene_id "ENSG00000277655.1_5"; gene_name "AC245407.1"; xloc "XLOC_000001"; ref_gene_id "ENSG00000277655.1_5"; cmp_ref "ENST00000614199.1_1"; class_code "="; tss_id "TSS1";
GL000192.1      HAVANA  exon    493155  493368  .       -       .       transcript_id "ENST00000614199.1_1"; gene_id "ENSG00000277655.1_5"; exon_number "1";
GL000192.1      HAVANA  exon    495329  495565  .       -       .       transcript_id "ENST00000614199.1_1"; gene_id "ENSG00000277655.1_5"; exon_number "2";
GL000193.1      HAVANA  transcript      49232   79916   .       -       .       transcript_id "ENST00000623554.3_1"; gene_id "ENSG00000280081.3_5"; gene_name "LINC01667"; xloc "XLOC_000003"; ref_gene_id "ENSG00000280081.3_5"; cmp_ref "ENST00000623554.3_1"; class_code "="; tss_id "TSS3";
GL000193.1      HAVANA  exon    49232   49618   .       -       .       transcript_id "ENST00000623554.3_1"; gene_id "ENSG00000280081.3_5"; exon_number "1";
GL000193.1      HAVANA  exon    50562   50619   .       -       .       transcript_id "ENST00000623554.3_1"; gene_id "ENSG00000280081.3_5"; exon_number "2";
GL000193.1      HAVANA  exon    59869   59956   .       -       .       transcript_id "ENST00000623554.3_1"; gene_id "ENSG00000280081.3_5"; exon_number "3";
GL000193.1      HAVANA  transcript      73408   81323   .       -       .       transcript_id "ENST00000623919.1_1"; gene_id "ENSG00000280081.3_5"; gene_name "LINC01667"; xloc "XLOC_000003"; ref_gene_id "ENSG00000280081.3_5"; cmp_ref "ENST00000623919.1_1"; class_code "="; tss_id "TSS4";
GL000193.1      HAVANA  exon    73408   74416   .       -       .       transcript_id "ENST00000623919.1_1"; gene_id "ENSG00000280081.3_5"; exon_number "1";
GL000193.1      HAVANA  transcript      74189   77037   .       -       .       transcript_id "ENST00000623933.1_1"; gene_id "ENSG00000280081.3_5"; gene_name "LINC01667"; xloc "XLOC_000003"; ref_gene_id "ENSG00000280081.3_5"; cmp_ref "ENST00000623933.1_1"; class_code "="; tss_id "TSS5";

Thanks for your support.

Best,

Joschka

HeyLifeHD commented 4 years ago

After reformatting with gffread gffread gffCompare.annotated.sorted.gff -T -o gffCompare.annotated.sorted_test.gff the gtf-to-bed12 script is working. I guess the additional information in column number 9 was bothering the program. Sorry for the inconvenience.

Best,

Joschka