secastel / phaser

phasing and Allele Specific Expression from RNA-seq
GNU General Public License v3.0
111 stars 36 forks source link

Error in using phaser_gene_ae.py #78

Open c-nabokov opened 2 years ago

c-nabokov commented 2 years ago

Dear developers,

I am very grateful to you for developing this software so that I can use it to explore the issues related to the expression of allelic genes in soybean hybrids.

However,this error "Traceback (most recent call last): File "/home/nabokov/software/phaser/phaser_gene_ae/phaser_gene_ae.py", line 240, in main(); File "/home/nabokov/software/phaser/phaser_gene_ae/phaser_gene_ae.py", line 111, in main mapped_reads = variant_feature_reads(row,feature); File "/home/nabokov/software/phaser/phaser_gene_ae/phaser_gene_ae.py", line 204, in variant_feature_reads hap_a_reads += str(row['aReads']).split(";")[xvar_index].split(","); IndexError: list index out of range" occurred when I calculated allelic gene expression levels (phaser_gene_ae.py). I can't understand whether this error is caused by my input file or software bug.

The following is a brief description of the code and files to run. I would like to ask you to help me solve this problem.

Thank you.

code: python /home/nabokov/software/phaser/phaser_gene_ae/phaser_gene_ae.py --haplotypic_counts test/F12_1_output.haplotypic_counts.txt --features genome/wm82a2.gene.bed --id_separator ';' --o F12_1.gene_ae.txt

bed.file: Chr01 27354 27655 . . - phytozomev10 three_prime_UTR . ID=Glyma.01G000100.1.Wm82.a2.v1.three_prime_UTR.1;Parent=Glyma.01G000100.1.Wm82.a2.v1;pacid=30544134 Chr01 27354 28320 . . - phytozomev10 gene . ID=Glyma.01G000100.Wm82.a2.v1;Name=Glyma.01G000100;ancestorIdentifier=Glyma01g00210.v1.1 Chr01 27354 28320 . . - phytozomev10 mRNA . ID=Glyma.01G000100.1.Wm82.a2.v1;Name=Glyma.01G000100.1;pacid=30544134;longest=1;Parent=Glyma.01G000100.Wm82.a2.v1 Chr01 27655 27824 . . - phytozomev10 CDS 1 ID=Glyma.01G000100.1.Wm82.a2.v1.CDS.3;Parent=Glyma.01G000100.1.Wm82.a2.v1;pacid=30544134 Chr01 27925 27991 . . - phytozomev10 CDS 1 ID=Glyma.01G000100.1.Wm82.a2.v1.CDS.2;Parent=Glyma.01G000100.1.Wm82.a2.v1;pacid=30544134 Chr01 28138 28218 . . - phytozomev10 CDS 0 ID=Glyma.01G000100.1.Wm82.a2.v1.CDS.1;Parent=Glyma.01G000100.1.Wm82.a2.v1;pacid=30544134 Chr01 28218 28320 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000100.1.Wm82.a2.v1.five_prime_UTR.1;Parent=Glyma.01G000100.1.Wm82.a2.v1;pacid=30544134 Chr01 58974 59111 . . - phytozomev10 three_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.three_prime_UTR.1;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 58974 67527 . . - phytozomev10 gene . ID=Glyma.01G000200.Wm82.a2.v1;Name=Glyma.01G000200 Chr01 58974 67527 . . - phytozomev10 mRNA . ID=Glyma.01G000200.1.Wm82.a2.v1;Name=Glyma.01G000200.1;pacid=30543475;longest=1;Parent=Glyma.01G000200.Wm82.a2.v1 Chr01 62005 62045 . . - phytozomev10 CDS 1 ID=Glyma.01G000200.1.Wm82.a2.v1.CDS.5;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 62566 62644 . . - phytozomev10 CDS 1 ID=Glyma.01G000200.1.Wm82.a2.v1.CDS.4;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 63065 63141 . . - phytozomev10 CDS 2 ID=Glyma.01G000200.1.Wm82.a2.v1.CDS.3;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 63333 63417 . . - phytozomev10 CDS 2 ID=Glyma.01G000200.1.Wm82.a2.v1.CDS.2;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 64054 64061 . . - phytozomev10 CDS 0 ID=Glyma.01G000200.1.Wm82.a2.v1.CDS.1;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 64061 64199 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.five_prime_UTR.5;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 65204 65261 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.five_prime_UTR.4;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 65453 65537 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.five_prime_UTR.3;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 66175 66320 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.five_prime_UTR.2;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 67327 67527 . . - phytozomev10 five_prime_UTR . ID=Glyma.01G000200.1.Wm82.a2.v1.five_prime_UTR.1;Parent=Glyma.01G000200.1.Wm82.a2.v1;pacid=30543475 Chr01 67769 67787 . . + phytozomev10 CDS 0 ID=Glyma.01G000300.1.Wm82.a2.v1.CDS.1;Parent=Glyma.01G000300.1.Wm82.a2.v1;pacid=30545121

secastel commented 2 years ago

Hmmm, not sure what is going on here. Could you provide both full input files (wm82a2.gene.bed, F12_1_output.haplotypic_counts.txt) so I can try to reproduce the error?

c-nabokov commented 2 years ago

Hmmm, not sure what is going on here. Could you provide both full input files (wm82a2.gene.bed, F12_1_output.haplotypic_counts.txt) so I can try to reproduce the error?

Of course. Thank you very much for your reply and help.

I have uploaded the file to the following link.

https://github.com/c-nabokov/Test-files

Thanks for your help again.

lintingyi2014 commented 7 months ago

I am also encountering this error

[lint6@cn4287 SCRIPTS]$ phaser_gene_ae.py --haplotypic_counts phaser_test.haplotypic_counts.txt --features gencode.v26.GRCh38.genes.bed.gz --o phaser_test_case_gene_ae.txt

################################################## Welcome to phASER Gene AE v1.2.0 Author: Stephane Castel (scastel@nygenome.org) ##################################################

1 Loading features...

Traceback (most recent call last): File "/usr/local/bin/phaser_gene_ae.py", line 241, in main(); File "/usr/local/bin/phaser_gene_ae.py", line 47, in main start = int(columns[1]); IndexError: list index out of range

lintingyi2014 commented 7 months ago

somehow it needed to be re-tabbed

It works with: zcat gencode.v26.GRCh38.genes.bed.gz | sed 's/\^I/\t/g' | awk 'BEGIN{OFS="\t"} {if(NF==4) print $0; else print "Invalid line:", $0}' > re_tabbed_file.bed

secastel commented 7 months ago

When I try this I am not seeing any lines come up as "Invalid". Are you? If so can you please share which? When you say it is now working, are you using the ".bed" file you generated with the above code as input, or are you gzipping it again?

netwon123 commented 4 months ago

Hellow, Could you solved this bug?

Soymethylation commented 2 months ago

Dear colleagues,

In my data, I solved this problem.

My mistake was that I didn't use the right bed file. After I changed the bed file to four columns, it worked fine.

As the software developer's tutorial said, the four columns of the bed file are chromosome, start position, end position and gene name.

For example: Chr01 27355 28320 Glyma.01G000100 Chr01 58975 67527 Glyma.01G000200 Chr01 67770 69968 Glyma.01G000300 Chr01 90152 95947 Glyma.01G000400 Chr01 90289 91197 Glyma.01G000500