How can I extract base on specific loci from all reads overlapping that?

Question

[Jeason Rad]

Hi, all. I'm using samtools to extract reads from BAM which I intrested. For example, if I want to extract all reads overlapping chr1_10000 and I know samtools view in.bam chr1:10000-10000 can do that. But how can I know which base on chr1_10000 in each reads? (especially when there has INDEL in reads)

Answer 1

[Pierre Lindenbaum]

using sam2tsv : http://lindenb.github.io/jvarkit/Sam2Tsv.html

$ java -jar dist/sam2tsv.jar ~/scf7180000354095.bam -R ~/scf_7180000354095.fasta 

#READ_NAME                              FLAG  CHROM             READ_POS  BASE  QUAL  REF_POS  REF  OP
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  0         A     D     3        A    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  1         C     B     4        C    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  2         A     D     5        A    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  3         T     D     6        T    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  4         T     D     7        T    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  5         G     E     8        G    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  6         A     H     9        A    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  7         A     H     10       A    M
HISEQ:348:H2YWCBCXX:1:1115:19355:88769  163   scf7180000354095  8         T     H     11       T    M

Answer 2

[cpad0112]

For four bases: ATG and C, gives counts base wise:

sh /<path to="" igvtools="">/igvtools count -w 1 --bases <filename>.bam <filename>.wig hg19