Question

Get genotype from reference fasta

1

Entering edit mode

9.9 years ago

win ▴ 990

Hi all,

In NGS analysis we use a reference fasta and I wanted to know if it would be possible to get genotype at a certain genomic location from that file?

fasta • 2.3k views

ADD COMMENT • link updated 2.7 years ago by Ram 44k • written 9.9 years ago by win ▴ 990

Ram · Answer 1 · 2015-01-15

3

Entering edit mode

9.9 years ago

Devon Ryan 104k

samtools faidx foo.fa chromosome:pos1-pos2

ADD COMMENT • link updated 2.7 years ago by Ram 44k • written 9.9 years ago by Devon Ryan 104k

0

Entering edit mode

And if you need to do it in a script or use Windows

pip install --user pyfaidx 
faidx foo.fa chr:pos1-pos2

ADD REPLY • link updated 2.7 years ago by Ram 44k • written 9.9 years ago by Matt Shirley 10k

0

Entering edit mode

This command produced a binary sai file. So I am still unsure how I could get genotype at pos1-pos2?

ADD REPLY • link updated 2.7 years ago by Ram 44k • written 9.9 years ago by win ▴ 990

1

Entering edit mode

It prints the sequence to the screen. It just has to build an index the first time you use it.

ADD REPLY • link 9.9 years ago by Devon Ryan 104k

1

Entering edit mode

I had to run samtools faidx file.fa and then your command. Now I can see the output in FASTA format.

ADD REPLY • link 9.9 years ago by win ▴ 990

1

Entering edit mode

Yes, an unfortunate behavior of samtools is that you must invoke faidx once to generate the index, then again to fetch sequence. This is something I did not emulate with pyfaidx.

ADD REPLY • link updated 2.7 years ago by Ram 44k • written 9.9 years ago by Matt Shirley 10k

0

Entering edit mode

This must have already gotten fixed, since the version of samtools I have (granted, I track the github repo) will create the index (if needed) and then load it in a single go.