This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How can I extract base and read name based on position from bam file?

Hello.

Sadly, I can't not use English well....

I want to make primer for many cultivar of my subject. In my primer position, there are two bialleic SNP.

ex) .....A/C..........C/G........

I think reason of this problem is that my subject is polyploid. I already checked that there are two bases on that posistion using bam file, IGV and bam-readcount.

What I want to know is base group(?).

ex)

A..........G (A,G)

C..........C (C,C)

or

A..........C (A,C)

C..........G (C,G)

If base group of reference is (A,G) and base groups that checked are (A,G) , (C,C), I can use primer that I maked using reference sequence.

I think it is similar work of haplotype phasing. Of course I'm not sure. I can't check information like this using bam-readcount. I can check it using IGV, but it takes too long. because I want check many position.

Is there any work that I check base on my target position and read name at the same time? Or is there a better way?

Thank you.

bam sam genotype
• 2,524 views • link
updated by fravadona • fravadona • written by gh • gh

Thank you! Unfortunately, I searched information but I don't know how can I see specific region that I want using sam2tsv. T-T

• 0 views • link
written by gh • gh

I'm trying to do this using this command line.

Of course It was unsuccessful yet...

samtools view /my_bam_file.bam scaffold158:794058-794078 | awk '{print $1"\t"$3"\t"$4"\t"$6"\t"$10}'

result : read_name contig_name read_start_position CIAGAR sequence

• 0 views • link
written by gh • gh

Ummm.. pair of base on the position that I want to see in per read. I try to explain it.

enter image description here

• 0 views • link
written by gh • gh

1 answer

This is an old thread but I think I understand what the problem is: for each read overlapping two SNP positions, you want to know which base a read has at each position, so you can determine which allele combinations ("base groups") co-occur on the same read.

I wrote a small CIGAR-aware Python utility called xumi (https://github.com/Fravadona/xumi) for region extraction that can be of help here.

Given two positions pos1 and pos2 on chr1:

xumi -a -r chr:pos1-pos1,chr:pos2-pos2 input.bam

The TSV output of the above command will yield one row per read with one column per position:

#qname    chr:pos1-pos1    chr:pos2-pos2
read1     A           G
read2     C           C
read3     A           C
read4     C           G

Extending this to as many SNP positions as you need comes down to add an other "region" to the -r option.

Note: reads must cover at least one of the specified "region" to appear in the output; reads covering only one position will have an empty values in the other columns.

• 0 views • link
written by fravadona • fravadona

Log in to answer this question.

More posts like this