ylab-hi / ScanExitron

A computational workflow for exitron splicing identification
MIT License
12 stars 6 forks source link

ScanExitron test.bam returned error #4

Closed zyh4482 closed 3 years ago

zyh4482 commented 3 years ago

I've installed all required packages for ScanExitron.py. When I give a try on the test file, it returns with error message as follow:

`tomas@tomas-virtual-machine:~/ScanExitron/tests$ python ~/ScanExitron/ScanExitron.py -i test.bam -r hg38 Checking for 'regtools': found /home/tomas/regtools/build/regtools

Checking for 'bedtools': found /home/tomas/anaconda3/bin/bedtools

Checking for 'samtools': found /home/tomas/anaconda3/bin/samtools

error! Command 'samtools view -q 0 -@ 1 -O BAM -o test.hq.bam test.bam && samtools index test.hq.bam' returned non-zero exit status 127. ` May I ask what's going on there? Did I make any mistake? Thank you.

dolittle007 commented 3 years ago

Can you check your SamTools version? If SamTools version < 1.0, the "-O" parameter will be not available.

zyh4482 commented 3 years ago

I got this problem solved. "Conda install" via bioconda channel will automatically install samtool v1.7. Even I forced conda to download version "1.10", it still reports many conflicts. Somehow samtools v1.7 seems not to be compatible and reported this error. I add another channel conda-forge, at this time it could be successfully installed.

tomas@tomas-virtual-machine:~/ScanExitron/aaa$ samtools --version samtools 1.10 Using htslib 1.10.2 Copyright (C) 2019 Genome Research Ltd.

BTW, to avoid such mistakes, I checked all the other packages. The versions of the other packages should work. (bedtools 2.30, pyfaidx 0.5.9, bedops 2.4.39, Regtools 0.4.2)

But another error occurs as shown below. The command 'samtools bedcov test.hq.position.bed test.bam -Q 0' returned error message.

tomas@tomas-virtual-machine:~/ScanExitron/aaa$ python /home/tomas/ScanExitron/ScanExitron.py -i test.bam -r hg38
Checking for 'regtools': found /home/tomas/regtools/build/regtools

Checking for 'bedtools': found /home/tomas/anaconda3/bin/bedtools

Checking for 'samtools': found /home/tomas/anaconda3/bin/samtools

MAPQ filtering finished!
regtools junctions extract -i 5 -I 10000000 test.hq.bam -o test.hq.bed
Calling junctions start
Calling junctions finished!
regtools junctions annotate /home/tomas/ScanExitron/aaa/test.hq.bed /home/tomas/human_genome/hg38.fa /home/tomas/human_genome/gencode.v21.annotation.gtf -o test.hq.janno
test.hq.janno generated!
test.hq.janno generated!
Reading test.hq.janno
bedtools intersect -s -wo -a 28005.junction.bed -b /home/tomas/human_genome/gencode.hg38.CDS.bed > 28005.overlap.bed
Junctions intersect with CDS
Junctions intersect with CDS finished!
Reading BAM file: test.bam
samtools bedcov test.hq.position.bed test.bam -Q 0
Error happend!: Command 'samtools bedcov test.hq.position.bed test.bam -Q 0' returned non-zero exit status 2.
b"[E::idx_find_and_load] Could not retrieve index file for 'test.bam'\nERROR: fail to open index BAM file 'test.bam'\n"
Traceback (most recent call last):
  File "/home/tomas/ScanExitron/ScanExitron.py", line 412, in <module>
    main()
  File "/home/tomas/ScanExitron/ScanExitron.py", line 406, in main
    percent_spliced_out(bam_file=args.input, src_exitron_file=src_exitron_file, position_bed_file=position_bed_file, ao_cutoff=args.ao, pso_cutoff=args.pso, mapq=args.mapq)
  File "/home/tomas/ScanExitron/ScanExitron.py", line 304, in percent_spliced_out
    five_prime_reads = depth_dict['{}\t{}'.format(chrm, start)] - ao
KeyError: 'chr17\t50071324'
dolittle007 commented 3 years ago

I think SamTools v1.7 will work, as well. The error message showed that you didn't have the INDEX file for the "test.bam". You need test.bam.bai or test.bai.

zyh4482 commented 3 years ago

My bad. I moved bam to a new folder but forgot the index.

After generating test.bam.bai using samtools index test.bam, it can output a "test.hq.exitron" file but contains no exitrons.

I don't know if it is the problem with src_exitron_file so that the function main() did not remove janno file and left "test.hq.bam", "test.hq.bai", "test.hq.janno" etc. Would you mind telling me how to debug this? Sorry to interrupt. Thank you again.

Checking for 'regtools': found /home/tomas/regtools/build/regtools

Checking for 'bedtools': found /home/tomas/anaconda3/bin/bedtools

Checking for 'samtools': found /home/tomas/anaconda3/bin/samtools

MAPQ filtering finished!
regtools junctions extract -i 5 -I 10000000 test.hq.bam -o test.hq.bed
Calling junctions start
Calling junctions finished!
regtools junctions annotate /home/tomas/ScanExitron/aaa/test.hq.bed /home/tomas/human_genome/hg38.fa /home/tomas/human_genome/gencode.v21.annotation.gtf -o test.hq.janno
test.hq.janno generated!
test.hq.janno generated!
Reading test.hq.janno
bedtools intersect -s -wo -a 3866.junction.bed -b /home/tomas/human_genome/gencode.hg38.CDS.bed > 3866.overlap.bed
Junctions intersect with CDS
Junctions intersect with CDS finished!
Reading BAM file: test.bam
samtools bedcov test.hq.position.bed test.bam -Q 0
Calculate PSO and PSI.
Finished reading BAM file: test.bam
dolittle007 commented 3 years ago

It is not a bug. You can set "-p 0", these two exitron splicing events will be there.

zyh4482 commented 3 years ago

Thank you so much. It works.