Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • sheibani
    replied
    ---------------------------------------------------
    Last edited by sheibani; 07-30-2015, 02:02 AM.

    Leave a comment:


  • dpryan
    replied
    I posted an example R script here. This was originally written in the context of making an input file for CQN, so you can comment out all of the GC% related lines if that really matters to you.

    Leave a comment:


  • dpryan
    replied
    I would say yes, the flatter the coverage the more likely it is to be real.

    Regarding RPKMs, you just need the gene (or other feature) lengths. I'm pretty sure I've posted a script to take a GTF and output union gene model lengths before. I'll have to look around for it (if nothing else, I'll just post it to github). You can then divide the normalized counts by those values and then divide by a million and you'll have RPKM (or FPKM if you used paired-end reads).

    Leave a comment:


  • kajot
    replied
    So the more "flat" coverage I get towards the center of the intron, the more reliable indication this is that I am dealing with a legitimate intron retention ? Is my understanding right ?

    And related to this, are you aware of any clean way of turning the HTSeq counts into RPKM ? I could write a script myself, but why re-invent the wheel.
    If there is nothing I could use then what would be the best value to normalize (the "M" part in RPKM) ? Million of aligned pairs ? Million of "counts" from HTSeq ?


    Have a great weekend!

    Piotr

    Leave a comment:


  • dpryan
    replied
    So you basically plotted an averaged RPKM variant versus fold-change. Statistical power ends up coming from the absolute number of counts in each sample and is independent of the length of an intron (or exon or gene). I realize that this may seem somewhat counter-intuitive, but that's the way the math works.

    One things to look at is the distribution of reads in an intron. If they tend to be clumpy or mostly clustered near the exon:intron bounds, then I would guess that they're not accurately representing actual exon inclusion events. I don't know of a simple metric for that, though I think RNAseqQC has a coverage metrics of some sort that can indicate bias...maybe that'd be useful in this case.

    Leave a comment:


  • kajot
    replied
    I am not sure if my quick coverage look-up for those introns makes sense, maybe this would be important - I took an average of counts from all conditions for any given intron and divided this value by the width of the intron (all based on values in DEXSeq output).
    I think it gives a rough estimate about the abundance of reads for a given event ?

    Below is a plot of those events. There is a cutoff for abs(log2FC) of 0.5. The y-axis is the value I just mentioned (rough estimate of coverage).

    Leave a comment:


  • dpryan
    replied
    What sort of fold-changes are you seeing for those questionable cases?

    Leave a comment:


  • kajot
    replied
    Hi

    I tried the script from Devon, thanks AGAIN for help It works nicely and I managed to run DEXSeq with the new intronic annotation. I decided to delete exons and rename introns to exonic_part, to fool DEXSeq into thinking that we are dealing with a normal everyday situation...

    I am looking at the result table, filtered with padj < 0.01 and discarding all entries with NA padj or log2FC. I am still dealing with around 270 reported hits. However, after inspecting some of those hits visually in IGV and also by looking at read count, I can see that the majority of events comes from low-coverage regions with most of the events having 0 - 100 counts.

    The dataset I am looking at doesn't have a staggering depth, it's a mammalian tissue with ~45M paired-end reads on average per sample (polyA selected) so simply looking at intron retention here might be not the most exciting thing to do, but it just makes me wonder:

    How prone is DEXSeq to report an intron as a hit (differential usage) if the read counts are overall low ? Doesn't the FDR go very high with the decreasing sequencing depth in case of splicing/intron retention analysis ? Since DEXSeq is fitting the model based on dispersion which is also affected by read count (as far as I understand it), is it also taking the low coverage into account while reporting the p.adjusted values ?
    Last edited by kajot; 11-20-2014, 09:00 AM. Reason: added RNAseq prep detail

    Leave a comment:


  • xguo
    replied
    Hi, Alejandro,
    Thanks for posting the perl script. It is simple and works great. Can we still generate visualization graph with plotDEXSeq if we use the annotation file with both exons and introns?

    thanks
    Xiang

    Leave a comment:


  • areyes
    replied
    The perl version of @dpryan's script (be careful, its perl!, :P)

    open(FILE, "Homo_sapiens.GRCh37.68.DEXSeq.gtf");
    my $previousGene = "";
    while(<FILE>){
    next if $_ =~ /aggregate\_gene/;
    $_ =~ /gene_id \"(\S+)\"/;
    $currentGene = $1;
    @lineinfo = split( /\t/, $_ );
    $currentStart = $lineinfo[3];
    $currentEnd = $lineinfo [4];
    $_ =~ /transcripts \"(\S+)\"/;
    $transcripts = $1;
    $_ =~ /exonic\_part\_number \"(\S+)\"/;
    $exonPart = $1;
    $_ =~ /gene\_id \"(\S+)\"/;
    $geneID = $1;
    if( $previousGene eq $currentGene ){
    if( $currentStart - $previousEnd > 1 ){
    $exonPart = $exonPart - 1;
    $exonPart = sprintf( "%3.3d", $exonPart );
    $nPart = $exonPart."i";
    $end = $currentStart - 1;
    $start = $previousEnd + 1;
    print "$lineinfo[0]\t$lineinfo[1]\t$lineinfo[2]\t$start\t$end\t.\t$lineinfo[6]\t.\ttranscripts \"$transcripts\"; exonic_part_number \"$nPart\"; gene_id \"$geneID\"\n";
    }
    }
    print $_;
    $previousGene = $currentGene;
    $previousEnd = $currentEnd;
    }
    close(FILE);

    Leave a comment:


  • xguo
    replied
    Dear Devon,
    Thanks a lot for sharing your code. I'll give it a try.

    regards
    -Xiang

    Leave a comment:


  • dpryan
    replied
    Here's an example script I wrote in R that takes an annotation GFF file prepared by DEXSeq and adds "intronic_part" records with associated gene_ids. You'll find both the original exonic bins and the newer intronic bins in the resulting file (you'd need to change the file names), which you may or may not want. You could also modify this script to make the intronic bins "exonic bins" and just look at exons and introns at once (I don't know how well that'd work) or remove the old exonic bins and simply label the intronic_parts exonic_parts, which would probably make your life easier in DEXSeq (I just modified DEXSeq). For the latter, just change the "asGFF2()" function.

    Code:
    #!/usr/bin/Rscript
    library(GenomicRanges)
    library(rtracklayer)
    library(parallel)
    
    gff <- import.gff2("Mus_musculus.GRCm38.71.DEXSeq.gff", asRangedData=F)
    #Add a level for "intronic_part"
    elementMetadata(gff)$type <- factor(elementMetadata(gff)$type, levels=c(levels(elementMetadata(gff)$type), "intronic_part"))
    #Fix the exonic_part_number to be a properly formatted character
    USE <- which(!is.na(elementMetadata(gff)$exonic_part_number))
    exonic_parts <- sprintf("%03i", elementMetadata(gff)$exonic_part_number[USE])
    elementMetadata(gff)$exonic_part_number <- as.character(elementMetadata(gff)$exonic_part_number)
    elementMetadata(gff)$exonic_part_number[USE] <- exonic_parts
    
    #Split by gene_id
    grl <- split(gff, elementMetadata(gff)$gene_id)
    #Add the introns as "intronic_parts
    add_introns <- function(gr) {
        exons <- gr[which(elementMetadata(gr)$type=="exonic_part"),]
        if(length(exons) > 1) {
            seqname <- seqnames(exons)[-1]
            starts <- end(exons)+1
            starts <- starts[-length(starts)]
            ends <- start(exons)-1
            ends <- ends[-1]
            bounds <- IRanges(start=starts, end=ends)
            strand <- strand(exons)[-1]
            introns <- GRanges(seqnames=seqname, ranges=bounds, strand=strand)
            intron_ids <- sprintf("%03i", c(1:length(introns)))
            #Remove 0-width introns
            DISCARD <- which(width(introns) <= 0)
            if(length(DISCARD) > 0) {
                introns <- introns[-DISCARD]
                intron_ids <- intron_ids[-DISCARD] #Set intron numbers so they follow their respective exonic parts
            }
            if(length(introns) > 0) {
                #create the meta-data
                df <- as.data.frame(elementMetadata(exons))
                nrows <- length(introns)
                metadf <- df[1:nrows,] #does this need to deal with gene_id and transcripts differently?
                metadf <- transform(metadf, gene_id=as.character(gene_id), transcripts=as.character(transcripts))
                metadf$transcripts <- as.character(c(rep(NA, nrows)))
                metadf$type <- factor(c(rep("intronic_part", nrows)), levels=levels(metadf$type))
                metadf$exonic_part_number <- intron_ids
                elementMetadata(introns) <- metadf
                #Merge the GRanges
                gr <- append(gr, introns)
                gr <- gr[order(start(gr), elementMetadata(gr)$type),] #resort
            }
        }
        return(gr)
    }
    with_introns <- endoapply(grl, add_introns) 
    #reorder things
    chroms <- sapply(with_introns, function(x) as.factor(seqnames(x))[1])
    starts <- sapply(with_introns, function(x) start(x)[1])
    o <- order(chroms, starts)
    with_introns2 <- with_introns[o]
    ##Merge into a GRange
    #with_introns2 <- unlist(with_introns2, use.names=F, recursive=T)
    #Create GFF formatted output
    asGFF2 <- function(x) {
        df <- as.data.frame(x)
        aggregates <- which(df$type == "aggregate_gene")
        meta <- character(nrow(df))
        meta[aggregates] <- sprintf("gene_id \"%s\"", df$gene_id[aggregates])
        #This gives introns a transcript "NA" field, which may not be ideal
        meta[-aggregates] <- sprintf("transcripts \"%s\"; exonic_part_number \"%s\"; gene_id \"%s\"", df$transcripts[-aggregates], df$exonic_part_number[-aggregates], df$gene_id[-aggregates])
        paste(df$seqnames, "dexseq_prepare_annotation.py", df$type, df$start, df$end, ".", df$strand, ".", meta, sep="\t")
    }
    outputGFF <- unlist(lapply(with_introns2, asGFF2))
    write.table(outputGFF, file="Mus_musculus.GRCm38.71.DEXSeq.introns.gff", row.names=F, col.names=F, quote=F)

    Leave a comment:


  • xguo
    started a topic DEXSeq for intron retention

    DEXSeq for intron retention

    Dear all,
    I'm trying to use DEXSeq to detect intron retention events in diseased human samples. Essentially we need to count how many reads map to each intron and then use DEXSeq to do statistical comparison. My question is how we shall prepare the intron annotation and count files for DEXSeq. One way is to retrieve intron gtf file from UCSC table browser, and collapse them into intronic counting bins. The other way is to extract intron coordinates based on exon gtf file prepared by dexseq_prepare_annotation.py included in DEXSeq.

    Does anyone have experience in identification of differential intron retention using DEXSeq? It'll be great if you can share your insights in the preparation of intron count files.

    thanks
    Xiang

Latest Articles

Collapse

  • seqadmin
    Recent Advances in Sequencing Technologies
    by seqadmin



    Innovations in next-generation sequencing technologies and techniques are driving more precise and comprehensive exploration of complex biological systems. Current advancements include improved accessibility for long-read sequencing and significant progress in single-cell and 3D genomics. This article explores some of the most impactful developments in the field over the past year.

    Long-Read Sequencing
    Long-read sequencing has seen remarkable advancements,...
    12-02-2024, 01:49 PM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 07:45 AM
0 responses
8 views
0 likes
Last Post seqadmin  
Started by seqadmin, Yesterday, 07:59 AM
0 responses
11 views
0 likes
Last Post seqadmin  
Started by seqadmin, 12-09-2024, 08:22 AM
0 responses
9 views
0 likes
Last Post seqadmin  
Started by seqadmin, 12-02-2024, 09:29 AM
0 responses
175 views
0 likes
Last Post seqadmin  
Working...
X