Request for strandless data implementation

[jkniehaus]

Hello,
Thanks for the tool! I appreciate having the choice to not extend into an overlapping gene.
I have a BAM from RNA-seq data that is not strand aware but is still 3' enriched.
Would it be possible to choose not to UTR extension if there's a gene on the antisense strand? I imagine this might get tricky if two genes overlap entirely, but this would be desirable in my case to make sure I'm quantifying expression in the appropriate gene.

Thanks!
Jesse

[zolotarovgl]

Hi Jesse! Indeed, the tool isn’t supposed to work with unstranded data, but it should be fairly easy to add this functionality. Can you drop me an email: ***@***.*** ? I’d need some sample of the data to see how it looks like and to try some things out. Thanks, Grisha

…

On Thu, 23 May 2024 at 22:49, Jesse Niehaus ***@***.***> wrote: Hello, Thanks for the tool! I appreciate having the choice to not extend into an overlapping gene. I have a BAM from RNA-seq data that is not strand aware but is still 3' enriched. Would it be possible to choose not to UTR extension if there's a gene on the antisense strand? I imagine this might get tricky if two genes overlap entirely, but this would be desirable in my case to make sure I'm quantifying expression in the appropriate gene. Thanks! Jesse — Reply to this email directly, view it on GitHub <#8>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGZDQJ5NODJKU4NU7NL365TZDZI53AVCNFSM6AAAAABIGL4A46VHI2DSMVQWIX3LMV43ASLTON2WKOZSGMYTGOBRHE4DAOI> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[jkniehaus]

Thanks for the quick response, Grisha!
My email is kylius0@gmail.com
I have can subset a bam to a few million reads generated from STAR. They're paired-end, which might be problematic.

I appreciate it!
Jesse

[jdrnevich]

Was the unstranded implementation ever added and if so, is there a flag for it? I also have unstranded data and do not see anything obvious in the help for it.

$ python /home/groups/hpcbio/apps/geneext/GeneExt/geneext.py -h
      ____                 _____      _
     / ___| ___ _ __   ___| ____|_  _| |_
    | |  _ / _ \ '_ \ / _ \  _| \ \/ / __|
    | |_| |  __/ | | |  __/ |___ >  <| |_
     \____|\___|_| |_|\___|_____/_/\_\\__|

          ______    ___    ______
    -----[______]==[___]==[______]===>----

    Gene model adjustment for improved single-cell RNA-seq data quantification


usage: geneext.py [-h] -g GENOME [-b BAM] [-p PEAKS] [-o OUTPUT] [-m MAXDIST] [-inf INF] [-ouf OUF] [-t T] [-tag TAG] [-v VERBOSE] [-j JOBS]
                  [--output_mode OUTPUT_MODE] [-l] [--clip_strand CLIP_STRAND] [--clip_5prime] [-orphan] [-orphan_maxdist ORPHAN_MAXDIST]
                  [-orphan_maxsize ORPHAN_MAXSIZE] [-peak_perc PEAK_PERC] [-nocluster] [-subsamplebam SUBSAMPLEBAM] [-keep] [--onlyfix]
                  [-force] [-rerun]

Program: GeneExt
Version: 1.0

optional arguments:
  -h, --help            show this help message and exit
  -g GENOME, --genome GENOME
                        Genome .gtf/.gff/.bed file.
  -b BAM, --bam BAM     Input .bam file.
  -p PEAKS, --peaks PEAKS
                        Peaks .bed file. Incompatible with -b.
                        If provided, extension is performed using specified coordinates in .bed format.
                        (Can be seful in cases of FLAM-seq / Nano-3P-seq data or when manual filtering of the peaks is needed.)
  -o OUTPUT, --output OUTPUT
                        Output annotation file.


                        ================ Additional arguments ================
  -m MAXDIST, --maxdist MAXDIST
                        Maximal distance for gene extension.
                        If not set, a median length of gene (genomic span!) is used.
  -inf INF              Input genes file format, if None, will be guessed from a file extension.
  -ouf OUF              Output file format, if not given, will be guessed from a file extension.
  -t T                  Temporary directory. [tmp_{output_file_prefix}]
  -tag TAG              Tag to be added to the fake gene source and IDs so these can be easily identified downstream. [GeneExt]
  -v VERBOSE, --verbose VERBOSE
                        Verbosity level. [0],1,2,3
  -j JOBS, --jobs JOBS  Number of parallel cores. [1]
  --output_mode OUTPUT_MODE
                        How to extend the gene (only for .gff/.gtf files) [new_transcript]
                                * new_transcript - creates a new transcript feature with the last exon extended
                                * new_exon - creates an extended last exon
  -l, --longest         Whether to select the longest isoform per gene.
  --clip_strand CLIP_STRAND
                        How to treat gene extension overlaps.
                        sense - default,restrict overlaps into downstream genes on the same strand
                        both - restrict overlaps regardless of the strand.
  --clip_5prime         Use this to clip overlaps between genes. The downstream gene will be clipped.
                        CAVE: use carefully if the genome contains many overlapping genes (e.g. mitochondrial genomes).

                        ================ Orphan peaks ================
                        Orpan peaks are peaks that could not be assigned to any gene due to the distance.
                        Since those can correspond to missing 3prime UTRs or missing genes, we suggest keeping them.
                        Orphan peaks are merged into orphan peak clusters by distance and then are kept as if they were genes.

  -orphan, --orphan     Whether to add orphan peaks
  -orphan_maxdist ORPHAN_MAXDIST, --orphan_maxdist ORPHAN_MAXDIST
                        Orphan peak merging: Maximum distance between orphan peaks to merge. [75-th percentile of intron length, bp]
  -orphan_maxsize ORPHAN_MAXSIZE, --orphan_maxsize ORPHAN_MAXSIZE
                        Orphan peak merging: Maximum size of an orphan peak cluster. Default: [median gene length, bp]
  -peak_perc PEAK_PERC, --peak_perc PEAK_PERC
                        Coverage threshold (percentile of macs2 genic peaks coverage). [1-99, 25 by default].
                        All peaks called with macs2 are required to have a coverage AT LEAST as N-th percentile of the peaks falling within genic regions.
                        This parameter allows to filter out the peaks based on the coverage BEFORE gene extension.
  -nocluster, --nocluster
                        Do not merge orphan peaks based on distance.


                        ================ Miscellaneous ================
  -subsamplebam SUBSAMPLEBAM, --subsamplebam SUBSAMPLEBAM
                        If set, will subsample bam to N reads before the peak calling. Useful for large datasets. Bam file should be indexed.
                        Default: None
  -keep, --keep_intermediate_files
                        Use this to keep .bam and other temporary files in the a temporary directory. Useful for troubleshooting.
  --onlyfix             If set, GeneExt will only try to fix the annotation, no extension is performed.
  -force, --force       If set, GeneExt will ignore previously computed files and will re-run everythng from scratch.
  -rerun, --rerun       If set, GeneExt will try to re-run the analysis using previously computed files it finds in the temporary directory (subsampled data, MACS2 results, peaks with coverage etc.).

[zolotarovgl]

Dear @jdrnevich, I haven't added the unstranded functionality yet but I will try to it ASAP. Out of curiosity, which kind of single-cell RNA-seq protocol generates an unstranded data?

[jdrnevich]

It is actually bulk, not single cell, using Takara's SMART-Seq mRNA library prep kit for very low input mRNA. These are coral larvae samples and the genome/gtf are definitely poorly annotated. Tons of primer-dimers and the need for aggressive trimming. Of those left, unique alignments were decent 70-85% but of those over 75% are intergenic. I've tried running GeneExt on one giant bam of all combined samples, but I got a warning about "Too many orphan peaks: 28,201 orphan peaks/clusters vs 24,564 input genes." I did count using the extended + orphan gtf file, but the intergenic rate barely dropped, and all I got was a slight increase in the ambiguous rate. It is very likely there are more issues than just poor annotation, as the SMART-Seq kit is known to grab anything remotely RNA-looking. Just trying to salvage anything I can!