2
$\begingroup$

And thanks in advance. I am debugging a text file which is the one the two inputs on the program RGmatch https://bitbucket.org/pfurio/rgmatch/src/master/.

The code works. Is that some records give division by zero.

I am doing it by hand.

But is there a way by making a bash script, to see which records show the problem?

I do not want to touch the python code.

python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b ../../test.bed -o myassoSort.txt

I have tried with classical debugging, but my actual question is how to automatize that?

Also there are no empty records. The bed file was created based on refseq.txt from the latest human assembly. And also the GTF file. On the refseq.txt I processed to obtain a bed file.

I have checked if this works on a single record, it does, but not in some specific ones. The error I get is:

  File "rgmatch.py", line 1311, in <module>
    main()
  File "rgmatch.py", line 287, in main
    run(gtf, dhs, outputfile)
  File "rgmatch.py", line 1153, in run
    pctg_area      = (float(region_overlap)/intron_length)*100
ZeroDivisionError: float division by zero

The record I found so far is this one.

chr19 14087850 14090751

I did grep on both GTF files and Refseq Files and in both files there is information about the record. Because the program search for areas between region of interest, If I put inside just the record only in GTF file and bed Files no error occurs.

grep "14090751" refseq.bed

692 NM_138352.3 chr19 - 14087850 14090751 14088416 14090420 6 14087850,14088670,14088879,14089067,14089784,14090084, 14088575,14088782,14088976,14089362,14090084,14090751, 0 SAMD1 cmpl cmpl 0,2,1,0,0,0,

grep "14090751" hg38.ncbiRefSeq.gtf

chr19   ncbiRefSeq.2021-02-10   exon    14090085    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";
chr19   ncbiRefSeq.2021-02-10   5UTR    14090421    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";
$\endgroup$
2
  • $\begingroup$ Please include the things I asked for over at Unix & Linux. We need to see the actual error message. Also, I call it "classical debugging" but I don't think anyone else will understand the term :) So, please edit your question, add the specific error message and then show us the commands that you want to automate. Have you checked it it works with just a single record? Does it work with just the first line of the bed file? $\endgroup$
    – terdon
    Jun 3 at 9:28
  • $\begingroup$ Also, if you still need help, feel free to ping me in the site chat room where we can directly debug the issue. $\endgroup$
    – terdon
    Jun 3 at 9:50
1
$\begingroup$

If you want to just start incrementally adding lines from the bed file until you hit the one that reproduces the problem, you can do something like this:

for i in $(seq 20 30 $(wc -l < ../../test.bed)); do
    head -n $i ../../test.bed > slice.bed
    if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.out
    then
        true
    else
        echo "Failed when running with the $i first lines."
        break
    fi
done
    

Explanation

seq 20 30 $(wc -l < ../../test.bed) will iterate over all numbers between 20 and the number of lines in your bed file with increments of 30. Next, the head -n $i ../../test.bed > slice.bed will take the first $i lines of the bed file. This way, the loop will run the rgmatch command on the first 20 lines of the file, then the first 50, then the first 80 and so on as long as the command runs successfully. When the rgmatch command fails, it will report the value of $i to help you find the line that caused the problem.

$\endgroup$
5
  • $\begingroup$ Thank you @terdon I will edit the answer and try this.. If not I will appreciate very much your help. Right now I have a work meeting. But after that, I will do everything you said to me. Thank you ! $\endgroup$ Jun 3 at 9:55
  • $\begingroup$ De res, @EdmondGeraudAguilar! And do feel free to ping me in chat if you need more help. $\endgroup$
    – terdon
    Jun 3 at 9:59
  • $\begingroup$ It works ! I mean, I would that the script makes like a try-catch then store the results.. but I will handle with this. Than you very much $\endgroup$ Jun 3 at 11:23
  • $\begingroup$ @EdmondGeraudAguilar it's very easy to extend to whatever you need. But this is why I asked you to edit and explain what commands you want to automate. If you want to keep the results of each run, just change the rgmatch output file: if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.$i.out; then ... $\endgroup$
    – terdon
    Jun 3 at 11:27
  • $\begingroup$ For those who may use this function be aware that some start coordinates of transcripts are missing in the GTF files that is the problem $\endgroup$ Jun 3 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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