About the dexseq_prepare_annotation.py?

[punching-samuel]

After running this script ‘dexseq_prepare_annotation.py’ with my_GTF, I obtained a file with the dexseq.gff suffix. Should I use this to run dexseq? Before using grsea, I had already run dexseq with my_GTF, it seems that the dexseq result of my previous calculation will give an error after running the script 'grsea.py':

(rmats) bash-4.2$ cat xxx/xxx/7_grsea/myouts/grase1860116.err
multiprocessing.pool.RemoteTraceback:
"""
Traceback (most recent call last):
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
return list(map(*args))
File "xxx/Software/GrASE/grase.py", line 168, in process_gene
g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir = get_gene_files(gene)
File "xxx/Software/GrASE/grase.py", line 78, in get_gene_files
return g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir
UnboundLocalError: local variable 'gff' referenced before assignment
"""

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "xxx/Software/GrASE/grase.py", line 1289, in
main()
File "xxx/Software/GrASE/grase.py", line 1275, in main
p.map(process_gene, genes)
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 367, in map
return self._map_async(func, iterable, mapstar, chunksize).get()
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 774, in get
raise self._value
UnboundLocalError: local variable 'gff' referenced before assignment

Could you help me? I will appreciate it.

[punching-samuel]

But my dexseq.gff is not like the common gtf format, and it can't be read by the makeTxDbFromGFF() function.

[jelardaquino]

Thank you for your question. I might need more information to help. Were you able to create the gene_files directory using creatingFilesByGene.sh before running grase.py? Additionally, are you successful in producing rMATS results using your gtf? Will you also please include an example script that you used to run the code.

When running DEXSeq, dexseq_prepare_annotation.py is used to generate all the exonic parts of your genes and all these exonic parts will be tested by the DEXSeq software. When running dexseq_prepare_annotation.py, make sure to include -r no as an argument. This will ensure that the dexseq.gff output will be in a non-aggregated form (also to ensure compatibility with GrASE). The rest of the DEXSeq pipeline can be completed as normal. Make sure to create a dexseq_output.txt file when the pipeline is finished. GrASE will use this output file in its analysis.

More information on how to use DEXSeq's dexseq_prepare_annotation.py can be found in their tutorial.

[punching-samuel]

Thank you for your question. I might need more information to help. Were you able to create the gene_files directory using creatingFilesByGene.sh before running grase.py? Additionally, are you successful in producing rMATS results using your gtf? Will you also please include an example script that you used to run the code.

When running DEXSeq, dexseq_prepare_annotation.py is used to generate all the exonic parts of your genes and all these exonic parts will be tested by the DEXSeq software. When running dexseq_prepare_annotation.py, make sure to include -r no as an argument. This will ensure that the dexseq.gff output will be in a non-aggregated form (also to ensure compatibility with GrASE). The rest of the DEXSeq pipeline can be completed as normal. Make sure to create a dexseq_output.txt file when the pipeline is finished. GrASE will use this output file in its analysis.

More information on how to use DEXSeq's dexseq_prepare_annotation.py can be found in their tutorial.

I don't know where gff is needed after running dexseq_prepare_annotation.py. I have now run to the step of creatingFilesByGene.sh, and the data structure also meets the requirements of your tutorial.

step1✅
r_script='xxx/Software/GrASE/SplicingGraphs.igraph.R'‘
gff_file= "xxx/correct.orfanage_full_attributes_filter.dexseq.gff'
gtf_file= "xxx/correct.orfanage_full_attributes_filter.gtf"
graphML_dir='/slurm/home/yrd/lulab/lulinrong/yangsen/PROJECT/ys/250207_database_splited/3_ReStart/7_grsea_retry_gff/graphML_dir'
args=(
"$gtf_file"
"Mus musculus" # 作为一个整体参数
"$graphML_dir"
)
Rscript "$r_script" "${args[@]}"
step2✅
sh_file=‘xxx/Software/GrASE/creatingFilesByGene.sh'
rmats_results='xxx/rMATs/myouts-thresds/output_rmats'
py_file=’xxx/DEXSeq/python_scripts/dexseq_prepare_annotation.py'

sh $sh_file -r -s $rmats_results -d $py_file -a $gtf_file -g $graphML_dir -p 8

[punching-samuel]

I encountered another problem when running grsea.py:

"""
Traceback (most recent call last):
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
return list(map(*args))
File "xxx/yangsen/Software/GrASE/grase.py", line 168, in process_gene
g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir = get_gene_files(gene)
File "xxx/yangsen/Software/GrASE/grase.py", line 78, in get_gene_files
return g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir
UnboundLocalError: local variable 'g' referenced before assignment
"""
Could you help me solve this error?

Thank you～

[punching-samuel]

I haven’t had the chance to check the complete script ‘grsea.py’. I only made changes to this part, but the same error still persists.

I encountered another problem when running grsea.py:

""" Traceback (most recent call last): File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 125, in worker result = (True, func(*args, **kwds)) File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar return list(map(*args)) File "xxx/yangsen/Software/GrASE/grase.py", line 168, in process_gene g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir = get_gene_files(gene) File "xxx/yangsen/Software/GrASE/grase.py", line 78, in get_gene_files return g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir UnboundLocalError: local variable 'g' referenced before assignment """ Could you help me solve this error?

Thank you～

[jelardaquino]

The dexseq.gff is needed when running the software DEXSeq in R. DEXSeq requires a flattened annotation file as input.

Here's an exmaple Rscript of how to run DEXSeq:

library(DEXSeq)

#count files had quotes, so removed with sed
countFiles = list.files("/path/to/dexseq/count/files", pattern="counts.txt$", full.names=TRUE)

flattenedFile = list.files("/path/to/dexseq/gff", pattern="gff$", full.names=TRUE)

sampleTable = data.frame(row.names = c("control1", "control2", "control3", "control4", "control5", "case1", "case2", "case3", "case4", "case5"),  condition = c("control", "control", "control", "control", "control", "case", "case", "case", "case", "case"))

dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleTable, design= ~ sample + exon + condition:exon, flattenedfile=flattenedFile)

dxd = estimateSizeFactors( dxd )
dxd = estimateDispersions( dxd )
dxd = testForDEU( dxd )
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults( dxd )
dxr1_df <- as.data.frame(dxr1)
dxr1_df$transcripts <- vapply(dxr1_df$transcripts, paste, collapse = " ", character(1L))
write.table(dxr1_df, "DEXSeq_Results.txt", sep="\t", quote=FALSE, row.names=FALSE)

After running DEXSeq, we use DEXSeq_Results.txt as the input to GrASE:

python grase.py \
                -g grase_results/gene_files \
                -s r \
                --rmats rmats_out \
                --dexseq DEXSeq_Results.txt \
                --nthread 8

Please, let me know if this fixes your issue.

[punching-samuel]

Regarding my first question: Where is the gff generated by the script 'dexseq_prepare_annotation.py' used?
I can run the code:
‘’‘
txdb = makeTxDbFromGFF("xxx.gtf")
flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE )
’‘’
I can get similar results as reading .gff file directly. Did you try the difference between these two methods?

Second question: The data structure I got now meets your tutorial requirement, but an error occurs when running 'grase.py'. Do you have any suggestions?

"""
Traceback (most recent call last):
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 125, in worker
result = (True, func(*args, **kwds))
File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar
return list(map(*args))
File "xxx/yangsen/Software/GrASE/grase.py", line 168, in process_gene
g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir = get_gene_files(gene)
File "xxx/yangsen/Software/GrASE/grase.py", line 78, in get_gene_files
return g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir
UnboundLocalError: local variable 'g' referenced before assignment
"""

Thank you ~

[punching-samuel]

It seems that in a multi-processing environment, the global variable g is not being properly updated across processes, causing it to be referenced before assignment？Just my guess, I'm not sure about that.

Regarding my first question: Where is the gff generated by the script 'dexseq_prepare_annotation.py' used? I can run the code: ‘’‘ txdb = makeTxDbFromGFF("xxx.gtf") flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE ) ’‘’ I can get similar results as reading .gff file directly. Did you try the difference between these two methods?

Second question: The data structure I got now meets your tutorial requirement, but an error occurs when running 'grase.py'. Do you have any suggestions? """ Traceback (most recent call last): File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 125, in worker result = (True, func(*args, **kwds)) File "xxx/Software/anaconda3/envs/rmats/lib/python3.10/multiprocessing/pool.py", line 48, in mapstar return list(map(*args)) File "xxx/yangsen/Software/GrASE/grase.py", line 168, in process_gene g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir = get_gene_files(gene) File "xxx/yangsen/Software/GrASE/grase.py", line 78, in get_gene_files return g, gene, gff, fromGTF_SE, fromGTF_RI, fromGTF_A3SS, fromGTF_A5SS, grase_output_dir UnboundLocalError: local variable 'g' referenced before assignment """

Thank you ~

[jelardaquino]

I just ran the entire GrASE pipeline on Mus musculus. I used these as the references:
gencode.vM36.primary_assembly.annotation.gtf
GRCm39.primary_assembly.genome.fa

The fastq files I used were from this project: PRJNA705172

I also ran DEXSeq using your method:

txdb = makeTxDbFromGFF("xxx.gtf")
flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE )

This is my R script for obtaining the DEXSeq results:

library(DEXSeq)
library(GenomicFeatures)  # for the exonicParts() function
library(txdbmaker)        # for the makeTxDbFromGFF() function
txdb = makeTxDbFromGFF("gencode.vM36.primary_assembly.annotation.gtf")

flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE )
names(flattenedAnnotation) =
  sprintf("%s:E%0.3d", flattenedAnnotation$gene_id, flattenedAnnotation$exonic_part)

bamFiles <- list.files("star_out", pattern="Aligned.sortedByCoord.out.bam$", full.names=TRUE)

se = summarizeOverlaps(
  flattenedAnnotation, bamFiles, singleEnd=FALSE,
  fragments=TRUE, ignore.strand=TRUE )

colData(se)$condition = factor(c("cc", "cc", "cc",  "lung", "lung", "lung"))
dxd = DEXSeqDataSetFromSE( se, design= ~ sample + exon + condition:exon )
dxd = estimateSizeFactors( dxd )
BPPARAM = MulticoreParam(60)
register(BPPARAM)
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
dxd = testForDEU( dxd )
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults( dxd )
dxr1_df <- as.data.frame(dxr1)
dxr1_df$transcripts <- vapply(dxr1_df$transcripts, paste, collapse = " ", character(1L))
write.table(dxr1_df, "DEXSeq_Results_mouse_cc_vs_lung.txt", sep="\t", quote=FALSE, row.names=FALSE)

When running GrASE, there are three steps

1. Creating graphml from Mus musculus using SplicingGraphs.igraph.R
2. Creating gene files using creatingFilesByGene.sh
3. Mapping rMATS and and DEXSeq results using grase.py

For step 1, this is the command I used for the mouse example:
Rscript SplicingGraphs.igraph.R gencode.vM36.primary_assembly.annotation.gtf Mus musculus grase_graphMLs/

For step 2, this is the command I used:
bash creatingFilesByGene.sh -r -s rmats_out/ -d dexseq_prepare_annotation.py -a gencode.vM36.primary_assembly.annotation.gtf -g grase_graphMLs/ -p 8

For step 3, this is the command I used:
python grase.py -g grase_results/gene_files -s r --rmats rmats_out --dexseq DEXSeq_Results_mouse_cc_vs_lung.txt --nthread 8

make sure for step 3, you provide the -g flag with the location of the gene_files directory that was created from step 2.

Please let me know the code you used to run DEXSeq and the command you used for grase.py.

[punching-samuel]

I just ran the entire GrASE pipeline on Mus musculus. I used these as the references: gencode.vM36.primary_assembly.annotation.gtf GRCm39.primary_assembly.genome.fa

The fastq files I used were from this project: PRJNA705172

I also ran DEXSeq using your method:

txdb = makeTxDbFromGFF("xxx.gtf")
flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE )

This is my R script for obtaining the DEXSeq results:

library(DEXSeq)
library(GenomicFeatures)  # for the exonicParts() function
library(txdbmaker)        # for the makeTxDbFromGFF() function
txdb = makeTxDbFromGFF("gencode.vM36.primary_assembly.annotation.gtf")

flattenedAnnotation = exonicParts( txdb, linked.to.single.gene.only=TRUE )
names(flattenedAnnotation) =
  sprintf("%s:E%0.3d", flattenedAnnotation$gene_id, flattenedAnnotation$exonic_part)

bamFiles <- list.files("star_out", pattern="Aligned.sortedByCoord.out.bam$", full.names=TRUE)

se = summarizeOverlaps(
  flattenedAnnotation, bamFiles, singleEnd=FALSE,
  fragments=TRUE, ignore.strand=TRUE )

colData(se)$condition = factor(c("cc", "cc", "cc",  "lung", "lung", "lung"))
dxd = DEXSeqDataSetFromSE( se, design= ~ sample + exon + condition:exon )
dxd = estimateSizeFactors( dxd )
BPPARAM = MulticoreParam(60)
register(BPPARAM)
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
dxd = testForDEU( dxd )
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults( dxd )
dxr1_df <- as.data.frame(dxr1)
dxr1_df$transcripts <- vapply(dxr1_df$transcripts, paste, collapse = " ", character(1L))
write.table(dxr1_df, "DEXSeq_Results_mouse_cc_vs_lung.txt", sep="\t", quote=FALSE, row.names=FALSE)

When running GrASE, there are three steps

1. Creating graphml from Mus musculus using SplicingGraphs.igraph.R
2. Creating gene files using creatingFilesByGene.sh
3. Mapping rMATS and and DEXSeq results using grase.py

For step 1, this is the command I used for the mouse example: Rscript SplicingGraphs.igraph.R gencode.vM36.primary_assembly.annotation.gtf Mus musculus grase_graphMLs/

For step 2, this is the command I used: bash creatingFilesByGene.sh -r -s rmats_out/ -d dexseq_prepare_annotation.py -a gencode.vM36.primary_assembly.annotation.gtf -g grase_graphMLs/ -p 8

For step 3, this is the command I used: python grase.py -g grase_results/gene_files -s r --rmats rmats_out --dexseq DEXSeq_Results_mouse_cc_vs_lung.txt --nthread 8

make sure for step 3, you provide the -g flag with the location of the gene_files directory that was created from step 2.

Please let me know the code you used to run DEXSeq and the command you used for grase.py.

Thank you～ I will retry the proces.

[punching-samuel]

I still have a question. If I discover that due to alternative splicing, the original exon undergoes a frameshift, resulting in a new transcript with a corresponding new exon region, but in DEXSeq, only the exonic regions after flattening are considered within the gene region, is there no way to distinguish between the overlapping old exon and the new exon?

[jelardaquino]

I still have a question. If I discover that due to alternative splicing, the original exon undergoes a frameshift, resulting in a new transcript with a corresponding new exon region, but in DEXSeq, only the exonic regions after flattening are considered within the gene region, is there no way to distinguish between the overlapping old exon and the new exon?

DEXSeq can only test exonic parts that are defined by the gff. If the frameshift and the new exon is part of the annotation (gff), then DEXSeq will be able to annotate that exonic part. If the frameshift and the new exon is not part of the annotation there is no way to identify it without modifying the gff and introducing the new exonic parts. You might have to create a new gff with the new exon to see if it is differentially used between your conditions.