6

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?

stack_learner
  • 1,262
  • 14
  • 26

1 Answers1

7

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
benn
  • 3,571
  • 9
  • 28
  • So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only. – stack_learner Dec 14 '17 at 10:39
  • Please show your output to these people. – benn Dec 14 '17 at 10:42
  • But you mean the bam file I have doesn't have any chr, read count and other details? I was thinking the above posted one are from header and there is another command to get the readcounts and other details. – stack_learner Dec 14 '17 at 10:46
  • It probably has read counts, but not per chromosome but per @SQ sequence (which are RNA sequences). With samtools idxstats you'll see the read counts per @SQ. – benn Dec 14 '17 at 10:49
  • Ok. I want to use featurecounts in R to get the read counts. I have the bam files and hg19.fasta and hg19.fasta.fai files. May I know how I should give the command to get the counts for all the samples. Thanq – stack_learner Dec 14 '17 at 11:00
  • 4
    Maybe you misunderstood, but you need to correctly align to hg19 first in order to do this in featureCounts. From your output it is clear that it is not aligned to hg19, but to a bunch of RNA sequences. Please discuss with your people that your data is not properly aligned to hg19 (show your output, this is proof), or try to align it yourself from fastq. When you have done that, you'll need a proper gtf annotation file from UCSC or ENSEMBL which is necessary for featureCounts. Good luck! – benn Dec 14 '17 at 11:05
  • Your options are to either align to the genome (as @b.nota said) or filter using samtools however you like (e.g., for MAPQ) and then get counts with samtools idxstats. – Devon Ryan Dec 14 '17 at 11:23
  • 3
    When you are told they are aligned to hg19, I guess they mean tat are aligned to the hg19 version of the transcriptome, not hg19 itself. In this case you done need featureCounts to get read counts (and indeed, it would be hard/impossible to do so). samtools idxstats should give you want you are looking for. – Ian Sudbery Dec 14 '17 at 12:53
  • Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats. – stack_learner Dec 14 '17 at 13:14
  • One question. Do you I need to sort the bam files? samtools idxstats what are the options I should give in the command to get read counts? – stack_learner Dec 14 '17 at 13:22
  • Yes you'll need to sort bam files first. samtools idxstats aln.sorted.bam will give you the results, see here. – benn Dec 14 '17 at 13:25
  • Yes, I saw the manual and tried like this and endup with an error.
    samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam. This gave me [W::bam_hdr_read] EOF marker is absent. The input is probably truncated; [E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes; [bam_sort_core] truncated file. Aborting.     – stack_learner Dec 14 '17 at 13:36
  • The tmp dir has 6 input.sorted.nnnn.bam files. – stack_learner Dec 14 '17 at 13:37
  • Please start a new question for that error. – benn Dec 14 '17 at 13:45