text: user1377: I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:
  samtools view -H IonXpress_input.bam
  
  Output:
  @HD VN:1.4  SO:coordinate   GO:none
  @SQ SN:NR_039978    LN:984
  @SQ SN:NM_130786    LN:1766
  @SQ SN:NM_138932    LN:9293
  @SQ SN:NM_033110    LN:2678
  @SQ SN:NM_000014    LN:4678
  @SQ SN:NM_144670    LN:5192
  @SQ SN:NR_040112    LN:1201
  @SQ SN:NM_017436    LN:2106
  @SQ SN:NM_016161    LN:1771
  @SQ SN:NM_015665    LN:1854
  
  I guess these are from header. After these I see something like following:
  @RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4
  
  Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?
  
  user939: Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case).
  For use in featureCounts my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).
  If your data is properly aligned to hg19 your output should be something like this:
  samtools view -H aligned2hg19.bam
  
  @SQ SN:chrM LN:16571
  @SQ SN:chr1 LN:249250621
  @SQ SN:chr2 LN:243199373
  @SQ SN:chr3 LN:198022430
  @SQ SN:chr4 LN:191154276
  @SQ SN:chr5 LN:180915260
  @SQ SN:chr6 LN:171115067
  @SQ SN:chr7 LN:159138663
  @SQ SN:chr8 LN:146364022
  @SQ SN:chr9 LN:141213431
  @SQ SN:chr10    LN:135534747
  @SQ SN:chr11    LN:135006516
  @SQ SN:chr12    LN:133851895
  @SQ SN:chr13    LN:115169878
  @SQ SN:chr14    LN:107349540
  @SQ SN:chr15    LN:102531392
  @SQ SN:chr16    LN:90354753
  @SQ SN:chr17    LN:81195210
  @SQ SN:chr18    LN:78077248
  @SQ SN:chr19    LN:59128983
  @SQ SN:chr20    LN:63025520
  @SQ SN:chr21    LN:48129895
  @SQ SN:chr22    LN:51304566
  @SQ SN:chrX LN:155270560
  @SQ SN:chrY LN:59373566
  
  
original_text: 1377: <p>I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:</p>
  
  <pre><code>samtools view -H IonXpress_input.bam
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>@HD VN:1.4  SO:coordinate   GO:none
  @SQ SN:NR_039978    LN:984
  @SQ SN:NM_130786    LN:1766
  @SQ SN:NM_138932    LN:9293
  @SQ SN:NM_033110    LN:2678
  @SQ SN:NM_000014    LN:4678
  @SQ SN:NM_144670    LN:5192
  @SQ SN:NR_040112    LN:1201
  @SQ SN:NM_017436    LN:2106
  @SQ SN:NM_016161    LN:1771
  @SQ SN:NM_015665    LN:1854
  </code></pre>
  
  <p>I guess these are from header. After these I see something like following:</p>
  
  <pre><code>@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4
  </code></pre>
  
  <p>Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?</p>
   939: <p>Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your <code>samtools view -H</code> code. With <code>samtools idxstats file.bam</code> you'll get the reads per chromosome (or per nucleotide sequence in your case).</p>
  
  <p>For use in <code>featureCounts</code> my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).</p>
  
  <p>If your data is properly aligned to hg19 your output should be something like this:</p>
  
  <pre><code>samtools view -H aligned2hg19.bam
  
  @SQ SN:chrM LN:16571
  @SQ SN:chr1 LN:249250621
  @SQ SN:chr2 LN:243199373
  @SQ SN:chr3 LN:198022430
  @SQ SN:chr4 LN:191154276
  @SQ SN:chr5 LN:180915260
  @SQ SN:chr6 LN:171115067
  @SQ SN:chr7 LN:159138663
  @SQ SN:chr8 LN:146364022
  @SQ SN:chr9 LN:141213431
  @SQ SN:chr10    LN:135534747
  @SQ SN:chr11    LN:135006516
  @SQ SN:chr12    LN:133851895
  @SQ SN:chr13    LN:115169878
  @SQ SN:chr14    LN:107349540
  @SQ SN:chr15    LN:102531392
  @SQ SN:chr16    LN:90354753
  @SQ SN:chr17    LN:81195210
  @SQ SN:chr18    LN:78077248
  @SQ SN:chr19    LN:59128983
  @SQ SN:chr20    LN:63025520
  @SQ SN:chr21    LN:48129895
  @SQ SN:chr22    LN:51304566
  @SQ SN:chrX LN:155270560
  @SQ SN:chrY LN:59373566
  </code></pre>
  