areyesq89 / dexseq Goto Github PK

I have a data set that has counts over our gene of interest using DESeq but no counts at all over the same region using DEXSeq. Do you know how this could be? Is there an internal filter in the DEXSeq call?

Hello all, I'm trying to use DEXseq to analyze usage of regions across different sample. My data is single cell RNA-seq, hence I sum all counts per sample and treat them as bulk RNA-seq data. It's due to DEXseq cannot perform on scRNA data due to memory problem. I additionally add 1 pseudo-count to count matrix (after summation per sample) since size factor estimation uses logarithm of product of counts.

After creating successfully DEXSeqDataSet and estimated size factor, I continue with dispersion estimation. However, It prompted the following error message and I don't know the reason for that

dxd.sct <- estimateDispersions(dxd.sct)
Error in `$<-.data.frame`(`*tmp*`, "dispersion", value = NA) : 
  replacement has 1 row, data has 0

I'm using R 4.1.3 and DEXseq 1.40.0.
Do you guys have any idea about it?

Thanks for any help.

Tried running:

python dexseq_count.py –s yes –p yes –f bam –r pos flat_file.gff alignment_file.bam DEXseq_exon.txt

And got this error:

dexseq_count.py: Error: Please provide three arguments

If I run it with the options at the end it does seem to work:

python dexseq_count.py flat_file.gff alignment_file.bam DEXseq_exon.txt –s yes –p yes –f bam –r pos

Any idea of what is going on? Did I misunderstand something? Thanks so much!

Looks like GenomicFeatures has deprecated disjointExons and will be removing it soon. It's replacement, exonicParts, returns the column "exon_rank" instead of "exonic_part"; "exon_rank" appears to be a transcript-wise exon ranking, instead of gene-wise, and contains multiple values where transcripts share an exon.

Hi,

The current version of the vignette (and the devel version) shows the abstract of the vignette without bolds or italics letters. You might be able to fix it before the next release.

" (\frac{\text{number of transcripts from the gene that contain this exon}} {\text{number of all transcripts from the gene}}).

"

I have experienced issues setting the correct group of samples as the control group. I understand that DEXSeq assigns a control condition based on the alphabetical order of the conditions (i.e. "Mutant" vs "WT", the "Mutant" samples would be set as the control group). I have attempted to correct for this by using 'factor' and/or 'relevel' to set the correct control group, which appears to work when I check the DEXSeq object structure. However when I visualize the results, the plots continue to show the incorrect comparison. Has anyone else experienced this issue? If so, how were you able to correct it?

Hi,

I am writing since I run the DEXSeq workflow to find differentially expressed exons. In the results I noticed that there are some rows in which the coefficients and the corresponding log2Fold changes are set to NA, the counts of the samples seem reasonable (See attached table). What could be the reason for these NA values?

Thanks,
Best,
Simona

Hi,

I used the DEXSeqResults function to generate a result object, transformed it into a dataframe and saved it. I was wondering if there is a way to import this table again as a DEXSeqResults object so that I can run the plotDEXSeq function and skip all the data generation beforehand.

Thanks for your help!

I have the same question as Aufiero: "How the pvalues and pAdj for those exons are then calculated if Log2FC is NA?" Could you clarify?

Originally posted by @pdarakjiank in #10 (comment)

Hi,

I used a dexseq_prepare_annotation.py and dexseq_count.py scripts on the human gtf annotation and RNA-Seq alignments. Some lines look like:
ENSG00000288258+ENSG00000182712+ENSG00000214827:001

Maybe I overlooked in the manual, what these "+" mean. Could you help me? Thank you!

Hi,

I'm running the dexseq_count.py script on paired bulk RNA-seq data, trying to get read counts on an exon level.
The first time I ran the script, I forgot to include the -p yes option for paired end sequencing. However when I ran it a second time, including the argument I got extremely low read counts, compared to the first one.

Here are the results from DEXseq count for some of the exons without -p yes

ENSG00000136531.16+ENSG00000236283.5:032                303
ENSG00000136531.16+ENSG00000236283.5:033                158
ENSG00000136531.16+ENSG00000236283.5:034                490
ENSG00000136531.16+ENSG00000236283.5:035                40
ENSG00000136531.16+ENSG00000236283.5:036                262
ENSG00000136531.16+ENSG00000236283.5:037                188
ENSG00000136531.16+ENSG00000236283.5:038                114

Here are the results from DEXseq count with -p yes

ENSG00000136531.16+ENSG00000236283.5:032                2
ENSG00000136531.16+ENSG00000236283.5:033                2
ENSG00000136531.16+ENSG00000236283.5:034                7
ENSG00000136531.16+ENSG00000236283.5:035                3
ENSG00000136531.16+ENSG00000236283.5:036                2
ENSG00000136531.16+ENSG00000236283.5:037                2
ENSG00000136531.16+ENSG00000236283.5:038                0

Here are the read counts from samtools flagstat for each of the exons above:

samtools flagstat ENSG00000236283.5:032.bam

605 + 0 in total (QC-passed reads + QC-failed reads)
1 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
605 + 0 mapped (100.00% : N/A)
604 + 0 paired in sequencing
308 + 0 read1
296 + 0 read2
604 + 0 properly paired (100.00% : N/A)
604 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:033.bam

274 + 0 in total (QC-passed reads + QC-failed reads)
1 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
274 + 0 mapped (100.00% : N/A)
273 + 0 paired in sequencing
158 + 0 read1
115 + 0 read2
273 + 0 properly paired (100.00% : N/A)
273 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:034.bam

972 + 0 in total (QC-passed reads + QC-failed reads)
1 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
972 + 0 mapped (100.00% : N/A)
971 + 0 paired in sequencing
490 + 0 read1
481 + 0 read2
971 + 0 properly paired (100.00% : N/A)
971 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:035.bam

291 + 0 in total (QC-passed reads + QC-failed reads)
1 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
291 + 0 mapped (100.00% : N/A)
290 + 0 paired in sequencing
156 + 0 read1
134 + 0 read2
290 + 0 properly paired (100.00% : N/A)
290 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:036.bam

709 + 0 in total (QC-passed reads + QC-failed reads)
4 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
709 + 0 mapped (100.00% : N/A)
705 + 0 paired in sequencing
360 + 0 read1
345 + 0 read2
705 + 0 properly paired (100.00% : N/A)
705 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:037.bam

420 + 0 in total (QC-passed reads + QC-failed reads)
2 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
420 + 0 mapped (100.00% : N/A)
418 + 0 paired in sequencing
188 + 0 read1
230 + 0 read2
418 + 0 properly paired (100.00% : N/A)
418 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ENSG00000236283.5:038.bam

276 + 0 in total (QC-passed reads + QC-failed reads)
3 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
276 + 0 mapped (100.00% : N/A)
273 + 0 paired in sequencing
114 + 0 read1
159 + 0 read2
273 + 0 properly paired (100.00% : N/A)
273 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Here are the bedtools intersect results for the same exons:

ENSG00000236283.5:032               605
ENSG00000236283.5:033               274
ENSG00000236283.5:034               972
ENSG00000236283.5:035               291
ENSG00000236283.5:036               709
ENSG00000236283.5:037               420
ENSG00000236283.5:038               276

It looks like the count for aligned R1 reads from samtools are somewhat close to those reported by DEXseq count without ‘-p yes`, but not quite the same. My instinct is just that the reads not counted by dexseq didn't pass the quality threshold. The counts reported by bedtools are the same as the total number of reads aligned reported by samtools.

Any ideas as why we're seeing such low counts with -p yes?

It seems I get an error from reading the htseq-count files with the "DEXSeqDataSetFromHTSeq" function due to the special counts at the end of the file:

"ENSG00000289641":"002"	0
"ENSG00000289642":"001"	16
"ENSG00000289643":"001"	0
"ENSG00000289643":"002"	0
"ENSG00000289644":"001"	0
_ambiguous	369156
_ambiguous_readpair_position	0
_empty	8838909
_lowaqual	0
_notaligned	0

The error I get is:

Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  : 
  line 662256 did not have 3 elements
Calls: DEXSeqDataSetFromHTSeq -> lapply -> lapply -> FUN -> read.table -> scan
Execution halted

The 662256th line matches the "_ambiguous" line in the counts. The counts files have been generated by the "dexseq_count.py" script.
Is the code taking into account these special lines from htseq-count? Should we remove these lines manually?

sessionInfo()
R version 4.0.5 (2021-03-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] DEXSeq_1.36.0               RColorBrewer_1.1-3         
 [3] AnnotationDbi_1.52.0        DESeq2_1.30.1              
 [5] SummarizedExperiment_1.20.0 GenomicRanges_1.42.0       
 [7] GenomeInfoDb_1.26.4         IRanges_2.24.1             
 [9] S4Vectors_0.28.1            MatrixGenerics_1.2.1       
[11] matrixStats_0.61.0          Biobase_2.50.0             
[13] BiocGenerics_0.36.0         BiocParallel_1.24.1        

loaded via a namespace (and not attached):
 [1] httr_1.4.2             bit64_4.0.5            splines_4.0.5         
 [4] assertthat_0.2.1       askpass_1.1            statmod_1.4.36        
 [7] BiocFileCache_1.14.0   blob_1.2.3             Rsamtools_2.6.0       
[10] GenomeInfoDbData_1.2.4 progress_1.2.2         pillar_1.7.0          
[13] RSQLite_2.2.8          lattice_0.20-45        glue_1.6.2            
[16] XVector_0.30.0         colorspace_2.0-3       Matrix_1.4-1          
[19] XML_3.99-0.9           pkgconfig_2.0.3        biomaRt_2.46.3        
[22] genefilter_1.72.1      zlibbioc_1.36.0        purrr_0.3.4           
[25] xtable_1.8-4           scales_1.1.1           tibble_3.1.6          
[28] openssl_2.0.0          annotate_1.68.0        generics_0.1.2        
[31] ggplot2_3.3.5          ellipsis_0.3.2         cachem_1.0.6          
[34] cli_3.2.0              survival_3.3-1         magrittr_2.0.3        
[37] crayon_1.5.1           memoise_2.0.1          fansi_1.0.3           
[40] xml2_1.3.3             hwriter_1.3.2.1        tools_4.0.5           
[43] prettyunits_1.1.1      hms_1.1.1              lifecycle_1.0.1       
[46] stringr_1.4.0          munsell_0.5.0          locfit_1.5-9.4        
[49] DelayedArray_0.16.3    Biostrings_2.58.0      compiler_4.0.5        
[52] rlang_1.0.2            grid_4.0.5             RCurl_1.98-1.6        
[55] rappdirs_0.3.3         bitops_1.0-7           gtable_0.3.0          
[58] curl_4.3.2             DBI_1.1.2              R6_2.5.1              
[61] dplyr_1.0.8            fastmap_1.1.0          bit_4.0.4             
[64] utf8_1.2.2             stringi_1.7.6          Rcpp_1.0.8.3          
[67] vctrs_0.4.0            geneplotter_1.68.0     dbplyr_2.1.1          
[70] tidyselect_1.1.2      

Regards,
Eric

Hi,

I noticed that the exon numbering for the reverse strand genes does not follow the commonly used exon numbering (5' to 3'). Instead, it forces exon numbering according to the forward strand, which can lead to misunderstandings. As an example, the exon numbering of the "FBgn0010909" gene shown in Figure 3, Chapter 6 of the vignette is reversed compared to ensemble.

I have investigated this issue and found two elements of the pipeline that cause this:

1. GenomicFeatures::exonicParts(txdb, linked.to.single.gene.only=TRUE) numbers exons solely according to the forward strand. I managed to manually reorder them with the simple self-written function shown below:

reorder.exons <- function(x){ #changing the order of the exons from the "-" strand
  output <- x  
  temp_tab <- as.data.frame(x)  
  for(i in unique(temp_tab$gene_id)){  
    if(unique(temp_tab$strand[temp_tab$gene_id == i])=="-"){  
      indexes <- which(temp_tab$gene_id == i)  
     output@elementMetadata@listData$exonic_part[indexes] <- rev(output@elementMetadata@listData$exonic_part[indexes])    
    }  
  }  
  return(output)  
}

2. Although I reordered it manually, your DEXSeqDataSetFromSE() wrapper tries to change it back in the further steps, causing an error. I used the DEXSeqDataSet() function directly as shown below and it seems to work fine now.

DEXSeqDataSet(countData = assay(se),
                            sampleData = as.data.frame(colData(se)), 
                            design = ~sample+exon+TREATMENT:exon, 
                            featureID = sprintf("E%3.3d", mcols(se)$exonic_part), 
                            groupID =  as.character(mcols(se)$gene_id),
                            featureRanges = rowRanges(se),
                            transcripts = as.list(mcols(se)$tx_name))

Based on the definition of the DEXSeqDataSetFromSE() function, I assume that the exon numbering according to the forward strand is intentional. Is there a reason for this? Please consider implementing reverse strand support in the next version of the package.

Kind regards,
Tomek

Hi there,

I am trying to use my 42 BAM files (aligned with Subjunc from RSubread package) to prepare a DEXSeqDataSet according to point 2 in the manual (https://bioconductor.org/packages/release/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.html#2_Preparations)

Following would be my script I copied and adjusted from the mentioned link:

txdb = makeTxDbFromGFF("/home/chuddy/bioinformatics/ref_genomes/mouse_38/annotation/ensembl/Mus_musculus.GRCm38.102.gtf.gz")


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("../BAM/subjunc/", full.names = T, include.dirs = F, pattern = "\\.BAM$")
bamFiles <- normalizePath(bamFiles) 

bamFiles = BamFileList(bamFiles)

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

I changed BamFileList(bamFiles) to bamFiles, because we used this function already on the bam files. But also with the original code snippet Rstudio crashed.
Since I have seven condition with each six replicates, I would love to use all BAM files for more complex analyses.
Are there any recommendations how to approach without the need of too much RAM (I have approx 60GB left over)?

Best
Thomas

I created the count files using dexseq_count.py according to the instructions. However, when I used DEXSeqDataSetFromHTSeq, I get the following error:
Error in DEXSeqDataSetFromHTSeq(countFiles, sampleData = sampleTable, :
Count files do not correspond to the flattened annotation file
Here is my code：

setwd("E:/1_科研/1_lncRNA/7_DTU&DEU/differential_exon_usage/HTSeq")
library("DEXSeq")
library("pasilla")

countFiles = list.files("E:/1_科研/1_lncRNA/7_DTU&DEU/differential_exon_usage/HTSeq", pattern=".txt$", full.names=TRUE)
flattenedFile = list.files("E:/1_科研/1_lncRNA/7_DTU&DEU/differential_exon_usage/HTSeq", pattern="gff$", full.names=TRUE)
sampleTable = data.frame(row.names = c( "Hai-05","Him-12","Him-13","Hai-16","Hai-19","Hai-20","Him-35","Him-36"),condition= c("hainanus", "himalayanus", "himalayanus", "hainanus", "hainanus", "hainanus", "himalayanus","himalayanus" ))

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

The following is my script for generating GTF files and txt files：
python /home/students/tools/DEXSeq/inst/python_scripts/dexseq_prepare_annotation.py /home/students/0_database/1_affinis_genemo_hisat2/affinis.genomic.gtf affinis.gff

python /home/students/tools/DEXSeq/inst/python_scripts/dexseq_count.py /home/students/3_qqli/7_lncRNA/nanopore_process/diffe_exon_usage/count/affinis.gff -p yes -f bam /home/students/3_qqli/1_affinis_him_hainan_mRNA/hisat2_featurecounts/min_length_50_hisat2/190505.bam 05.txt

Txt files and GTF files are listed below：
###txt
Parent=evm.model.Contig1.100:001 997
Parent=evm.model.Contig1.100:002 1001
Parent=evm.model.Contig1.100:003 1074
Parent=evm.model.Contig1.100:004 625
Parent=evm.model.Contig1.100:005 1199
Parent=evm.model.Contig1.100:006 1111
Parent=evm.model.Contig1.100:007 1124
Parent=evm.model.Contig1.100:008 952
Parent=evm.model.Contig1.100:009 1413
Parent=evm.model.Contig1.101:001 0
Parent=evm.model.Contig1.101:002 0
Parent=evm.model.Contig1.101:003 0
Parent=evm.model.Contig1.102:001 14
Parent=evm.model.Contig1.102:002 17
Parent=evm.model.Contig1.102:003 16

###gtf
Contig1 dexseq_prepare_annotation.py aggregate_gene 82977 99686 . + . gene_id "Parent=evm.model.Contig1.1"
Contig1 dexseq_prepare_annotation.py exonic_part 82977 82982 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "001"
Contig1 dexseq_prepare_annotation.py exonic_part 84818 84934 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "002"
Contig1 dexseq_prepare_annotation.py exonic_part 85276 85374 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "003"
Contig1 dexseq_prepare_annotation.py exonic_part 85954 86096 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "004"
Contig1 dexseq_prepare_annotation.py exonic_part 86527 86665 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "005"
Contig1 dexseq_prepare_annotation.py exonic_part 90369 90512 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "006"
Contig1 dexseq_prepare_annotation.py exonic_part 91229 91965 . + . gene_id "Parent=evm.model.Contig1.1"; transcripts "Parent=evm.model.Contig1.1"; exonic_part_number "007"

I map the RNA-seq to the genome by hisat2 and produced bam files. The number of "exonic_part_number" in GTF and txt files were consistent, I feel that there is no problem with the files. I wonder if you can help me determine where the problem is, thank you very much!

Hi
I'm having the same error of the author of this post, but my GTF file is not corrupted or whatever, and it has the same number of exonic regions as the number of lines of the count files (644354).
I am also following step by this tutorial.
Here, the code I used:

str(sampleTable)

'data.frame': 9 obs. of 2 variables: $ condition: chr "control"
"control" "control" "control" ... $ libType : chr "paired-end"
"paired-end" "paired-end" "paired-end" ...

flattenedFile

[1] "/Users/admin/GRCh37/genes.gff"

countFiles

[1] "16.EMATO-15.finalCount.txt" "16.EMATO-26.finalCount.txt" [3]
"16.EMATO-27.finalCount.txt" "16.EMATO-28.finalCount.txt" [5]
"16.EMATO-58.finalCount.txt" "16.EMATO-64.finalCount.txt" [7]
"16.GR2108-1.finalCount.txt" "16.GR2108-2.finalCount.txt" [9]
"16.GR2108-3.finalCount.txt"

To perform de analysis:

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

I really can't figure out what's going on!
Thank you in advance,
Giuseppe

Hello, could the DEXseq use customized reference for analysis? Thanks

I am wanting to run a DTU analysis on a differentiation series with three samples iPSC, NSC, and Neuron. I know I can run DEXseq in parallel for each of the three possible pairwise combinations but to make life easiere I just wanted to input everything at once into a DEXseq object and then use something similar like contrasts that is used with DEseq2 to extract the DTU between each of the comparisons.

Is this possible? How would one extract the dxr results for each of the comparisons?

Hello, DEXseq team,

I am currently using DEXseq for alternative splicing analysis. After reading the official documents and my own practical application, I found that when the sample size exceeds 20, the calculation will always stop at the DEXSeq() function. At present, BPPARAM = MultiCoreParam() has been used according to the document and genes with too many 0 values ​​have been deleted.

May I ask, can the BatchJobsParam() mentioned in the documentation for use with large samples provide instructions or examples for adapting to the DEXseq project?

Thanks in advance for any suggestions.