Answers.Body: <p>Given the high level of multimapping in this region, you'll need to use the <code>-M --primary</code> options if you want to keep many of the alignments. I would be very hesitant to use these numbers as input for DESeq2 or similar programs, since it's fairly questionable whether one should fully trust the "randomness" of the aligner's assignments. I'm more comfortable using something like <a href="http://salmon.readthedocs.io/en/latest/" rel="nofollow noreferrer">salmon</a> for such cases.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-19T21:16:01.000Z
Answers.LastActivityDate: 2017-07-19T21:16:01.000Z
Answers.UserDisplayName: Devon Ryan
Body: <p>I have surprisingly low counts when running <code>featureCounts</code> on some (single-end) RNA-seq data mapped on <em>C. elegans</em> genome using <code>hisat2</code>.</p>
  
  <p>To more easily show the problem, I generated a small subset of the bam file and of the annotation file I'm using. Here is what I can see when loading these two files on IGV:</p>
  
  <p><a href="https://i.sstatic.net/iijez.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/iijez.png" alt="View of the his-16 -&gt; his-12 region in IGV"></a></p>
  
  <p>The coverage for the "his-11" and "his-15" genes peak at around 3600. For "his-12" and "his-16", it peaks at more than 1000</p>
  
  <p>Here are the content of my annotation file:</p>
  
  <pre><code>$ cat his_11-16.gtf 
  II  ensembl transcript  13817759    13818143    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001890"; gene_name "his-16"; gene_source "ensembl"; gene_version "1"; p_id "P8354"; transcript_biotype "protein_coding"; transcript_id "ZK131.10"; transcript_name "ZK131.10"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS11539";
  II  ensembl transcript  13818371    13818739    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001889"; gene_name "his-15"; gene_source "ensembl"; gene_version "1"; p_id "P8185"; transcript_biotype "protein_coding"; transcript_id "ZK131.9"; transcript_name "ZK131.9"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52436";
  II  ensembl transcript  13819626    13819937    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001888"; gene_name "his-14"; gene_source "ensembl"; gene_version "1"; p_id "P21890"; transcript_biotype "protein_coding"; transcript_id "ZK131.8"; transcript_name "ZK131.8"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS50308";
  II  ensembl transcript  13820210    13820620    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001887"; gene_name "his-13"; gene_source "ensembl"; gene_version "1"; p_id "P2149"; transcript_biotype "protein_coding"; transcript_id "ZK131.7"; transcript_name "ZK131.7"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS35106";
  II  ensembl transcript  13821198    13821581    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001886"; gene_name "his-12"; gene_source "ensembl"; gene_version "1"; p_id "P26082"; transcript_biotype "protein_coding"; transcript_id "ZK131.6"; transcript_name "ZK131.6"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS23693";
  II  ensembl transcript  13821778    13822314    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001885"; gene_name "his-11"; gene_source "ensembl"; gene_version "1"; p_id "P10535"; transcript_biotype "protein_coding"; transcript_id "ZK131.5"; transcript_name "ZK131.5"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52792";
  </code></pre>
  
  <p>The <code>featureCounts</code> run:</p>
  
  <pre><code>$ featureCounts -a his_11-16.gtf -o his_11-16_counts.txt -t transcript -g gene_name -O his_11-16_sorted.bam
  
          ==========     _____ _    _ ____  _____  ______          _____  
          =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
            =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
              ====      \___ \| |  | |  _ &lt;|  _  /|  __|   / /\ \ | |  | |
                ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
          ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
        v1.5.2
  
  //========================== featureCounts setting ===========================\\
  ||                                                                            ||
  ||             Input files : 1 BAM file                                       ||
  ||                           S his_11-16_sorted.bam                           ||
  ||                                                                            ||
  ||             Output file : his_11-16_counts.txt                             ||
  ||                 Summary : his_11-16_counts.txt.summary                     ||
  ||              Annotation : his_11-16.gtf (GTF)                              ||
  ||      Dir for temp files : ./                                               ||
  ||                                                                            ||
  ||                 Threads : 1                                                ||
  ||                   Level : meta-feature level                               ||
  ||              Paired-end : no                                               ||
  ||         Strand specific : no                                               ||
  ||      Multimapping reads : not counted                                      ||
  || Multi-overlapping reads : counted                                          ||
  ||   Min overlapping bases : 1                                                ||
  ||                                                                            ||
  \\===================== http://subread.sourceforge.net/ ======================//
  
  //================================= Running ==================================\\
  ||                                                                            ||
  || Load annotation file his_11-16.gtf ...                                     ||
  ||    Features : 6                                                            ||
  ||    Meta-features : 6                                                       ||
  ||    Chromosomes/contigs : 1                                                 ||
  ||                                                                            ||
  || Process BAM file his_11-16_sorted.bam...                                   ||
  ||    Single-end reads are included.                                          ||
  ||    Assign reads to features...                                             ||
  ||    Total reads : 35037                                                     ||
  ||    Successfully assigned reads : 1849 (5.3%)                               ||
  ||    Running time : 0.00 minutes                                             ||
  ||                                                                            ||
  ||                         Read assignment finished.                          ||
  ||                                                                            ||
  || Summary of counting results can be found in file "his_11-16_counts.txt"    ||
  ||                                                                            ||
  \\===================== http://subread.sourceforge.net/ ======================//
  </code></pre>
  
  <p>And the resulting counts file:</p>
  
  <pre><code>$ cat his_11-16_counts.txt
  # Program:featureCounts v1.5.2; Command:"featureCounts" "-a" "his_11-16.gtf" "-o" "his_11-16_counts.txt" "-t" "transcript" "-g" "gene_name" "-O" "his_11-16_sorted.bam" 
  Geneid  Chr Start   End Strand  Length  his_11-16_sorted.bam
  his-16  II  13817759    13818143    -   385 869
  his-15  II  13818371    13818739    +   369 3
  his-14  II  13819626    13819937    -   312 953
  his-13  II  13820210    13820620    +   411 23
  his-12  II  13821198    13821581    -   384 423
  his-11  II  13821778    13822314    +   537 1
  </code></pre>
  
  <p><strong>Why such a low assignment rate, why so little counts?</strong></p>
  
  <p>I see that <code>hisat2</code> decided to split an important proportion of the reads between orthologs, like "his-16" and "his-12", but my understanding is that the default values for the <code>minOverlap</code> and <code>fracOverlap</code> parameters should ensure that such split reads will be counted:</p>
  
  <pre><code>  --minOverlap &lt;int&gt;  Minimum number of overlapping bases in a read that is
                        required for read assignment. 1 by default. Number of
                        overlapping bases is counted from both reads if paired
                        end. If a negative value is provided, then a gap of up
                        to specified size will be allowed between read and the
                        feature that the read is assigned to.
  
    --fracOverlap &lt;float&gt; Minimum fraction of overlapping bases in a read that is
                        required for read assignment. Value should be within range
                        [0,1]. 0 by default. Number of overlapping bases is
                        counted from both reads if paired end. Both this option
                        and '--minOverlap' option need to be satisfied for read
                        assignment.
  </code></pre>
  
  <p>What did I get wrong?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-19T16:44:37.167Z
Comments.Text: How many of those are labeled as multimappers? Those are excluded by default.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-19T17:04:22.747Z
Comments.Text: About 62% of the reads appear to have the "secondary alignment" flag bit, and almost 95% have the value for the NH tag higher than 1.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-19T17:44:26.967Z
Comments.Text: well there's the cause.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-20T11:45:36.420Z
Comments.Text: From the screenshot it seems that his-11 and his-15 have *virtually identical* (not merely similar) peak shapes. That should give pause.
Comments.UserDisplayName: Konrad Rudolph
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-20T11:54:16.223Z
Comments.Text: @KonradRudolph To me, this indicates that the mapper does a "good" random assignment of reads when multiple equally-scoring mapping locations are possible (his-15 and his-11 have identical sequences). I'm more worried about the reads that are split between two his-* genes. This is most likely wrong.
Comments.UserDisplayName: bli
Tags: <read-mapping><coverage>
Title: Counts obtained by featureCounts seem much less than observed coverage
ThreadText: # Counts obtained by featureCounts seem much less than observed coverage
  
  (Asked by: username_0 on 2017-07-19)
  
  I have surprisingly low counts when running `featureCounts` on some (single-end) RNA-seq data mapped on *C. elegans* genome using `hisat2`.
  
  To more easily show the problem, I generated a small subset of the bam file and of the annotation file I'm using. Here is what I can see when loading these two files on IGV:
  
  [![View of the his-16 -> his-12 region in IGV](https://i.sstatic.net/iijez.png)](https://i.sstatic.net/iijez.png)
  
  The coverage for the "his-11" and "his-15" genes peak at around 3600. For "his-12" and "his-16", it peaks at more than 1000
  
  Here are the content of my annotation file:
  
  ```
  $ cat his_11-16.gtf 
  II  ensembl transcript  13817759    13818143    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001890"; gene_name "his-16"; gene_source "ensembl"; gene_version "1"; p_id "P8354"; transcript_biotype "protein_coding"; transcript_id "ZK131.10"; transcript_name "ZK131.10"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS11539";
  II  ensembl transcript  13818371    13818739    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001889"; gene_name "his-15"; gene_source "ensembl"; gene_version "1"; p_id "P8185"; transcript_biotype "protein_coding"; transcript_id "ZK131.9"; transcript_name "ZK131.9"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52436";
  II  ensembl transcript  13819626    13819937    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001888"; gene_name "his-14"; gene_source "ensembl"; gene_version "1"; p_id "P21890"; transcript_biotype "protein_coding"; transcript_id "ZK131.8"; transcript_name "ZK131.8"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS50308";
  II  ensembl transcript  13820210    13820620    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001887"; gene_name "his-13"; gene_source "ensembl"; gene_version "1"; p_id "P2149"; transcript_biotype "protein_coding"; transcript_id "ZK131.7"; transcript_name "ZK131.7"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS35106";
  II  ensembl transcript  13821198    13821581    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001886"; gene_name "his-12"; gene_source "ensembl"; gene_version "1"; p_id "P26082"; transcript_biotype "protein_coding"; transcript_id "ZK131.6"; transcript_name "ZK131.6"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS23693";
  II  ensembl transcript  13821778    13822314    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001885"; gene_name "his-11"; gene_source "ensembl"; gene_version "1"; p_id "P10535"; transcript_biotype "protein_coding"; transcript_id "ZK131.5"; transcript_name "ZK131.5"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52792";
  
  ```
  
  The `featureCounts` run:
  
  ```
  $ featureCounts -a his_11-16.gtf -o his_11-16_counts.txt -t transcript -g gene_name -O his_11-16_sorted.bam
  
          ==========     _____ _    _ ____  _____  ______          _____  
          =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
            =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
              ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
                ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
          ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
        v1.5.2
  
  //========================== featureCounts setting ===========================\\
  ||                                                                            ||
  ||             Input files : 1 BAM file                                       ||
  ||                           S his_11-16_sorted.bam                           ||
  ||                                                                            ||
  ||             Output file : his_11-16_counts.txt                             ||
  ||                 Summary : his_11-16_counts.txt.summary                     ||
  ||              Annotation : his_11-16.gtf (GTF)                              ||
  ||      Dir for temp files : ./                                               ||
  ||                                                                            ||
  ||                 Threads : 1                                                ||
  ||                   Level : meta-feature level                               ||
  ||              Paired-end : no                                               ||
  ||         Strand specific : no                                               ||
  ||      Multimapping reads : not counted                                      ||
  || Multi-overlapping reads : counted                                          ||
  ||   Min overlapping bases : 1                                                ||
  ||                                                                            ||
  \\===================== http://subread.sourceforge.net/ ======================//
  
  //================================= Running ==================================\\
  ||                                                                            ||
  || Load annotation file his_11-16.gtf ...                                     ||
  ||    Features : 6                                                            ||
  ||    Meta-features : 6                                                       ||
  ||    Chromosomes/contigs : 1                                                 ||
  ||                                                                            ||
  || Process BAM file his_11-16_sorted.bam...                                   ||
  ||    Single-end reads are included.                                          ||
  ||    Assign reads to features...                                             ||
  ||    Total reads : 35037                                                     ||
  ||    Successfully assigned reads : 1849 (5.3%)                               ||
  ||    Running time : 0.00 minutes                                             ||
  ||                                                                            ||
  ||                         Read assignment finished.                          ||
  ||                                                                            ||
  || Summary of counting results can be found in file "his_11-16_counts.txt"    ||
  ||                                                                            ||
  \\===================== http://subread.sourceforge.net/ ======================//
  
  ```
  
  And the resulting counts file:
  
  ```
  $ cat his_11-16_counts.txt
  # Program:featureCounts v1.5.2; Command:"featureCounts" "-a" "his_11-16.gtf" "-o" "his_11-16_counts.txt" "-t" "transcript" "-g" "gene_name" "-O" "his_11-16_sorted.bam" 
  Geneid  Chr Start   End Strand  Length  his_11-16_sorted.bam
  his-16  II  13817759    13818143    -   385 869
  his-15  II  13818371    13818739    +   369 3
  his-14  II  13819626    13819937    -   312 953
  his-13  II  13820210    13820620    +   411 23
  his-12  II  13821198    13821581    -   384 423
  his-11  II  13821778    13822314    +   537 1
  
  ```
  
  **Why such a low assignment rate, why so little counts?**
  
  I see that `hisat2` decided to split an important proportion of the reads between orthologs, like "his-16" and "his-12", but my understanding is that the default values for the `minOverlap` and `fracOverlap` parameters should ensure that such split reads will be counted:
  
  ```
    --minOverlap <int>  Minimum number of overlapping bases in a read that is
                        required for read assignment. 1 by default. Number of
                        overlapping bases is counted from both reads if paired
                        end. If a negative value is provided, then a gap of up
                        to specified size will be allowed between read and the
                        feature that the read is assigned to.
  
    --fracOverlap <float> Minimum fraction of overlapping bases in a read that is
                        required for read assignment. Value should be within range
                        [0,1]. 0 by default. Number of overlapping bases is
                        counted from both reads if paired end. Both this option
                        and '--minOverlap' option need to be satisfied for read
                        assignment.
  
  ```
  
  What did I get wrong?
  
  **Tags:** read-mapping, coverage
  
  > **Comments:**
  > **username_1**: How many of those are labeled as multimappers? Those are excluded by default.  
  > **username_0**: About 62% of the reads appear to have the "secondary alignment" flag bit, and almost 95% have the value for the NH tag higher than 1.  
  > **username_1**: well there's the cause.  
  > **username_2**: From the screenshot it seems that his-11 and his-15 have *virtually identical* (not merely similar) peak shapes. That should give pause.  
  > **username_0**: @username_2 To me, this indicates that the mapper does a "good" random assignment of reads when multiple equally-scoring mapping locations are possible (his-15 and his-11 have identical sequences). I'm more worried about the reads that are split between two his-* genes. This is most likely wrong.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2017-07-19. Score: 4 | Accepted answer
  
  Given the high level of multimapping in this region, you'll need to use the `-M --primary` options if you want to keep many of the alignments. I would be very hesitant to use these numbers as input for DESeq2 or similar programs, since it's fairly questionable whether one should fully trust the "randomness" of the aligner's assignments. I'm more comfortable using something like [salmon](http://salmon.readthedocs.io/en/latest/) for such cases.
  
  
  ---
  
  