Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • DEXseq - very low numbers of counts

    Hello,

    I wanted to ask the community for help with my problem related to using DEXseq, particularily generating count tables with the Python script. My problem is as follows:

    I have paired-end, strand specific (Illumina TruSeq stranded) RNA-seq data, 101bp read length, multiplexed. The rawdata was aligned using Tophat2 2.0.8 with Bowtie2, no novel junctions. The annotation file used was GRCm38 from Ensembl (mouse) - both the full genome fasta and .gtf annotation file.

    The bam files are around 16GB each, the example .bam I am trying to find solution on has 176609632 reads. Bam files were sorted by name using samtools (samtools sort -n) As next, I flattened the .gtf file using the Python annotation preparation script from DEXSeq which produced a .gff flat. A quick inspection of the formed .gff shows it is "healthy".

    Now the problem begins, even though my dataset is pretty large, most of the exons in the count table have 0 counts. Some of them have few, there are single occurences of exons with 2000 counts, but they are very rare. Naturally, DEXseq returned NA for all the p-values since it is impossible to do anything with just zeroes.

    I loaded the .bam file into IGV to visually inspect if what the count table tells me is true, but it is absolutely not right. As an example of Rbm20 gene, in IGV I can see the exon14 having high coverage whereas it seems almost empty in the count file. As you can see it's not a matter of the script not working at all, it is working only partially.

    The exon 14 of Rbm20 is exon bin 19 for HTSeq, here is the counting part fo this gene:

    ENSMUSG00000043639:001 0
    ENSMUSG00000043639:002 0
    ENSMUSG00000043639:003 0
    ENSMUSG00000043639:004 0
    ENSMUSG00000043639:005 0
    ENSMUSG00000043639:006 2
    ENSMUSG00000043639:007 1
    ENSMUSG00000043639:008 0
    ENSMUSG00000043639:009 0
    ENSMUSG00000043639:010 0
    ENSMUSG00000043639:011 0
    ENSMUSG00000043639:012 4
    ENSMUSG00000043639:013 1
    ENSMUSG00000043639:014 1
    ENSMUSG00000043639:015 1
    ENSMUSG00000043639:016 0
    ENSMUSG00000043639:017 0
    ENSMUSG00000043639:018 0
    ENSMUSG00000043639:019 15

    Now look at that last exon on my IGV screenshot in the attachment.

    And here are small excerpts from my .gtf and my .bam:

    .GTF:

    19 protein_coding exon 53859392 53859513 . + . exon_id "ENSMUSE00001178201"; exon_number "13"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000043639"; gene_name "Rbm20"; p_id "P13408"; transcript_id "ENSMUST00000164202"; transcript_name "Rbm20-201"; tss_id "TSS23511";
    19 protein_coding CDS 53864080 53864187 . + 0 exon_number "14"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000043639"; gene_name "Rbm20"; p_id "P13408"; protein_id "ENSMUSP00000129447"; transcript_id "ENSMUST00000164202"; transcript_name "Rbm20-201"; tss_id "TSS23511";
    19 nonsense_mediated_decay exon 53864080 53867080 . + . exon_id "ENSMUSE00001167670"; exon_number "7"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000043639"; gene_name "Rbm20"; p_id "P15978"; transcript_id "ENSMUST00000161856"; transcript_name "Rbm20-002"; tss_id "TSS56270";
    19 protein_coding exon 53864080 53867080 . + . exon_id "ENSMUSE00001107658"; exon_number "14"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000043639"; gene_name "Rbm20"; p_id "P13408"; transcript_id "ENSMUST00000164202"; transcript_name "Rbm20-201"; tss_id "TSS23511";
    19 protein_coding stop_codon 53864188 53864190 . + 0 exon_number "14"; gene_biotype "protein_coding"; gene_id "ENSMUSG00000043639"; gene_name "Rbm20"; p_id "P13408"; transcript_id "ENSMUST00000164202"; transcript_name "Rbm20-201"; tss_id "TSS23511";

    Part of flattened .gff:

    1 dexseq_prepare_annotation.py aggregate_gene 3054233 3054733 . + . gene_id "ENSMUSG00000090025"
    1 dexseq_prepare_annotation.py exonic_part 3054233 3054733 . + . transcripts "ENSMUST00000160944"; exonic_part_number "001"; gene_id "ENSMUSG00000090025"
    1 dexseq_prepare_annotation.py aggregate_gene 3102016 3102125 . + . gene_id "ENSMUSG00000064842"
    1 dexseq_prepare_annotation.py exonic_part 3102016 3102125 . + . transcripts "ENSMUST00000082908"; exonic_part_number "001"; gene_id "ENSMUSG00000064842"
    1 dexseq_prepare_annotation.py aggregate_gene 3205901 3671498 . - . gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3205901 3206522 . - . transcripts "ENSMUST00000162897"; exonic_part_number "001"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3206523 3207317 . - . transcripts "ENSMUST00000162897+ENSMUST00000159265"; exonic_part_number "002"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3213439 3213608 . - . transcripts "ENSMUST00000159265"; exonic_part_number "003"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3213609 3214481 . - . transcripts "ENSMUST00000162897+ENSMUST00000159265"; exonic_part_number "004"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3214482 3215632 . - . transcripts "ENSMUST00000070533+ENSMUST00000162897+ENSMUST00000159265"; exonic_part_number "005"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3215633 3216344 . - . transcripts "ENSMUST00000162897+ENSMUST00000070533"; exonic_part_number "006"; gene_id "ENSMUSG00000051951"
    1 dexseq_prepare_annotation.py exonic_part 3216345 3216968 . - . transcripts "ENSMUST00000070533"; exonic_part_number "007"; gene_id "ENSMUSG00000051951"


    and my .bam (first few lines):

    DHCDZDN1:177:C31G2ACXX:3:1101:1206:88736 163 MT 6628 50 101M = 6668 141 CAGGAATACCACGACGCTACTCAGACTACCCAGATGCTTACACCACATGAAACACTGTCTCTTCTATAGGATCATTTATTTCACTAACAGCTGTTCTCATC ?;=?D:BBDD3:CEAEEDE?3?<DD?D@DDCCDDII@DEADEBDBDII@@@DCC=CCDDDCACDD@ADDD;;@;ADDEDAA:>>AAAA>A?########## AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:101 YT:Z:UU XS:A:+ NH:i:1
    DHCDZDN1:177:C31G2ACXX:3:1101:1206:88736 83 MT 6668 50 101M = 6628 -141 CACCACATGAAACACTGTCTCTTCTATAGGATCATTTATTTCACTAACAGCTGTTCTCATCATGATNTTTATAATTTGAGAGGCCTTTGCTNCAAAACGAG B?BBABEDA?@>A@CCCCBB@;@;B=CE?3FCIEACCEFFGFFIIIGFFIEFCFBGFEFEEGDB:1#GCFFFIIIIIIFFGGCFEGFC<22#DD1;AB;@? AS:i:-2 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:66C24T9 YT:Z:UU XS:A:+ NH:i:1
    DHCDZDN1:177:C31G2ACXX:3:1101:1207:49279 163 13 55133172 50 101M = 55133267 196 GACCTCCTCAAGAAGATTCTGGGGGAACTAAAAAGGAATTTAGAACTGGTGCAGGTTCTACAGAGAGCAATAGGGTCCCCAGACCCTAACCTGAGAAATTC @@@DBDDEDDHHHEIGIGGIIIIIIBHGHIGIE@GIICGIIIGHIIAFIBCECGGCHGGIIHHGGEE=DCC@C@B>CABD?BDDBDDDCBBDDCDCDDCD@ AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:101 YT:Z:UU XS:A:+ NH:i:1
    DHCDZDN1:177:C31G2ACXX:3:1101:1207:49279 83 13 55133267 50 101M = 55133172 -196 AAATTCCTTCTGGGTTTGTGTCCTGGGATCCACTGCAGTGGGTTTATACATAACTGTGGATGGAGCNGGTCTTCCAGGAAGTATTCCTAAANCTTGTAGAG DCCCDB?DDDDD@DDDDDCDDDFEDDCEGHHEHGGIIGGGGCEIGGIHHICIHJJIEJGHCGFD?1#GHDIIGJJJIGJJJIGGJJIHFA2#HFDD?FC@@ AS:i:-2 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:66A24G9 YT:Z:UU XS:A:+ NH:i:1


    and the tail of the .counts file:

    _ambiguous 147010
    _empty 90546556
    _lowaqual 0
    _notaligned 0


    I searched a bit online, but it seemed like most people having problems with empty count tables had them completely empty. Any suggestions ? Thanks a lot in advance!
    Attached Files
    Last edited by kajot; 04-01-2014, 07:10 AM.

  • #2
    Can you post a reproducible example? The simplest way to do this would be to:

    Code:
    samtools sort -@ 4 namesorted_sample.bam sorted_sample
    samtools index sorted_sample.bam
    samtools view sorted_sample.bam 19:53864080-53867080 > subset.sam
    and then post "subset.sam" somewhere (an attachment here, copy.com, dropbox, where ever). That should correspond to exonic_bin number 19 (at least it does on GRCm38.71). Then we can look at exactly why things are wonky.

    Comment


    • #3
      Hi kajot,

      I have the feeling that could be related to the protocol that you are using. In the most recent protocol of the illumina runs (at least the latest used by the sequencing facility at EMBL), you get strand specific data, but the reads map opposite to the strand where they come from. For example, for a gene that is in the forward strand, you will get all the reads mapping to the reverse strand. This is a it confusing!

      So, what I recommend is that in your IGV browser you could colour your reads depending on the strand where they map to, and then see if what I describe above is the case. If so, I think the parameter "-s reverse" parameter of the dexseq_count.py script should do the trick.

      Alejandro

      Comment


      • #4
        Originally posted by areyes View Post
        ...For example, for a gene that is in the forward strand, you will get all the reads mapping to the reverse strand. This is a [b]it confusing!

        So, what I recommend is that in your IGV browser you could colour your reads depending on the strand where they map to, and then see if what I describe above is the case. If so, I think the parameter "-s reverse" parameter of the dexseq_count.py script should do the trick.

        Alejandro
        More accurately, for paired end reads, read 1 will map to the anti-sense strand while read 2 maps to the sense strand, so in this case you will get reads mapping to both strands. You need to distinguish both mapping strand and read (1 or 2). The attached image from IGV shows mapping of TruSeq Stranded, paired end reads. In this case the gene is transcribed from the forward strand of the reference (left to right). Reads are colored by strand, blue == minus, red == plus (in this highly zoomed in view the reads glyphs also have arrows indicating their direction). The alignments are then grouped by "first-in-pair". The top group is read 1 and the bottom group read 2.

        Alejandro is correct that for TruSeq Stranded Kit libraries you need to specify "-s reverse" for dexseq_count.py or htseq-count if you are counting reads to genes for DESeq(2)
        Attached Files

        Comment


        • #5
          Thank you for all your rapid replies! I was swamped today with regular wet-lab work so I had little time to check all your suggestions, I will try to re-run counting script with -p reverse tomorrow morning.

          I only managed to get the subset of reads mapping to aforementioned exon 14 of Rbm20. This was the exon that was covered by multiple reads in IGV and received only around 50 counts in DEXseq.

          The subset for exon 14 is here:



          I looked at IGV and displayed reads joined together when they are a pair. What I can see is that for Rbm20 which is on + strand, I have read 1 mapped to - strand, and read 2 mapped to + strand, pair orientation F2R1. As I understand this means I have the "Illumina protocol bug" and have to run counting with -s reverse option ?




          --------------------

          Update: I just ran the counting script with -s reverse and it all seems to work. I have now almost 6 million reads that are empty, 1.5 milion ambigous and the rest (around 72.5 million) counted. Since I specified fr-firststrand in Tophat, isn't it somehow affecting my alignment ? Shouldn't it be fr-secondstrand then ? Or am I mixing something up here now ?
          Last edited by kajot; 04-02-2014, 08:54 AM.

          Comment


          • #6
            Originally posted by kajot View Post

            Update: I just ran the counting script with -s reverse and it all seems to work. I have now almost 6 million reads that are empty, 1.5 milion ambigous and the rest (around 72.5 million) counted. Since I specified fr-firststrand in Tophat, isn't it somehow affecting my alignment ? Shouldn't it be fr-secondstrand then ? Or am I mixing something up here now ?
            fr-firststrand is the correct setting to use in Tophat for TruSeq Stranded libraries so you are fine. What 'fr-firststrand' means is the the first read of the pair (or only read for single end reads) matches the 'first strand of cDNA synthesized', which is anti-sense to the RNA.

            Comment


            • #7
              Originally posted by areyes View Post

              I have the feeling that could be related to the protocol that you are using. In the most recent protocol of the illumina runs (at least the latest used by the sequencing facility at EMBL), you get strand specific data, but the reads map opposite to the strand where they come from. For example, for a gene that is in the forward strand, you will get all the reads mapping to the reverse strand. This is a it confusing!

              So, what I recommend is that in your IGV browser you could colour your reads depending on the strand where they map to, and then see if what I describe above is the case. If so, I think the parameter "-s reverse" parameter of the dexseq_count.py script should do the trick.

              Alejandro
              Thank you!

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Strategies for Sequencing Challenging Samples
                by seqadmin


                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                03-22-2024, 06:39 AM
              • seqadmin
                Techniques and Challenges in Conservation Genomics
                by seqadmin



                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                Avian Conservation
                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                03-08-2024, 10:41 AM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 06:37 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, Yesterday, 06:07 PM
              0 responses
              9 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-22-2024, 10:03 AM
              0 responses
              49 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 03-21-2024, 07:32 AM
              0 responses
              67 views
              0 likes
              Last Post seqadmin  
              Working...
              X